Hi,
I’m worked through the Rethinking book and I’m a beginner with STAN and also with Gaussian processes. I’ve generated simulated data and now want to estimate the parameters used for the simulation. With the help of this forum, the simulation is working. Now I struggling with getting good fit. I wondering if someone can point me to the mistake I did. Here is my Code:
library(rethinking)
library(MASS)
# Generate distance matrix
set.seed(5)
N_PARTICLE=10
particle_positions <- matrix(runif(2*N_PARTICLE, 0, 1),nrow=N_PARTICLE, ncol=2)
particle_distance_matrix <- as.matrix(dist(particle_positions))
# These are the parameters I want to get back:
eta_sq = 3 # Maximum covariance between any two particles
rho_sq = 10 # rate of decline ofcorrelation with distance
delta = 0.01 # don't now what it is, but this is used in the book
# Generate the covariance matrix
covariance_matrix <- matrix(, nrow = N_PARTICLE, ncol = N_PARTICLE)
# Using gaussian kernel
for (i in 1:(N_PARTICLE-1)) {
covariance_matrix[i, i] <- eta_sq + delta;
for (j in (i+1):N_PARTICLE) {
covariance_matrix[i,j] = eta_sq * exp(-rho_sq*(particle_distance_matrix[i,j])^2)
covariance_matrix[j,i] = covariance_matrix[i,j]
}
}
covariance_matrix[N_PARTICLE, N_PARTICLE] = eta_sq + delta;
# Generate samples
N_SAMPLES <- 50
shifts <- mvrnorm(n=N_SAMPLES,mu=rep(0, N_PARTICLE), Sigma=covariance_matrix)
shifts <- as.vector(t(shifts)) # Flatten the matrix
particle_id <- rep(1:N_PARTICLE,N_SAMPLES)
# Fit the gaussian process
dat_list <- list(
S = shifts,
N = N_PARTICLE,
particle = particle_id,
Dmat = particle_distance_matrix )
m14.8 <- ulam(
alist(
S <- k[particle],
vector[N]:k ~ multi_normal( 0 , SIGMA ),
matrix[N,N]:SIGMA <- cov_GPL2( Dmat , etasq , rhosq , 0.01 ),
etasq ~ dexp( 2 ),
rhosq ~ dexp( 0.5 )
), data=dat_list , chains=4 , cores=4 , iter=2000 )
precis(m14.8)
The estimated parameters are:
Which is obviously different from what I’ve used during simulation. What did I wrong?