#' Wrapper function for Gibbs sampler
#'
#' Wrapper function for Gibbs sampler
#' @param df.sptmod model data
#' @param method method to be used
#' @param theta.init prior mean of the parameters
#' @param mixture boolean. Is the model a mixture model.
#' @param N.run Number of runs of the sampler.
#' @param beta.range Range of beta for Potts.
#' @param print.res logical. printing or not the steps values.
#' @keywords estimation
#' @export
#' @examples
#' data("wq_analysis_week2")
#' SPTMData(wq.raw.obs, frequency = "quarter")
estimGibbs = function(df.sptmod, method = 'simple', priors = NULL, mixture = TRUE, N.run = 10000, beta.range = c(0.5, 3), print.res = FALSE, debug=FALSE, meanFunc = NULL, Xcov = NULL){
if(debug){browser()}
p.comp = df.sptmod$lu
season = df.sptmod$data$season
basis = df.sptmod$data$basis
season.sites=season
season.lu = updateSeason(season.sites,p.sites= p.comp)
y = df.sptmod$data$obs.data
if(method == 'simple'){
theta.mean = rnorm(ncol(p.comp), mean = rep(mean(y$obs, na.rm=T),ncol(p.comp)),
sd = rep(mean(abs(y$obs), na.rm=T),ncol(p.comp))/10)
theta.season = rnorm(ncol(season.lu), mean = 1, sd = 0.1)
#theta.sigma = rnorm(ncol(p.comp), mean = 0.5, sd = 0.05)
if(is.null(priors)){
priors = list(theta = list(m0 = theta.init[1:(ncol(season.lu))], s0 = abs(theta.init[1:(ncol(season.lu))]) / 10),
sigma2 = list(a0 = rep(5, ncol(p.comp)), b0 = rep(2, ncol(p.comp))))
}else{priors = priors}
res = gibbSampler(priors, y, p.sites = p.comp, season.sites = season.sites,basis = basis,
ID = levels(y$ID), N.run = N.run, mixture = mixture, meanFunc = meanFunc, debug, print.res = print.res)
p.comp = apply(res$z, 1:2, mean)
RET = list(Result = list(theta = colMeans(res$theta), latent = p.comp),
seasonSites = season.sites,
GibbsOut = res,
data = df.sptmod)
}
if(method == 'spat'){
dd = deldir(df.sptmod$coord[,1], df.sptmod$coord[,2],suppressMsge=TRUE)
tile.dd = tile.list(dd)
Bds = comPts(tile.dd)
theta.mean = rnorm(ncol(p.comp), mean = rep(mean(y$obs, na.rm=T),ncol(p.comp)),
sd = rep(mean(abs(y$obs), na.rm=T),ncol(p.comp))/10)
theta.season = rnorm(ncol(season.lu), mean = 1, sd = 0.1)
theta.sigma = rnorm(ncol(p.comp), mean = 0.5, sd = 0.05)
if(is.null(priors)){
priors = list(theta = list(m0 = theta.init[1:(ncol(season.lu))], s0 = theta.init[1:(ncol(season.lu))] / 10),
sigma2 = list(a0 = rep(5, ncol(p.comp)), b0 = rep(2, ncol(p.comp))))
}else{priors = priors}
res = gibbSPTM(priors, y, p.sites = p.comp, season.sites = season.sites,basis = basis,
Coords=NULL ,Bds=Bds, ID = levels(y$ID), N.run = N.run, beta.range = beta.range,
meanFunc = meanFunc, debug, print.res = print.res)
p.comp = apply(res$z, 1:2, mean)
RET = list(Result = list(theta = colMeans(res$theta), latent = p.comp),
seasonSites = season.sites,
GibbsOut = res,
data = df.sptmod)
}
if(method == 'mixedeffect'){
dd = deldir(df.sptmod$coord[,1], df.sptmod$coord[,2],suppressMsge=TRUE)
tile.dd = tile.list(dd)
Bds = comPts(tile.dd)
theta.mean = rnorm(ncol(p.comp), mean = rep(mean(y$obs, na.rm=T),ncol(p.comp)),
sd = rep(mean(abs(y$obs), na.rm=T),ncol(p.comp))/10)
theta.season = rnorm(ncol(season.lu), mean = 1, sd = 0.1)
theta.sigma = rnorm(ncol(p.comp), mean = 0.5, sd = 0.05)
if(is.null(priors)){
priors = list(theta = list(m0 = theta.init[1:(ncol(season.lu))], s0 = theta.init[1:(ncol(season.lu))] / 10),
sigma2 = list(a0 = rep(5, ncol(p.comp)), b0 = rep(2, ncol(p.comp))))
}else{priors = priors}
attr(basis, "indexSeason") <- df.sptmod$data$list.idx
res = gibbSPTMixed(priors, y, p.sites = p.comp, season.sites = season.sites,basis = basis,
Coords=NULL ,Bds=Bds, ID = levels(y$ID), N.run = N.run, beta.range = beta.range,
meanFunc = meanFunc, debug, print.res = print.res, Xcov = df.sptmod$covariates)
p.comp = apply(res$z, 1:2, mean)
RET = list(Result = list(theta = colMeans(res$theta), latent = p.comp),
seasonSites = season.sites,
GibbsOut = res,
data = df.sptmod)
}
class(RET) <- 'sptmRes'
return(RET)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.