# R/utilsFrancesco.R In BayesianTools: General-Purpose MCMC and SMC Samplers and Tools for Bayesian Statistics

#### Defines functions as_mcmc.matrixcalc_mpsrfmpsrf

```as_mcmc.matrix = function(x, names, par = NULL) {
dimnames(x)=list(NULL, names)
if (is.null(par)) {
attr(x, 'mcpar') = c(1, nrow(x), 1)
} else {
attr(x, 'mcpar') = par
}
class(x) = 'mcmc'
return(x)
}

calc_mpsrf = function(x, end) {
x = window(x, start(x), end, 1)
Niter = niter(x)
Nchain = nchain(x)
Nvar = nvar(x)
x = lapply(x, as.matrix)
S2 = array(sapply(x, var, simplify=TRUE), dim=c(Nvar,Nvar,Nchain))
W = apply(S2, c(1,2), mean)
xbar = matrix(sapply(x, apply, 2, mean, simplify=TRUE), nrow=Nvar, ncol=Nchain)
B = Niter * var(t(xbar))
CW = chol(W)
emax = eigen(backsolve(CW, t(backsolve(CW, B, transpose=TRUE)), transpose=TRUE), symmetric=TRUE, only.values=TRUE)\$values
res = sqrt( (1 - 1/Niter) + (1 + 1/Nvar) * emax/Niter )
return(res)
}

mpsrf = function(x, step = 50, ...) {
if (nchain(x) < 2 | nvar(x) == 1)
stop("You need at least two chains and two parameters.")
z = seq(start(x)-1+step, end(x), by = step)
res = sapply(z, function(i) calc_mpsrf(x,i))
x = data.frame(z,MPSRF = res)
return(x)
}
```

## Try the BayesianTools package in your browser

Any scripts or data that you put into this service are public.

BayesianTools documentation built on Dec. 10, 2019, 1:08 a.m.