BaPreStoPro (version 0.1)

estimate,mixedRegression-method: Estimation for the hierarchical (mixed) regression model

Description

Bayesian estimation of the parameter of the hierarchical regression model $y_{ij} = f(\phi_j, t_{ij}) + \epsilon_{ij}, \phi_j\sim N(\mu, \Omega), \epsilon_{ij}\sim N(0,\gamma^2\widetilde{s}(t_{ij}))$.

Usage

"estimate"(model.class, t, data, nMCMC, propSd, adapt = TRUE, proposal = c("normal", "lognormal"))

Arguments

model.class

class of the hierarchical regression model including all required information, see mixedRegression-class

t

list or vector of time points

data

list or matrix of observation variables

nMCMC

length of Markov chain

propSd

vector of proposal variances for $\phi$

adapt

if TRUE (default), proposal variance is adapted

proposal

proposal density: "normal" (default) or "lognormal" (for positive parameters)

References

Hermann, S., K. Ickstadt, and C. H. Mueller (2016). Bayesian Prediction of Crack Growth Based on a Hierarchical Diffusion Model. Applied Stochastic Models in Business and Industry, DOI: 10.1002/asmb.2175.

Examples

mu <- c(10, 5); Omega <- c(0.9, 0.01)
phi <- cbind(rnorm(21, mu[1], sqrt(Omega[1])), rnorm(21, mu[2], sqrt(Omega[2])))
model <- set.to.class("mixedRegression",
                 parameter = list(phi = phi, mu = mu, Omega = Omega, gamma2 = 0.1),
                 fun = function(phi, t) phi[1]*t + phi[2], sT.fun = function(t) 1)
t <- seq(0, 1, by = 0.01)
data <- simulate(model, t = t, plot.series = FALSE)
est <- estimate(model, t, data[1:20,], 1000)
plot(est)