Nothing
#'@importFrom terra rasterize
#'@importFrom stats dgamma dpois rnorm runif rgamma
#Occurrence raster
.OccRast <- function(x, ras){
rast <- terra::rasterize(x, ras, fun = "count")
return(rast)
}
#get number of decimal places
.DecimalPlaces <- function(x) {
if ((x %% 1) != 0) {
nchar(strsplit(sub('0+$', '',
as.character(x)),
".", fixed = TRUE)[[1]][[2]])
} else {
return(0)
}
}
get_lambda_ij <- function(q, w, X) {
if (is.numeric(X)) {
return(q * exp(-w * X))
} else {
return(as.vector(q * exp(-(w %*% t(X)))))
}
}
get_poi_likelihood <- function(Xcounts, lambdas) {
return(dpois(Xcounts, lambdas, log = TRUE))
}
multiplier_proposal <- function(i, d = 1.2) {
u <- runif(length(i))
l <- 2 * log(d)
m <- exp(l * (u - 0.5))
ii <- i * m
hastings_ratio <- sum(log(m))
return(list(ii, hastings_ratio))
}
get_post_rate <- function(w,alpha){
a0 <- 1
b0 <- 0.001
rate <- rgamma(1, a0 + length(w) * alpha, rate = b0 + sum(w))
return(rate)
}
.RunSampBias <- function(x,
rescale_counts = 1,
rescale_distances = 1000,
iterations = 1e+05,
burnin = 0,
post_samples = NULL,
outfile = NULL,
prior_q = c(1, 0.01),
prior_w = c(1, 1),
run_null_model = FALSE,
use_hyperprior = TRUE,
verbose = FALSE) {
indx <- c(3:ncol(x))
# X is the distance matrix with 1 row for each cell and 1 col for each predictor
X <- x[, indx] / rescale_distances
# if (length(indx) > 1) {
# X <- t(X)
# }
Xcounts <- x$record_count
# coefficient, i.e. Poisson rate at zero distance from all
qA <- mean(Xcounts)
# init weigths of predictors
wA <- abs(rnorm(length(indx), 0.01, 0.01))
if (run_null_model == TRUE) {
wA = rep(0, length(indx))
}
if (use_hyperprior & length(indx) > 1) {
# hp only used with more than one bias factor
hp_update_freq = 0.01
} else {
hp_update_freq = 0
}
lambdas <- get_lambda_ij(qA, wA, X)
likA <- sum(get_poi_likelihood(Xcounts, lambdas))
priorA <- dgamma(qA, prior_q[1], rate = prior_q[2], log = TRUE)
if (!run_null_model) {
priorA <- priorA + sum(dgamma(wA, prior_w[1],
rate = prior_w[2], log = TRUE))
}
names_b <- paste("w", names(x[, indx]), sep = "_")
out <- data.frame()
if (!is.null(outfile)) {
file_name <-
paste("mcmc", paste(names_b[indx - 2], collapse = "_"), ".log", sep = "_")
outfile <- file.path(outfile, file_name)
cat(
c("it", "likA", "priorA", "q", names_b[indx - 2], "hp_rate", "\n"),
file = outfile,
sep = "\t"
)
}
# print to screen
if (verbose) {
message(paste(c("it", "likA", "priorA", "q", names_b[indx - 2]),
collapse = " "))
}
for (it in 1:iterations) {
w <- wA
q <- qA
mh <- TRUE
if (runif(1) < hp_update_freq & run_null_model == FALSE) {
prior_w[2] <- get_post_rate(wA, prior_w[1])
mh <- FALSE
hastings <- 0
} else if (runif(1) < 0.3 | run_null_model) {
update <- multiplier_proposal(qA, d = 1.1)
q <- update[[1]]
hastings <- update[[2]]
} else {
update <- multiplier_proposal(wA, d = 1.05)
w <- update[[1]]
hastings <- update[[2]]
}
lambdas <- get_lambda_ij(q, w, X)
if (mh) {
lik <- sum(get_poi_likelihood(Xcounts, lambdas))
} else {
lik <- likA
}
prior <- dgamma(q, prior_q[1], rate = prior_q[2], log = TRUE)
if (!run_null_model) {
prior <- prior + sum(dgamma(w, prior_w[1], rate = prior_w[2],
log = TRUE))
}
if (mh == FALSE | (lik + prior) - (likA + priorA) + hastings >= log(runif(1))) {
likA <- lik
priorA <- prior
wA <- w
qA <- q
}
if (it > burnin & it %% 100 == 0) {
if (verbose) {
message(paste(round(c(it, likA, priorA, qA, wA), 3), collapse = " "))
}
out <- rbind(out, c(it, likA, priorA, qA, wA, prior_w[2]))
if (!is.null(outfile)) {
cat(
c(it, likA, priorA, qA, wA, prior_w[2], "\n"),
file = outfile,
sep = "\t",
append = TRUE
)
}
}
}
names(out) <- (c("it", "likA", "priorA", "q", names_b[indx - 2],"hp_rate"))
return(out)
}
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.