Reparametrization for parameter restriction

Modeling
1

(corrected )
I am using a model very similar to Noisy Categorical Measurement Model in the Stan reference 14.4.
I have a parameter vector with variant range (theta), so I transformed it to alpha to restrict the range with mu > 0. Mu is the poisson parameter. But when sampling, I got the error:

Error in sampler$call_sampler(c(args, dotlist)) :
Exception: poisson_lpmf: Rate parameter is -1.05972, but must be >= 0! (in ‘modele8162f3d0378_categorical’ at line 20)

How can this happen??
Here is my code:

data {
  int N;
  int K;
  int Y[N];
}

parameters {
  vector[K] theta;
  simplex[K] pi;
  vector<lower=0>[K] alpha;
}

transformed parameters {
  vector<lower=0>[K] mu;
  vector[K] log_q_z[N];
  for(n in 1:N){
    log_q_z[n] = log(pi);
    for(k in 1:K){
      mu[k] = 2^k-theta[k];
      log_q_z[n,k] = log_q_z[n,k]  + poisson_lpmf(Y[n] | mu[k]);
    }
  }
}

model {
  pi ~ dirichlet(alpha);
  theta ~ cauchy(0,3);
    
  for(n in 1:N)
    target += log_sum_exp(log_q_z[n]);
}

Then R code:

model_categorical <- stan_model("~/Project/montagna_mar_collab/Shuonan/categorical.stan")
Y <- data_9_Y;
data_temp <- list(N = length(Y), K = 4, Y = Y)
fit_9_Y <- sampling(model_categorical, data = data_temp,
                   iter = 1000, chains = 4, cores = 4, control = list(adapt_delta = 0.9))

I read the reparameterization chapter in the stan reference (20.4) but I am not familiar with Jacobian matrix so would really appreciate if theres a alternative approach :)
Thank you !

2

When k is 1 and theta[k] is greater than 2, then mu[k] = 2^k-theta[k] is negative, hence the error.

3

Sorry I had to correct my code. It should not have that restriction on theta. But I have restriction on mu >0 (not alpha, sorry for confusion!) . I had some typo issue with the original description.

4

Specifying

does not make it so; it just means the inequality restriction on each element of mu is checked at the end of the transformed parameters block. You still need a construction of mu[k] = 2^k-theta[k] that makes the restriction hold, such as vector<upper=1>[K] theta;

5

I am not sure if I get it… so mu is checked for after the construction of mu[k] = 2^k-theta[k] , and if this is <0, then it is rejected. So how I can still get mu <0 in the sampling?

6

No, you don’t want to ever be in the position of doing rejections. That will cause Stan to reduce to a random walk and slow to a crawl and potentially even fail to fit in finite time.

What you want to do is construct variables that satisfy the constraints. If you require mu = 2^k - theta[k] > 0 then you require 2^k > theta[k]. There’s a section of the manual that explains how to deal with arrays of parameters with varying constraints, but in this case, it’s simple:

vector<upper = 0>[K] theta_raw;

theta[k] = theta_raw + 2^k;

Then you know theta[k] < 2^k and you’re good to go. You don’t need a Jacobian as you’re only adding a constant.

7

Although in this case you’ll probably just want to exponentiate your current calculation of the rate,

poisson_lpmf(Y[n] | exp(mu[k]))

which can more conveniently be written as

poisson_log_lpmf(Y[n] | mu[k])
1 Like
8

Thank you! But I was trying to use mu to do the same thing I think??

mu[k]= 2^k-theta[k] > 0
and
theta_raw[k] = theta[k]-2^k < 0

9

But actually I didnt get error on the parameter range using your method. What is the differences on these two??

10

I’m not sure what your approach was—what you wrote down in the last reply isn’t Stan code.

11

Yes,
I used mu = 2^k - theta[k] and restricted mu to be > 0.
But I think you used theta_raw = theta[k] - 2^k and restricted to <0.
When sampling, the first approach gives error that poisson parameter is < 0 (even though I set it to be > 0)

12

You want to make sure that every value of parameters that satisfy the declared constraints leads to a finite log density. You do not want to try to guess and reject.

13

Got it. Very clear.
Thank you so much!!!

14

Thank you this is brilliant.
Thank you so much for all your help!!