## New family object for the normal distribution,
## can be used by function opt_bfit() and sam_GMCMC().
normal_bamlss <- function(...) {
f <- list(
"family" = "normal",
"names" = c("mu", "sigma"),
"links" = c("identity", "log"),
"d" = function(y, par, log = FALSE) {
dnorm(y, mean = par$mu, sd = par$sigma, log = log)
},
"score" = list(
"mu" = function(y, par, ...) {
drop((y - par$mu) / (par$sigma^2))
},
"sigma" = function(y, par, ...) {
drop(-1 + (y - par$mu)^2 / (par$sigma^2))
}
),
"hess" = list(
"mu" = function(y, par, ...) {
drop(1 / (par$sigma^2))
},
"sigma" = function(y, par, ...) {
rep(2, length(y))
}
)
)
class(f) <- "family.bamlss"
return(f)
}
if (FALSE) ## Test on simulated data.
d <- GAMart()
b <- bamlss(num ~ s(x1) + s(x2) + s(x3),
data = d, family = "normal")
plot(b)
## Compute the log-likelihood using the family object.
f <- family(b)
sum(f$d(y = d$num, par = f$map2par(fitted(b)), log = TRUE))
## For using JAGS() more details are needed.
norm4JAGS_bamlss <- function(...) {
f <- normal_bamlss()
f$bugs <- list(
"dist" = "dnorm",
"eta" = BUGSeta,
"model" = BUGSmodel,
"reparam" = c(sigma = "1 / sqrt(sigma)")
)
return(f)
}
## Now with opt_bfit() and sam_JAGS().
b <- bamlss(num ~ s(x1) + s(x2) + s(x3), data = d,
optimizer = opt_bfit, sampler = sam_JAGS, family = "norm4JAGS")
plot(b)
Run the code above in your browser using DataLab