Skip to content

Instantly share code, notes, and snippets.

@florianhartig
Last active November 9, 2021 11:07
Show Gist options
  • Save florianhartig/9078258 to your computer and use it in GitHub Desktop.
Save florianhartig/9078258 to your computer and use it in GitHub Desktop.
A modification of code posted originally by Petr Keil http://www.petrkeil.com/?p=1910, to illustrate some comments I made in response to this blog post
library(mvtnorm) # to draw multivariate normal outcomes
library(R2jags) # JAGS-R interface
# function that makes distance matrix for a side*side 2D array
dist.matrix <- function(side)
{
row.coords <- rep(1:side, times=side)
col.coords <- rep(1:side, each=side)
row.col <- data.frame(row.coords, col.coords)
D <- dist(row.col, method="euclidean", diag=TRUE, upper=TRUE)
D <- as.matrix(D)
return(D)
}
# function that simulates the autocorrelated 2D array with a given side,
# and with exponential decay given by lambda
# (the mean mu is constant over the array, it equals to global.mu)
cor.surface <- function(side, global.mu, lambda)
{
D <- dist.matrix(side)
# scaling the distance matrix by the exponential decay
SIGMA <- exp(-lambda*D)
mu <- rep(global.mu, times=side*side)
# sampling from the multivariate normal distribution
M <- matrix(nrow=side, ncol=side)
M[] <- rmvnorm(1, mu, SIGMA)
return(M)
}
# parameters (the truth) that I will want to recover by JAGS
side = 10
global.mu = 0
lambda = 0.3 # let's try something new
# simulating the main raster that I will analyze as data
M <- cor.surface(side = side, lambda = lambda, global.mu = global.mu)
image(M)
mean(M)
# simulating the inherent uncertainty of the mean of M:
#test = replicate(1000, mean(cor.surface(side = side, lambda = lambda, global.mu = global.mu)))
#hist(test, breaks = 40)
# preparing the data for JAGS
y <- as.vector(as.matrix(M))
my.data <- list(N = side * side, D = dist.matrix(side), y = y)
cat("
model
{
# priors
lambda ~ dgamma(1, 0.1)
global.mu ~ dnorm(0, 0.01)
for(i in 1:N)
{
# vector of mvnorm means mu
mu[i] <- global.mu
}
# derived quantities
for(i in 1:N)
{
for(j in 1:N)
{
# turning the distance matrix to covariance matrix
D.covar[i,j] <- exp(-lambda*D[i,j])
}
}
# turning covariances into precisions (that's how I understand it)
D.tau[1:N,1:N] <- inverse(D.covar[1:N,1:N])
# likelihood
y[1:N] ~ dmnorm(mu[], D.tau[,])
}
", file="mvnormal.txt")
fit <- jags(data=my.data,
parameters.to.save=c("lambda", "global.mu"),
model.file="mvnormal.txt",
n.iter=10000,
n.chains=3,
n.burnin=5000,
n.thin=5,
DIC=FALSE)
plot(as.mcmc(fit))
pairs(as.matrix(as.mcmc(fit)))
@MushagalusaCz
Copy link

Dear Florian Greetings!
You've created a function that simulates an autocorrelated 2D array with a specific side and a lambda-based exponential decay.
I'm curious if lambda and the Moran autocorrelation index have any relationship.
What is the nature of this connection, if it exists?
What can I do to change the function such that lambda remains constant across the 2D array?
Kind regards

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment