Nothing
mvsn.mcmc <- function(Y, prior.Mu0=NULL, prior.Sigma0=NULL,
prior.muDelta0=NULL, prior.sigmaDelta0=NULL,
prior.H0=NULL, prior.P0=NULL,
nmcmc=10000, nburn=nmcmc/10, nthin=1, seed=100){
if(!is.matrix(Y) || ncol(Y) < 2 || nrow(Y) < 3){
stop("The input observations of the 'mvsn.mcmc' has to be a matrix
of more than two rows and one column.")
}
if(!is.null(prior.Mu0) && (!is.vector(prior.Mu0) || length(prior.Mu0)!=ncol(Y))){
stop("The mean vector of multivariate normal prior of the
parameter 'mu' of mvsn has to be a vector of the same
length as the number of columns of observations.")
}
if(!is.null(prior.Sigma0) && (!is.matrix(prior.Sigma0) || any(eigen(prior.Sigma0)$values<=0) || any(dim(prior.Sigma0)!=ncol(Y)))){
stop("The variance matrix of multivariate normal
prior of the parameter 'mu' of mvsn has to be a positive definite
matrix with the number of rows and columns equal to
the number of columns of observations.")
}
if(!is.null(prior.muDelta0) && (!is.vector(prior.muDelta0) || length(prior.sigmaDelta0)!=ncol(Y))){
stop("The vector of the means of normal prior
of parameter 'D' of mvsn has to be a vector of the same
length as the number of columns of observations.")
}
if(!is.null(prior.sigmaDelta0) && (!is.vector(prior.sigmaDelta0) || length(prior.sigmaDelta0)!=ncol(Y) || any(prior.sigmaDelta0 <= 0))){
stop("The vector of the standard deviations of normal prior
of parameter 'D' of mvsn has to be a vector of the same
length as the number of columns of observations and
its components have to be positive.")
}
if(!is.null(prior.H0) && (!is.matrix(prior.H0) || any(eigen(prior.H0)$values<=0) || any(dim(prior.H0)!=ncol(Y)))){
stop("The inverse of scale matrix of Wishart prior of the inverse of
parameter 'Sigma' of mvsn has to be a positive definite
matrix with the number of rows and columns equal to
the number of columns of observations.")
}
if(!is.null(prior.P0) && (!is.vector(prior.P0) || length(prior.H0)!=1 || prior.H0 < ncol(Y))){
stop("The degrees of freedom of Wishart prior of the inverse of
parameter 'Sigma' of mvsn has to be positive scalar and no less than
the number of columns of observations.")
}
P <- ncol(Y)
N <- nrow(Y)
if(is.null(prior.Mu0)){
Mu0 <- apply(Y, 2, median)
}
else{
Mu0 <- prior.Mu0
}
if(is.null(prior.Sigma0)){
prior.Sigma0 <- 100*var(Y)
}
Tau0 <- solve(prior.Sigma0)
if(is.null(prior.muDelta0)){
prior.muDelta0 <- rep(0, P)
}
muDelta0 <- prior.muDelta0
if(is.null(prior.sigmaDelta0)){
prior.sigmaDelta0 <- sqrt(diag(var(Y)))
}
tauDelta0 <- 1/(prior.sigmaDelta0^2)
if(is.null(prior.P0)){
P0 <- P
}
else{
P0 <- prior.P0
}
if(is.null(prior.H0)){
H0 <- P*var(Y)
}
else{
H0 <- prior.H0
}
Delta.ini <- mvrnorm(1, rep(0, P), diag(diag(var(Y)/(1-2/pi))))
data=list("Y", "N", "P", "muDelta0", "tauDelta0", "Mu0", "Tau0", "H0", "P0")
inits = function(){list(Delta=Delta.ini)}
cpath <- getwd()
setwd(cpath)
cat(
"model {
for(i in 1:N){
Y[i,1:P] ~ dmnorm (MuY[i, 1:P], TauY[1:P, 1:P])
MuY[i, 1:P] <- Mu + Dz[i, 1:P]
}
for(i in 1:N){
for(j in 1:P){
Dz[i,j] <- Delta[j]*Z[i,j]
}
}
for(j in 1:P){
Delta[j] ~ dnorm(muDelta0[j], tauDelta0[j])
}
for(i in 1:N){
for(j in 1:P){
Z[i, j] ~ dnorm(0, 1)T(0, )
}
}
Mu[1:P] ~ dmnorm(Mu0[1:P], Tau0[1:P, 1:P])
TauY[1:P , 1:P] ~ dwish(H0[1:P, 1:P], P0)
Sigma[1:P, 1:P] <- inverse(TauY)
}",
file="mvsn.jags")
set.seed(seed)
jagsfit <- jags(data, inits,
parameters.to.save=c("Mu", "Delta", "Sigma"),
model.file="mvsn.jags",
n.chains = 1, n.iter=nmcmc, n.burnin=nburn, n.thin=nthin)
unlink("mvsn.jags")
Mu <- jagsfit$BUGSoutput$sims.list$Mu
Sigma <- jagsfit$BUGSoutput$sims.list$Sigma
Delta <- jagsfit$BUGSoutput$sims.list$Delta
DIC=jagsfit$BUGSoutput$DIC
list(Mu=Mu, Delta=Delta, Sigma=Sigma, DIC=DIC)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.