R/loglikelihood.matrix.R

In HMSC: Hierarchical modelling of species community

loglikelihood.matrix <-
function(mat,param,means,likelihoodtype="sigma",...){
#### F. Guillaume Blanchet - Septembre 2013
##########################################################################################
	if(likelihoodtype=="sigma"){
		symmetric<-isSymmetric.matrix(mat, check.attributes = FALSE)
		if(!symmetric){
			### This correct for rounding error and makes mat symmetric
			mat <- (mat+t(mat))/2
		}
		imat<-chol2inv(chol(mat))
		paramCt<-param-means
		resFinal <- drop(-paramCt%*%imat%*%paramCt/2)
	}
	
	if(likelihoodtype=="RandomCov"){
		symmetric<-isSymmetric.matrix(mat, check.attributes = FALSE)
		if(!symmetric){
			### This correct for rounding error and makes mat symmetric
			mat <- (mat+t(mat))/2
		}
		resFinal <- -(1/2)*log(det(mat))
	}
	
	if(likelihoodtype=="R"){
		nsites<-nrow(mat)
		nsp<-ncol(mat)
		R<-matrix(NA,nsp,nsp)
		for(i in 1:nsp){
			for(j in 1:nsp){
				R[i,j]<-param[i,j]/sqrt(param[i,i]*param[j,j])
			}
		}
		
		### Écrire dans la documentation !
		iparam<-solve(param)
		prior<- -(nsp+1)*log(det(param))-sum(diag(iparam))/2
		
		### likelihood
		matinvRmat <- apply(mat, 1, function(x) x %*% solve(R) %*% x)
		likeres <- -(nsites/2)*log(det(R))-sum(matinvRmat)/2+prior
		
		resFinal<-list(likelihood=likeres,R=R)
	}
	
	if(likelihoodtype=="wishart"){
		### Wishart prior (not yet implemented and keep it as comment when it is coded to make sure it is "easily" available when needed)
#		
		### Special inverse-Wishart prior leading to marginally uniform densities for all off-diagonal values of R
		
	}
	
	return(resFinal)
}