R/scamH.HMSCdata.R

In HMSC: Hierarchical modelling of species community

scamH.HMSCdata <-
function(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, June 2013, February 2014
##########################################################################################
	#========================
	#### Initiation of object
	#========================
	nparamX<-ncol(start$paramX)
	nsp<-ncol(data$Y)
	
	paramXCurrent<-start$paramX
	means<-start$means
	sigma<-as.matrix(start$sigma)
	
if(TRUE){
	#### Calculate log-likelihood - mean
	likeCurrentMean<-numeric()
	if(!any(class(data)=="HMSCdataRandom")){
		for(i in 1:nsp){
			likeCurrentMean[i]<-loglikelihood(data$Y[,i],data$X,paramXCurrent[i,],likelihoodtype=likelihoodtype)
		}
	}else{
		for(i in 1:nsp){
			likeCurrentMean[i]<-loglikelihood(data$Y[,i],data$X,paramXCurrent[i,],data$Random,start$paramRandom[i,],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)
	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
			if(!any(class(data)=="HMSCdataRandom")){
				likeNewMean<-loglikelihood(data$Y[,k],data$X,paramXNew[k,],likelihoodtype=likelihoodtype)
			}else{
				likeNewMean<-loglikelihood(data$Y[,k],data$X,paramXNew[k,],data$Random,start$paramRandom[k,],likelihoodtype=likelihoodtype)
			}
			#### Calculate log-likelihood - sigma
			likeNewSigma <- loglikelihood(sigma,paramXNew[k,],means,likelihoodtype="sigma")
			
			accept[k,j,2]<-accept[k,j,2]+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
			}
		}
	}
}	
if(TRUE){
	#--------------------------
	#### 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)/(kappaSp)*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(TRUE){
	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)
}