## main sampler:
IsingOccu_Murray_sampler_Ising_det <- function(X,detmat,detX # list with each entry as a period
,mcmc.iter = 10000, burn.in = 10
,vars_prop = list( beta_occu = rep(1e-5,2 * ncol(X))
,beta_det = rep(1e-5,2 * (ncol(detX[[1]][[1]]) + ncol(X)) )
,eta_intra = rep(1e-5,nspp)
,eta_inter = rep(1e-5,nspp*(nspp-1)/2)
,d_intra=rep(1e-5,nspp)
,d_inter = rep(1e-5,nspp)
,spp_mat = 1e-5
,spp_mat_det = 1e-5)
,para_prior = list( beta_occu = rep(1000,2 * ncol(X))
,beta_det = rep(1000,2 * (ncol(detX[[1]][[1]]) + ncol(X)) )
,eta_intra = rep(1e-1,nspp)
,eta_inter = rep(1000,nspp*(nspp-1)/2)
,d_intra=rep(1000,nspp)
,d_inter = rep(1000,nspp)
,spp_mat = 1000
,spp_mat_det = 1000)
,uni_prior = T
,Zprop_rate = 1
,distM,link_map
,dist_mainland , link_mainland
,int_range_intra="nn",int_range_inter="nn"
,seed = 42,ini,thin.by = 100,report.by=100,nIter=100, Gibbs = T,method="CFTP"){ # ini has same formate of theta
cat("Initializing...\n\n")
require(coda)
require(Matrix)
require(RcppArmadillo)
source("./R/misc.R")
source("./R/Murray_ratio.R")
source('./R/importance_Z_helper.R')
Rcpp::sourceCpp("./src/IsingCpp_CFTP_sparse.cpp")
set.seed(seed)
if(uni_prior) getlogprior <- getlogprior_uniform
else getlogprior <- getlogprior_normal
vars_prop <- vars_prop[names(ini)]
para_prior <- para_prior[names(ini)]
cat("Setting for imperfect observation and missing sites:\n")
if(Zprop_rate==0) cat(" Perfect observation, given by Z\n")
else cat(" Imperfect observation, Z propose with rate",Zprop_rate,"\n")
cat("MCMC reported every",report.by,"iterations and thinned by",thin.by,"\n\n")
nsite <- (nrow(distM))
theta <- ini
spp_neig <- 1 *( spp_mat!=0 )
nspp <- nrow(spp_mat)
nrep <- length(detmat)
#theta_tuta=ini
ini <- lapply(ini,as.matrix)
theta.mcmc <- list()
for(i in 1:length(ini)){
theta.mcmc[[i]] <- mcmc(matrix(nrow = floor(mcmc.iter/thin.by),ncol = length(ini[[i]])),thin = thin.by)
}
names(theta.mcmc) <- names(ini)
Z.mcmc <- mcmc(matrix(nrow = floor(mcmc.iter/thin.by),ncol = nsite*nspp*nrep),thin = thin.by)
#Z_absolute = (sapply(detmat,rowSums)>0) * 2 - 1
detmat_nona <- lapply(detmat,function(mat){
mat[is.na(mat)] <- -1
return(mat)
})
Z_absolute <- (sapply(detmat_nona,function(detmat_i){rowSums((detmat_i+1)/2)>0})) * 2 - 1
rm(detmat_nona)
Z_curr <- Z_absolute
Z_temp <- Z_absolute
MRF_curr <- getMRF(ini,envX,distM,link_map,dist_mainland,link_mainland,int_range_intra,int_range_inter)
theta_curr <- ini
cat("Algorithm start...\n")
accept_Z <- 0
accept_Z_missing_obs <- 0
accept_theta_occu <- 0
accept_theta_det <- 0
low_acc_Z <- 0
low_acc_Z_missing_obs <- 0
low_acc_theta_occu <- 0
low_acc_theta_det <- 0
propose_Z_num <- 0
propose_Z_missing_obs <- 0
timing <- proc.time()
n_para_group <- length(theta_curr)
constrains <- as.vector( Z_absolute)
constrains[constrains==-1] <- NA
for(i in 1:(burn.in+mcmc.iter)){# to burn in
#propose theta
theta_prop <- theta_curr
for(j in c(1:length(theta_curr))[-c(2,n_para_group)]){ # no detection proposing
theta_prop[[j]] <- matrix( rnorm(length(theta_curr[[j]]),mean = 0,sd = sqrt(vars_prop[[j]])),nrow(theta_curr[[j]]),ncol(theta_curr[[j]]) )+ theta_curr[[j]]
}
theta_prop$spp_mat <- theta_prop$spp_mat * spp_neig # theta_prop$spp_mat=theta_prop$spp_mat * spp_neig
theta_prop$spp_mat <- .5*(theta_prop$spp_mat + t( theta_prop$spp_mat)) # must be sym
# MH ratio
MRF_prop <- getMRF(theta_prop,envX,distM,link_map,dist_mainland,link_mainland,int_range_intra,int_range_inter)
Z_temp <- rIsingOccu_multi(MRF_prop,n=1,method = method,nIter)
Murray_ratio <- Murray_ratio_occu_theta(MRF_curr ,MRF_prop, getlogprior(theta_prop,theta_curr,para_prior)
,Z_curr
,Z_temp
,para_prior
,distM,link_map
,dist_mainland,link_mainland
,int_range_intra,int_range_inter)
r <- runif(1)
if(is.na(Murray_ratio)){
Murray_ratio <- 0
#cat("occuNA\n")
}
if(Murray_ratio<exp(-10)) low_acc_theta_occu <- low_acc_theta_occu + 1
if(r<=Murray_ratio){
theta_curr <- theta_prop
#Z_temp = Z_temp
MRF_curr <- MRF_prop
accept_theta_occu <- accept_theta_occu + 1
}
theta_prop <- theta_curr
theta_prop[[2]] <- matrix( rnorm(length(theta_curr[[2]]),mean = 0,sd = sqrt(vars_prop[[2]])),nrow(theta_curr[[2]]),ncol(theta_curr[[2]]) )+ theta_curr[[2]]
theta_prop[[n_para_group]] <- matrix( rnorm(length(theta_curr[[n_para_group]]),mean = 0,sd = sqrt(vars_prop[[n_para_group]])),nrow(theta_curr[[n_para_group]]),ncol(theta_curr[[n_para_group]]) )+ theta_curr[[n_para_group]]
theta_prop$spp_mat_det <- theta_prop$spp_mat_det * spp_neig
theta_prop$spp_mat_det <- .5*(theta_prop$spp_mat_det + t( theta_prop$spp_mat_det)) # must be sym
Murray_ratio <- MH.ratio.Ising_det(theta_curr ,theta_prop,getlogprior(theta_prop,theta_curr,para_prior)
,Z_curr
,detmat
,vars_prior
,envX, detX
)
r <- runif(1)
if(is.na(Murray_ratio)) {
Murray_ratio <- 0
#cat("detNA\n")
}
if(Murray_ratio<exp(-10)) low_acc_theta_det <- low_acc_theta_det + 1
if(r<=Murray_ratio){
theta_curr <- theta_prop
accept_theta_det <- accept_theta_det + 1
}
Z_prop <- Z_curr
if(Gibbs){
Z_curr <- Gibbs_Z_rep(theta_curr, envX, detX,detmat,Z_curr,Z_absolute,MRF_curr)
}
else{
Z_prop <- propose_Z_plain(Z_curr,Z_absolute,Zprop_rate)
MH_ratio <- MH_ratio_Z_plain(theta_curr, MRF_curr,Z_curr, Z_prop
,detmat,envX, detX)
r <- runif(1)
if(is.na(MH_ratio)) {
MH_ratio <- 0
#cat("ZNA\n")
}
if(MH_ratio<exp(-10)) low_acc_Z <- low_acc_Z + 1
if(r<=MH_ratio){
accept_Z <- accept_Z + 1 - (sum(Z_prop==Z_curr)==length(Z_prop))
Z_curr <- Z_prop
}
}
if(i%%report.by == 0) {
if(i<=burn.in) cat("Burn in iteration",i-report.by+1,"to",i,":\n\n")
else cat("Sampling iteration",i-burn.in-report.by+1,"to",i-burn.in,"\n\n")
#cat(" # of Z proposed for imperfect detection: ",propose_Z_num,"\n")
if(!Gibbs){
cat(" # of Z acceptance: " , accept_Z,"\n")
cat(" # of Z acceptance ratio <exp(-10): ",low_acc_Z,"\n\n")
}
cat(" # of occupancy theta acceptance: " , accept_theta_occu,"\n")
cat(" # of occupancy acceptance ratio <exp(-10): ",low_acc_theta_occu,"\n\n")
cat(" # of detection theta acceptance:" , accept_theta_det,"\n")
cat(" # of detection acceptance ratio <exp(-10): ",low_acc_theta_det,"\n\n")
timing = proc.time()- timing
cat("Time used in this" ,report.by,":",timing[1],"s\n")
cat("\n\n")
accept_Z <- 0
accept_Z_missing_obs <- 0
accept_theta_occu <- 0
accept_theta_det <- 0
low_acc_Z <- 0
low_acc_Z_missing_obs <- 0
low_acc_theta_occu <- 0
low_acc_theta_det <- 0
propose_Z_num <- 0
propose_Z_missing_obs <- 0
timing <- proc.time()
}
sample_iter <- i-burn.in
## saving result when in sample
if( sample_iter > 0 & sample_iter %% thin.by==0){
for(j in 1:length(theta_curr)){
theta.mcmc[[j]][sample_iter/thin.by,] <- as.vector( theta_curr[[j]])
} # saving the results
Z.mcmc[sample_iter/thin.by,] <- Z_curr
}
}
theta.mean <- lapply(theta.mcmc,function(thetaa){ apply(thetaa,2,mean,na.rm=T)})
res <- list(theta.mcmc = theta.mcmc
,means = theta.mean
,Z.mcmc = Z.mcmc,vars=vars_prop
,distM = distM
,dist_mainland = dist_mainland
,linkmap = link_map
,link_mainland = link_mainland
,interaction.range =list( inter = int_range_inter,intra = int_range_intra)
,envX=X, detX = detX,detmat = detmat,uni_prior=uni_prior,para_prior = para_prior,method = method)
class(res) <- "IsingOccu_samples"
return(res)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.