Reparametrization for parameter restriction


#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])

#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!!