Error in generating a geometric random variable

General
1

I am trying to generate a Geometric random variable. My code is:

generated quantities{

real theta[I];
int lifetime[I];

for (i in 1:I) {
theta[i] = beta_rng(a, b);
lifetime[i] = neg_binomial_rng(1 ,theta[i]/(1-theta[i]) ) ;
}

}

But I get the following error:

Exception: neg_binomial_rng: Random number that came from gamma distribution is 5.58759e+09, but must be less than 1.07374e+09

I believe Stan is trying to generate a really large number. Is there a way to preempt this error by telling Stan to give some number (that I am ok with) when this issue happens? Thanks.

EDIT: I extracted theta and I see that some values are very small. Is there a way to circumvent those small values by placing a lower limit (say, 1e-7)?

2

You can constrain variables using for instance real<lower=0, upper=1> theta[L] (unless I’m mistaken that includes 0 and 1, so you may use something small/close to 1 instead), since the variables you are having trouble with seem to be only in the generated quantities block it won’t even affect inference and should solve your problem.

3

Hi, I tried your advice by using:

real<lower=1e-7, upper=1> theta[I];

But I still got this error:

[1] "Error in sampler$call_sampler(args_list[[i]]) : "
[2] " Exception: modelbdb84c3dc631_sBG_namespace::write_array: theta[k0__] is 3.1914e-08, but must be greater than or equal to 1e-07 "

Instead, I use a brute-force method after I sample theta:

for (i in 1:I) {
theta[i] = beta_rng(a, b);
if (theta[i] < 1e-7) theta[i]=1e-7;
lifetime[i] = neg_binomial_rng(1 ,theta[i]/(1-theta[i]) ) ;
}

I got no more error. I believe that’s the way to go. But thank you for chiming in and giving me that suggestion.

However, could someone tell me whether my code could be vectorized?

4

That doesn’t make sense to me, since you clearly restricted the variable to values above 1e-7. I don’t know if there’s something else going on there, maybe a developer can comment.

Stan is basically C++, so vectorizing usually means simply writing a function with the loop you need. There is probably no difference for most functions.

5

If the parameter constraint is placed in the generated quantities block then it really only acts as a means of error-checking, rather than constraining the number generation itself.

For example, given the original block:

generated quantities{
  real theta[I];
  int lifetime[I];

  for (i in 1:I) {
    theta[i] = beta_rng(a, b);
    lifetime[i] = neg_binomial_rng(1 ,theta[i]/(1-theta[i]) ) ;
  }
}

Putting a <lower=xx> constraint on theta won’t affect the generated number, because that is solely determined by beta_rng(a, b). All it does is throw an error if the generated number does not meet the constraint (which is what happened above).

1 Like
7

Are you sure it happens when theta is small? You can only get huge numbers of \alpha = \theta / (1 - \theta) if \theta is near 1.

Just multiply the beta_rng by 1 - \epsilon to bound from above

...
theta[i] =  (1 - 1e-7) * beta_rng(a, b);
8

Thank you for everyone’s response.

I solved the problem by writing my own geometric rng:

real geometric_rng(real theta) {
real u = uniform_rng(0, 1);
real x = floor( log(1-u)/log(1-theta)) ;
return x;
}

This works without having to place any lower bound.