Recovering parameters of fake model

Modeling
1

Hello Stan users,

I have created 640 observations (data is attached in CSV format) in the following form:

Y = alpha - u + error.

where alpha is intercept with value of 0.10
u comes from a truncated normal distribution [0, ] with mean (mu_u) of 0.20 and standard deviation (sigma_u) of 0.15
and error comes from a normal distribution of mean zero and standard deviation (sigma) of 0.13.

I am trying to recover the values of alpha, mu_u, sigma_u and sigma.

However I am having issues in recovering the model parameters with divergences and recovered parameters way off from what I created.

Here is the output of my recovery model.

image
image1233×385 18 KB

Here is the Stan code I am using


data {
int<lower=0> N; // Number of obs
vector[N] Y;
}

parameters {
real mu_u;
real alpha;
real<lower=0> sigma;
real<lower=0> sigma_u;
vector<lower=0>[N] u; 
}


transformed parameters{
vector[N] yhat;
yhat = alpha - u;
}

model {

alpha ~ normal(0,1);

sigma ~ cauchy(0, 1);


target += normal_lpdf(u | mu_u, sigma_u);

for (n in 1:N) 
if (u[n] < 0)
 target += negative_infinity();
else
target += normal_lpdf(mu_u |0,1);
target += cauchy_lpdf(sigma_u | 0,1);
target += -normal_lccdf(0 | mu_u, sigma_u);





Y ~ normal(yhat, sigma);

}


generated quantities {

}

Can someone please help with this recovery process. I have spend couple of days working on this with no success.

Thanks
simulated_data_set.csv (10.2 KB)

2

To double-check things I recreated the synthetic data according to the specifications in your post:

u \sim N(\mu_u,\sigma_u) \text{ T[0,]} \\ \mu_u = 0.20, \sigma_u = 0.15 \\ \varepsilon \sim HalfNormal(\sigma) \\ Y = \alpha -u + \varepsilon

Which in R is:

N <- 640 
alpha <- 0.1
mu_u <- 0.20
sigma_u <- 0.15
lower_bound <- 0
sigma <- 0.13

u <- truncnorm::rtruncnorm(n = N, a = lower_bound, mean = mu_u, sd = sigma_u)
error <- extraDistr::rhnorm(n = N, sigma = sigma)

Y = alpha - u + error

Because your generating parameters are so close to 0, but are also lower-bounded at 0, the sampler is going to have a difficult time with the posterior.

To make things a bit easier, we can use the non-centered parameterisation for u (an explanation of how to specify the bounds is in this forum post), and use tighter priors for the parameters:

data {
  int<lower=0> N;
  vector[N] Y;
}

parameters {
  real<lower=0> mu_u;
  real alpha;
  real<lower=0> sigma;
  real<lower=0> sigma_u;
  vector<lower=-mu_u/sigma_u>[N] u_raw;
}

transformed parameters {
  // Implies u ~ normal(mu_u, sigma_u) with lower bound of 0
  vector<lower=0>[N] u = mu_u + u_raw * sigma_u;
  // Calculate density adjustment for lower truncation at 0 for
  // every observation
  real trunc_adjustment = -N * normal_lccdf(0 | mu_u, sigma_u);
}

model {
  alpha ~ normal(0, 0.5);
  mu_u ~ normal(0, 0.5);
  sigma_u ~ normal(0, 0.5);
  sigma ~ normal(0, 0.5);

  u_raw ~ std_normal();
  target += trunc_adjustment;

  Y ~ normal(alpha - u, sigma);
}

When running the model, you’ll also need to increase the adapt_delta and max_treedepth settings, so that the sampler has to operate with more precision when traversing the posterior:

mod <- cmdstan_model("trunc.stan")
samp <- mod$sample(data = list(N = N, Y = Y),
                   adapt_delta = 0.999,
                   max_treedepth = 12)

This appears to recover the parameters reasonably well:

> samp$summary(c("mu_u", "alpha", "sigma", "sigma_u"), c("mean", "rhat", "ess_bulk"))
# A tibble: 4 x 4
  variable  mean  rhat ess_bulk
  <chr>    <dbl> <dbl>    <dbl>
1 mu_u     0.123  1.74     6.12
2 alpha    0.135  1.61     6.69
3 sigma    0.113  1.43     8.37
4 sigma_u  0.114  1.04   127.

The chains haven’t completely converged, so you may need to increase the adapt_delta even further

2 Likes