Hi Ben,
I tried going down in cores and using a single core (which crashed the R session completely). I also tried increasing memory and monitoring it’s usage, it picks up when the program starts running but comes nowhere near the limit which makes me think there must be something wrong with the model itself.
I’m trying to run build a hierarchy into a random walk where I have confirmed accurate measurements at the start/end of the period and then noisy, irregularly-spaced measurements in between. I also have confirmed measurements at a lower district level at the start and end of the period, but no measurements in between. I’m trying to model a random walk for each district between the start and end measurements where a weighted average of the district measurements should equal the region estimate for each day and the region random walk should move from the starting measurement to the finishing measurement subject to the noisy measurements. Also there are two distinct time periods where the whole process repeats itself (new start/finish, etc).
Is there anything particularly egregious or bad practice that I’m doing here? Thanks for any input you can give me!
data{
//read in dimensions
int N; //number of measurments
int p; //number of random walks to be estimated - measuring p distinct quantities for now uncorrelated
int N_k; //number of distinct time periods - in this case there are 2
int station_k_id; //number of measurement stations over N_k periods (if same station is used twice it is counted twice)
int s[station_k_id]; //position of measurement in time, used to match measurements to mu
int total_days; //number of days in all the time periods
int n_days[N_k]; //number of days in each time period
int day_index[N]; //index for matching measurements to correct place in mu
int N_station; //number of unique stations
int station_id[station_k_id]; //unique station id
int last_day_index[N_k]; //index in mu for end of distinct time period
matrix[N_k,p] mu_start; //starting values
matrix[N_k,p] mu_finish; //finishing values
matrix[N,p] y; //measurements
matrix[N,p] y_se; //se of measurements
int R; //number of districts
row_vector[R] weights; //weights of districts
int district_day_index[R,total_days]; //index for district measurements in mu_district - row for each district
matrix[N_k*R,p] mu_start_district; //confirmed starting values district
matrix [N_k*R,p] mu_finish_district; //confirmed finishing values district
int mu_finish_district_index[R*N_k];
}
parameters{
matrix[N_station,p] systematic_error_std;
matrix[total_days,p] innovation;
matrix[total_days*R,p] innovation2;
real<lower=0> sigma;
real<lower=0> sigma2;
real<lower=0> sd_inflator[N_station];
real<lower=0> fin_sd_std;
real<lower=0> fin_sd_std2;
}
transformed parameters{
matrix[total_days,p] mu; //regional estimate
matrix[total_days*R,p] mu_district; //district estimate - only have starting and finishing values
matrix[N_station,p] systematic_error = systematic_error_std/20; //systematic error in measurement stations
real<lower=0> fin_sd = fin_sd_std/100;
real<lower=0> fin_sd2 = fin_sd_std2/100;
{
int pos = 1;
for (r in 0:(R-1)){
int st = 1;
for (k in 1:N_k) {
mu_district[pos,] = mu_start_district[st+r,]; // fixing starting period of estimate to be confirmed starting estimate
for (t in (pos+1):(pos+n_days[k]-1)){
for (i in 1:p){
mu_district[t,i] = mu_district[t-1,i] + innovation2[t,i]*sigma2; //district random walk
}
}
pos = pos + n_days[k];
st = st + R;
}
}
}
{
int pos = 1;
for (k in 1:N_k) {
mu[pos,] = mu_start[k,]; // fixing starting period of estimate to be confirmed starting estimate
for (t in (pos+1):(pos+n_days[k]-1)){
for (i in 1:p){
mu[t,i]=dot_product(weights,mu_district[district_day_index[,t-1],i])+innovation[t,i]*sigma; //regional estimate is random walk of weighted average of district estimates
}
}
pos = pos + n_days[k];
}
}
}
model{
to_vector(innovation) ~ student_t(4, 0, 1);
to_vector(innovation2) ~ student_t(4, 0, 1);
fin_sd_std~normal(0,1);
fin_sd_std2~normal(0,1);
to_vector(systematic_error_std) ~ normal(0,1);
sigma ~ normal(0,1);
sigma2 ~ normal(0,1);
sd_inflator ~ normal(1,5);
// fixing final period of estimate to be very close to actual confirmed final estimate
for (i in 1:p){
mu_finish[,i] ~ normal(mu[last_day_index,i], fin_sd);
mu_finish_district[,i] ~ normal(mu_district[mu_finish_district_index,i], fin_sd2);
}
// Measurements at the regional level
{
int pos2= 1;
for (k in 1:station_k_id) {
for (i in 1:p){
segment(y[,i], pos2, s[k]) ~ normal(mu[segment(day_index, pos2, s[k]),i]+systematic_error[station_id[k],i], sd_inflator[station_id[k]]*segment(y_se[,i], pos2, s[k]));
}
pos2 = pos2 + s[k];
}
}
}