 # Error in generating a geometric random variable

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

generated quantities{

real theta[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)?

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.

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

But I still got this error:

 "Error in sampler\$call_sampler(args_list[[i]]) : "
 " 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?

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.

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];

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

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);


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.