Nothing
scamH.formula <-
function(formula,data,start,priors,rotationVal,rotationVect,wParamX,wOuterParamX,sdvalue,accept,likelihoodtype="logit",burning=TRUE,rotation=TRUE,iter,...){
#### F. Guillaume Blanchet - Janvier 2013, March 2013, April 2013
##########################################################################################
#========================
#### Initiation of object
#========================
nparamX<-ncol(start$paramX)
nsp<-ncol(data$Y)
paramXCurrent<-start$paramX
means<-start$means
sigma<-start$sigma
#### Calculate log-likelihood - mean
likeCurrentMean<-numeric()
for(i in 1:nsp){
dataLL<-data.frame(data$Y[,i],data$X)
colnames(dataLL)[1]<-attributes(data)$Ypattern
likeCurrentMean[i]<-loglikelihood.formula(formula,paramXCurrent[i,],dataLL,likelihoodtype=likelihoodtype)
}
#### Calculate log-likelihood - sigma
likeCurrentSigma<-numeric()
for(i in 1:nsp){
likeCurrentSigma[i]<-loglikelihood(sigma,paramXCurrent[i,],means,likelihoodtype="sigma")
}
#=============
#### Algorithm
#=============
#----------------------------------------
### Update means (each one independently)
#----------------------------------------
paramXNew<-matrix(NA,nrow=nsp,ncol=nparamX)
colnames(paramXNew)<-colnames(paramXCurrent)
for(k in 1:nsp){
for(j in 1:nparamX){
### Sample from a normal distribution and weight the sample by the eigenvalue of the covariance matrix
paramXNew[k,]<-paramXCurrent[k,]+rnorm(1,mean=0,sd=sdvalue[k,j])*sqrt(rotationVal[k,j])*rotationVect[,j,k]
#### Calculate log-likelihood - mean
dataLL<-data.frame(data$Y[,k],data$X)
colnames(dataLL)[1]<-attributes(data)$Ypattern
likeNewMean<-loglikelihood.formula(formula,paramXNew[k,],dataLL,likelihoodtype=likelihoodtype)
#### Calculate log-likelihood - sigma
likeNewSigma <- loglikelihood(sigma,paramXNew[k,],means,likelihoodtype="sigma")
accept[k,j,2]<-accept[k,j,2]+1
### Reject if the likelihood is -Inf
if(!(is.infinite(likeNewMean) & sign(likeNewMean)==-1)){
### Accept or reject the new paramX
if(runif(1) <= exp(likeNewMean+likeNewSigma-likeCurrentMean[k]-likeCurrentSigma[k])){
paramXCurrent[k,]<-paramXNew[k,]
likeCurrentMean[k]<-likeNewMean
likeCurrentSigma[k]<-likeNewSigma
accept[k,j,1]<-accept[k,j,1]+1
}
}
}
}
#--------------------------
#### means and sigma update
#--------------------------
meansparamX<-colMeans(paramXCurrent)
SigmaparamX<-crossprod(sweep(paramXCurrent,2,meansparamX,FUN="-"))
MeansSp<-(priors$kappa0/(priors$kappa0+nsp))*priors$means0+(nsp/(priors$kappa0+nsp))*meansparamX
kappaSp<-priors$kappa0+nsp
nuSp<-priors$nu0+nsp
LambdaSp<-priors$Lambda0 + SigmaparamX + (priors$kappa0*nsp)/(priors$kappa0+nsp)*outer(meansparamX-priors$means0,meansparamX-priors$means0)
iLambdaSp<-solve(LambdaSp)
isigma <- rWishart(1,nuSp,iLambdaSp)[,,1]
sigma <- solve(isigma)
sigmaKappa<-sigma/kappaSp
means<-mvrnorm(1,MeansSp,sigmaKappa)
if(burning){
wAccept<-1-0.1*exp(-iter/500) # w
### Greediness during the adaptation period
greedAdapt<-1+exp(-iter/500) # q
### Calculate the accept ratio
acceptRatio<-accept[,,1]/accept[,,2]
### Modify the standard deviation
sdvalue<-sdvalue*greedAdapt^(acceptRatio-0.44)
### Weight the accept ratio
accept<-accept*wAccept
### Weight the current ParamX
wParamX<-wParamX+wAccept*paramXCurrent # s1
### Weight the outer product of the current ParamX
for(j in 1:nsp){
wOuterParamX[,,j]<-wOuterParamX[,,j]+wAccept*outer(paramXCurrent[j,],paramXCurrent[j,]) # s2
}
if(rotation){
### sum the weight
weight<-sum(1-0.1*exp(-(1:iter)/500))
for(j in 1:nsp){
### Construct a weighted covariance matrix
covMat<-(wOuterParamX[,,j]-outer(wParamX[j,],wParamX[j,])/weight)/(weight-1)+diag(sqrt(.Machine$double.eps),nparamX)
rotationBase<-eigen(covMat)
rotationVal[j,]<-rotationBase$values
rotationVect[,,j]<-rotationBase$vectors
}
}
}
#### Result object
colnames(sigma)<-names(means)
rownames(sigma)<-names(means)
startNew<-list(paramX=paramXCurrent,sigma=sigma,means=means)
rotNew<-list(values=rotationVal,vectors=rotationVect,wParamX=wParamX,wOuterParamX=wOuterParamX)
res<-list(start=startNew,rotation=rotNew,sdvalue=sdvalue,accept=accept)
return(res)
}
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.