Nothing
ttest.Gibbs.AR = function(before,after,iterations=1000,areaNull=c(-.2,.2),treat=NULL,r.scale=1,alphaTheta=1,betaTheta=5,sdMet=.3, progress=TRUE,return.chains=FALSE)
{
y = c(before,after)
N = length(y)
if(is.null(treat))
{
treat = c(rep(-0.5,length(before)),rep(0.5,length(after)))
}else{
if(length(treat)!=length(y))
{
stop("Invalid condition vector: treat.")
}
}
iterations = as.integer(iterations)
if(progress){
progress = round(iterations/100)
pb = txtProgressBar(min = 0, max = as.integer(iterations), style = 3)
}else{
pb=NULL
}
pbFun = function(samps){ if(progress) setTxtProgressBar(pb, samps)}
miss = as.integer(which(is.na(y)) - 1)
Nmiss = as.integer(length(miss))
yimp = y
yimp[miss+1] = mean(y, na.rm =TRUE)
chains = .Call("RgibbsTwoSampleAR", as.numeric(yimp), N, treat, r.scale, alphaTheta, betaTheta,
areaNull[1], areaNull[2], iterations, sdMet, miss, Nmiss, progress,
pbFun, new.env(), package="BayesSingleSub")
if(progress) close(pb)
dim(chains) = c(iterations,7)
chains = data.frame(chains)
colnames(chains) = c("mu0","beta","ldens","sig2","g","rho","nullArea")
chains$delta = chains$beta/sqrt(chains$sig2)
logdens = logMeanExpLogs(chains$ldens)
nulllogdens = dcauchy(0,log=TRUE) - log(r.scale)
logbf = logdens - nulllogdens
logbfArea = log(mean(chains$nullArea)) - log(1-mean(chains$nullArea)) -
(log(diff(pcauchy(areaNull,scale=r.scale))) - log(1 - diff(pcauchy(areaNull,scale=r.scale))))
acc = mean(diff(chains$rho)!=0)
cat("\n","rho acceptance rate:",acc,"\n")
if(return.chains)
{
return(list(logbf=logbf,chains=mcmc(chains[,c(1,8,4,5,6)]),acc=acc,logbfArea=logbfArea))
}else{
return(logbf)
}
}
ttest.MCGQ.AR = function(before,after,iterations=1000,treat=NULL,method="MC",r.scale=1,alphaTheta=1,betaTheta=5)
{
y = c(before,after)
miss = as.integer(which(is.na(y)))
Nmiss = as.integer(length(miss))
if(Nmiss > 0){
y = y[-miss]
}
N = length(y)
distMat = abs(outer(1:(N+Nmiss),1:(N+Nmiss),'-'))
if(Nmiss > 0){
distMat = distMat[-miss,-miss]
}
oneVec = matrix(1,nrow=N)
Jn = matrix(1,N,N)
if(is.null(treat))
{
treat = c(rep(-0.5,length(before)),rep(0.5,length(after)))
}else{
if(length(treat)!= (length(y) + Nmiss))
{
stop("Invalid condition vector: treat.")
}
}
if(Nmiss > 0){
treat = treat[-miss]
}
iterations = as.integer(iterations)
vmlike0 = Vectorize(mlike.null.AR,"theta")
mlike0.gq = log(integrate(vmlike0,lower=0,upper=1,y=as.numeric(y),N=N,tr=treat,oneVec=oneVec,Jn = Jn, alphaTheta=alphaTheta, betaTheta=betaTheta, distMat=distMat)[[1]])
if(method=="MC")
{
logbf = .Call("RMCTwoSampleAR", as.numeric(y), N, as.integer(distMat), treat, r.scale, alphaTheta, betaTheta, iterations, package="BayesSingleSub") - mlike0.gq
}else if(method=="GQ")
{
vmlike.alt.gtheta.AR = Vectorize(mlike.alt.gtheta.AR,"theta")
vmlike.alt.g.AR = Vectorize(mlike.alt.g.AR,"g")
mlike1.gq = log(integrate(vmlike.alt.g.AR,lower=0,upper=Inf,y=y,N=N,tr=treat,oneVec=oneVec,Jn = Jn,alphaTheta=alphaTheta,betaTheta=betaTheta,fun1=vmlike.alt.gtheta.AR,distMat=distMat)[[1]])
logbf = mlike1.gq - mlike0.gq
}else{
stop(paste("Invalid method specified:",method))
}
return(-logbf)
}
mlike.null.AR = function(theta,y,tr,N=length(y),oneVec=matrix(y*0+1,ncol=1),Jn = oneVec%*%t(oneVec),alphaTheta=1,betaTheta=5,...)
{
psifunc = function(theta,distMat) theta^distMat/(1-theta^2)
invpsi = solve(psifunc(theta,...))
invpsi0 = invpsi - invpsi %*% Jn %*% invpsi/as.vector(t(oneVec)%*%invpsi%*%oneVec)
ret = - 0.5 * log(t(oneVec)%*%invpsi%*%oneVec) +
-(N-1)/2 * log(t(y)%*%invpsi0%*%y) + 0.5 * determinant(invpsi,logarithm=TRUE)$modulus + dbeta(theta,alphaTheta,betaTheta,log=TRUE)
exp(ret)
}
mlike.alt.gtheta.AR = function(theta,g,y,tr,N=length(y),oneVec=matrix(y*0+1,ncol=1),Jn = oneVec%*%t(oneVec),alphaTheta=1,betaTheta=5,rscale=1,...)
{
psifunc = function(theta,distMat) theta^distMat/(1-theta^2)
invpsi = solve(psifunc(theta,...))
opsio = as.vector(t(oneVec)%*%invpsi%*%oneVec)
invpsi0 = invpsi - invpsi%*%Jn%*%invpsi/opsio
tpsi0t = as.vector(t(tr)%*%invpsi0%*%tr) + 1/g
invpsi1 = invpsi0 - invpsi0%*%tr%*%t(tr)%*%invpsi0/tpsi0t
devs = t(y)%*%invpsi1%*%y
exp(
-((N-1)/2)*log(devs) +
(-0.5)*log(opsio) + (-0.5)*log(tpsi0t) +
0.5*determinant(invpsi,logarithm=TRUE)$modulus +
-0.5*log(g) + dbeta(theta,alphaTheta,betaTheta,log=TRUE) + log(dinvgamma(g,0.5,rscale^2/2))
)
}
mlike.alt.g.AR = function(g,y,tr,N=length(y),psifunc,oneVec=matrix(y*0+1,ncol=1),Jn = oneVec%*%t(oneVec),alphaTheta=1,betaTheta=5,rscale=1,fun1,...)
{
integrate(fun1,lower=0,upper=1,g=g,y=y,N=N,tr=tr,oneVec=oneVec,Jn = oneVec%*%t(oneVec),rscale=rscale,alphaTheta=alphaTheta,betaTheta=betaTheta,...)[[1]]
}
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.