Fitting data as a series of univariate Normal distributions

Modeling
1

Dear Stan users,

I’m trying to fit a model that describes my data as a series of univariate Normal distributions with varying mean and SD depending on a parameter SNR.

I hypothesize that for each SNR level (ranging from 2 to 5), the corresponding data would be distributed according to a Normal distribution with its own mean and SD, and that these means and SDs would vary monotonously with SNR within bounds.

This is the model I wrote.

There are strong convergence issues which have me guessing I wrote something poorly, or plain wrong.

I added the output of my model diagnostics below the Stan code here (default fitting options besides adapt_delta 0.99).

Could you point me in the right direction? I’m not sure what to do next.


functions{
}

data {
  int Nobs; // number of data obervations
  real mu_low; //prior knowledge from external data
  real sigma_low; //prior knowledge from external data
  vector[Nobs] dat; // data, one per row
  int<lower = 2, upper = 6> SNR[Nobs]; // experimental factor
}


parameters {
  vector[5] b;// unconstrained parameters
  real<lower = 0> pos_b; // postive parameters
}

transformed parameters{
real mu;
real sigma;

for(i in 1:Nobs){
     mu = b[1] + (b[2]-b[1])/(1+exp(-b[3]*(SNR[i] - b[4])));
     sigma = pos_b/(1+exp(-b[5]*(SNR[i] - b[4])));
}

}
model {
  
  

  b[1] ~ normal(mu_low,30); // min
  b[2] ~ normal(mu_low,30); // max
  b[3] ~ normal(0,10); // slope
  b[4] ~ normal(4,10); // threshold
  b[5] ~ normal(0,10); // slope sigma 
  pos_b ~ normal(sigma_low,30); // max sigma
 
  
  //LLH
  for(i in 1:Nobs){
    target += normal_lpdf(dat[i] | mu,sigma);
    }
}

generated quantities{
  vector[Nobs] log_lik;

    for (i in 1:Nobs){
    log_lik[i] = normal_lpdf(dat[i]|mu, sigma);
    }
}

Warning messages:
1: There were 34 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup 
2: There were 492 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded 
3: Examine the pairs() plot to diagnose sampling problems

>print(uni_stan, pars = c("b"), probs = c(0.025, 0.975),digits=3)

Inference for Stan model: unimodal.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

       mean se_mean     sd    2.5%  97.5% n_eff  Rhat
b[1] 23.534   4.957 37.228 -57.559 57.219    56 1.050
b[2] 21.662   5.535 35.624 -56.190 56.867    41 1.109
b[3]  1.692   1.235  9.941 -18.292 21.061    65 1.041
b[4]  1.462   4.955  9.834 -16.486 21.333     4 1.605
b[5]  1.200   4.911  9.773 -18.194 19.233     4 1.584

Samples were drawn using NUTS(diag_e) at Wed Mar 25 16:05:13 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
2

I believe you want mu and sigma to be vectors, not a single number each, so you would do:

transformed parameters{
    vector[Nobs] mu;
    vector[Nobs] sigma;
    for(i in 1:Nobs){
        mu[i] = b[1] + (b[2]-b[1])/(1+exp(-b[3]*(SNR[i] - b[4])));
        sigma[i] = pos_b/(1+exp(-b[5]*(SNR[i] - b[4])));
    }
}

And in the model block you can do just:

    target += normal_lpdf(dat | mu,sigma);

(i.e. no for loop in the model block)

3

That sounds fair, thank you. I made the modification.

Unfortunately this doesn’t change much…

Here are the trace and pairs plots for the b parameters.

Rplot03
Rplot03813×796 243 KB

Rplot04
Rplot04813×796 205 KB

I’m a bit lost with model diagnostics beyond figuring out that thing are going wrong. But in this case it seems that the upper and lower bounds for mu are interchangeable, and that the slope switches sign depending on which value is attributed to which parameter (either around 80 or 0).

I assume this issue comes from my definition of the generalized logistic function for mu?
I switched

mu[i] = b[1] + (b[2]-b[1])/(1+exp(-b[3]*(SNR[i] - b[4])));

for

mu[i] = b[1] + (b[2]/(1+exp(-b[3]*(SNR[i] - b[4]))));

But still no luck.

4

Looks like an identifiability issue with parameters b[1], b[2], and b[3]. Is the math for the underlying model something you developed yourself, or is there any prior work we could look at to get a better sense of what you’re trying to achieve?

1 Like
5

The math is an expression of a sigmoid function with a varying slope b[3], x50 b[4] and bounds b[1] and b[2].
One reference where it is used in a similar manner to mine is this one for example (search for “nonlinear regression” in the text and you’ll find it more easily) but it’s also available on Wikipedia in its full form.

Maybe this isn’t the correct way to write this in Stan?

6

Alright, so:

Since it appears that b[1] >b[2] with a positive slope is the same as having b[2] > b[1] with a negative slope in terms of model fit, I simply got things worked out by constraining my slope to be positive. Now my b parameters are no longer interchangeable, and can only take the value of one of the two bounds, depending on whether mu increases or decreases in the data.

There are no identifiability issues left, and this doesn’t impair my ability to interpret the parameters of the model.

Let me know if that sounds fishy to you, otherwise everything is good on my end!
Thank you for the help!

1 Like
7

Sounds good to me. Happy you worked it out!