issue with neg_binomial_rng

[Linas Mockus]

Hi,

It looks like that neg_binomial_rng is still creating problems.  I have compared with equivalent R code and numbers don't match - neg_binomial_rng gives lots of zerors while R gives reasonable numbers.  The workaround poisson_rng (gamma_rng (shape, inv_scale)) provided expected results.

Linas

[Daniel Lee]

Hi Linas,

We'll look into it when we have a chance. I've put up an issue here:
https://github.com/stan-dev/stan/issues/591

Do you have a small example that shows the behavior? I guess I'm interested in values of alpha and beta that are causing issues. That'll help us pinpoint faster.

Also, what version of Stan are you using? And is it in RStan, CmdStan, PyStan, MStan?

--
You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Linas Mockus]

Daniel,

I am using Rstan < 2.2.0 (something from December last year).

The code is:
    data {
        int<lower=0> N;
        int<lower=0> NL;
        int<lower=0> lot [N];
        int<lower=0> counts [N];
    }

    transformed data {
    }

    parameters {
        real alpha0;
        real<lower=0> sigma2_alpha;
        real<lower=0> sigma2_gamma;
        real alpha [NL];
        real<lower=0> shape;
        real gamma [N];
    }

    transformed parameters {
        real<lower=0> sigma_alpha;
        real<lower=0> sigma_gamma;
        vector[N] log_lik;
        real mu [N];

        sigma_alpha <- sqrt (sigma2_alpha);
        sigma_gamma <- sqrt (sigma2_gamma);
        for (i in 1 : N) {
            mu [i] <- exp (alpha0 + alpha [lot [i]] + gamma [i]);
            log_lik [i] <- neg_binomial_log (counts [i], shape,
                shape / mu [i]);
        }
    }

    model {
        alpha0 ~ normal (0, 100);
        sigma2_alpha ~ cauchy (0, 5);
        alpha ~ normal (0, sigma_alpha);
        gamma ~ normal (0, sigma_gamma);
        shape ~ cauchy (0, 5);
        increment_log_prob (sum (log_lik));
    }
    generated quantities {
        real alpha_rep;
        real gamma_rep [5];
        real y_rep [5];
        real y_mean;
        real mu_rep [5];

        alpha_rep <- normal_rng (0, sigma_alpha);
        for (i in 1 : 5) {
            gamma_rep [i] <- normal_rng (0, sigma_gamma);
            mu_rep [i] <- exp (alpha0 + alpha_rep + gamma_rep [i]);
            y_rep [i] <- poisson_rng (gamma_rng (shape,
                shape / mu_rep [i]));
        }
        y_mean <- mean (y_rep);
    }

Linas

[Daniel Lee]

Thanks.

Can you let me know what some of the values are for:

- shape

- shape / mu_rep[i]

I really don't need the rest of the model to hunt it down.

[Linas Mockus]

shape
[1] 80.02484 139.35417 374.26470 69.21295 328.44291 59.54400
[7] 77.25161 131.09119 162.37569 253.16368 196.94005 150.54016
[13] 173.35594 50.64902 95.52137 192.47620 200.44689 34.56073
[19] 108.97727 54.90816 50.29049 174.58498 99.90834 211.67859
[25] 119.34693 98.97055 116.58883 45.58940 82.03280 47.81400
[31] 44.86904 33.33081 33.59745 120.16425 38.88848 100.14836
[37] 77.96474 130.17035 47.83731 59.12910 102.48124 230.63175
[43] 53.47698 46.76760 256.13068 56.42462 69.89229 413.19436
[49] 119.45114 43.67303 149.06586 143.65760 72.00611 89.23581
[55] 55.81904 50.33290 61.07071 70.66281 49.39308 70.80331
[61] 1032.27123 109.99165 62.45677 44.99347 520.74048 46.20570
[67] 86.76055 75.16100 56.82449 55.54931 143.70130 35.96519
[73] 172.98580 177.64758 40.12357 46.03435 204.98821 58.46726
[79] 42.56416 73.43106 413.07517 59.54960 69.93423 93.64422
mu[,1]
[1] 50880.36 57035.17 50782.69 55098.73 51243.42 57492.07 57932.20
58003.16
[9] 47981.74 53153.36 50550.82 46628.83 42436.63 54503.56 50632.47
46575.60
[17] 49846.08 57097.84 46527.47 58215.96 63401.36 56723.12 52193.80
49048.57
[25] 48042.59 50427.96 54367.64 59929.64 51685.83 73599.43 55302.54
58229.82
[33] 55296.21 57613.69 57786.18 48594.53 51415.92 51531.23 53722.19
58203.74
[41] 52971.13 52122.97 54019.32 56929.22 53744.41 57938.09 61334.14
49663.89
[49] 48713.33 56967.25 52545.77 52810.82 58053.29 55292.09 63808.33
53059.92
[57] 45798.07 54879.23 70894.16 52032.56 48187.63 52962.99 58913.42
60566.13
[65] 47048.43 51972.81 56347.79 47325.26 60547.63 54822.56 49920.05
73323.24
[73] 53249.80 54065.75 56536.99 60094.83 48781.26 52094.53 54998.14
47422.96
[81] 48050.97 54309.81 45126.12 65160.73 54749.13 44417.60 48162.06
59595.60
[89] 50300.42 47064.74 54504.89 60039.92 50821.99 55994.00 51647.36
48089.80

[Linas Mockus]

mu_rep[,1]
[1] 52815.52 46725.17 61551.10 55025.72 64081.78 61435.34 57657.41
[8] 58321.70 65835.66 70798.10 50653.66 61981.91 43356.51 64246.17
[15] 53283.04 58579.65 68806.42 61680.05 55046.92 60659.20 54415.93
[22] 69949.37 60422.99 54622.61 57185.91 73470.99 61242.68 58583.11
[29] 62977.18 58251.86 59542.92 62168.12 63858.80 57786.28 55716.93
[36] 51219.75 55783.80 110076.04 58059.43 59909.71 67627.66 66719.20
[43] 84076.06 76679.97 69603.53 56866.39 62308.10 49782.20 62342.84
[50] 59126.71 79301.26 63274.85 54982.22 63694.55 43583.13 67814.96
[57] 61423.22 65836.70 81815.56 59272.33 68133.80 57695.95 62230.10
[64] 63242.33 86610.44 57283.06 60493.90 57017.10 59115.34 54262.59
[71] 60091.80 56619.94 69551.97 48791.12 61940.03 59442.99 62303.81
[78] 52749.81 62528.53 54123.53 46170.60 52273.64 68233.80 56483.48
[85] 67623.09 70118.37 64010.85 63597.00 51093.40 47937.24 53049.47
[92] 124562.62 58900.35 51291.35 52365.57 62565.32 106879.13 73072.87
[99] 61827.42 55060.30

[Daniel Lee]

Mind dumping that to an R file? So it looks something like "shape <- c(...)"

You can either use stan_rdump() or dump(). (It's hard to "read" that format straight in.)

Daniel

--
You received this message because you are subscribed to the Google Groups "stan users mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

[Linas Mockus]

[Avraham Adler]

I'm not sure this has anything to do with the problem, but I think the definition of the NB as listed in the manual is incorrect, unless a different parameterization is being used.

I believe the combinatorial prefix should be:

Choose(n + a - 1, a)

instead of:

Choose(n + a - 1, a - 1)

Thus, the logs should be:

lgamma(a + n) - lgmma(a + 1) - lgamma(n)

instead of:

lgamma(a + n) - lgmma(n + 1) - lgamma(a)

and the errors propagate to the gradient definitions.

If this is what is being programmed, that may be the source of some issues.

Avi

On Friday, March 7, 2014 1:12:52 PM UTC-5, Linas Mockus wrote:

[Avraham Adler]

On further review, it may well be that this could be due to R and Stan using different parameterizations. If Stan is actually following Gelman et al (2014 p. 44), then it is self consistent, but has defined successes and failures opposite that of R (check ?rbinom). Stan seems to be using alpha is the number of set failures for which we are finding the probability of y successes. R, on the other hand, has x as the random successes before a fixed number of n failures. In other words, "p" and (1-p) are reversed so that the exponent to the 1 / 1+ b and b / 1 + b have to be switched between the two.

Compare Wikipedia with Wolfram for an example where the exponents are "swapped" because of the definitions.

Linus, using your Stan parametrization, did you generate variates in R using size and mu?

[Linas Mockus]

Yes, there are parametrization challenges.  I have used size and mu parametrization in R.

Thank you,
Linas
765-494-0344 

On 3/7/2014 2:07 PM, Avraham Adler wrote:

On further review, it may well be that this could be due to R and Stan using different parameterizations. If Stan is actually following Gelman et al (2014 p. 44), then it is self consistent, but has defined successes and failures opposite that of R (check ?rbinom). Stan is using the same formulation as found in Wikipedia, in that alpha is the number of set failures for which we are finding the probability of y successes. R, on the other hand, has x as the random successes before a fixed number of n failures. In other words, "p" and (1-p) are reversed so that the exponent to the 1 / 1+ b and b / 1 + b have to be switched between the two.

Linus, using your Stan parametrization, did you generate variates in R using size and mu?


On Friday, March 7, 2014 1:33:17 PM UTC-5, Avraham Adler wrote:

--
You received this message because you are subscribed to a topic in the Google Groups "stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/8ILE6JLcENU/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Avraham Adler]

In R, size, if I recall correctly, is beta, = (1 - p) p for the R/Wolfram/Klugman Panjer Wilmot/Degroot etc. parametrization, which is also the beta from the Gamma prior. Were you using Stan's Gamma's beta for size? You may need its reciprocal if my mental algebra is correct.

[Linas Mockus]

shape=alpha=r
mu=(1-p)/p=beta
thus in stan neg_binomial(shape,shape/mu) while in R rnbinom(size=shape,mu=mu)

Thank you,
Linas
765-494-0344 

[Bob Carpenter]

Stan is indeed using the parameterization of Gelman et al.'s Bayesian
Data Analysis. If that's not what it is, there's a bug somewhere.

- Bob

On Mar 7, 2014, at 8:07 PM, Avraham Adler <avraha...@gmail.com> wrote:

> On further review, it may well be that this could be due to R and Stan using different parameterizations. If Stan is actually following Gelman et al (2014 p. 44), then it is self consistent, but has defined successes and failures opposite that of R (check ?rbinom). Stan is using the same formulation as found in Wikipedia, in that alpha is the number of set failures for which we are finding the probability of y successes. R, on the other hand, has x as the random successes before a fixed number of n failures. In other words, "p" and (1-p) are reversed so that the exponent to the 1 / 1+ b and b / 1 + b have to be switched between the two.
>
> Linus, using your Stan parametrization, did you generate variates in R using size and mu?
>
>
> On Friday, March 7, 2014 1:33:17 PM UTC-5, Avraham Adler wrote:

> --
> You received this message because you are subscribed to the Google Groups "stan users mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Daniel Lee]

I think there's a bit of a problem in the random number generator.

I'll try to verify and fix it soon.


Daniel

[Marco Inacio]

Hi, the parametrization would be:

alpha=r=shape (indeed)

and

mu = alpha/beta -> beta = alpha/mu



In your code (workaround) you use:

y_rep [i] <- poisson_rng (gamma_rng (shape,
shape / mu_rep [i]));

which means
y_rep [i] <- poisson_rng (gamma_rng (alpha,
beta));



If your code is correct, than stan's isn't:

neg_binomial_rng(const double alpha,
const double beta,
RNG& rng) {
(...)

return stan::prob::poisson_rng(stan::prob::gamma_rng(alpha, 1.0 /
beta,
rng),rng);

>>> <https://groups.google.com/d/topic/stan-users/8ILE6JLcENU/unsubscribe>.

>>> To unsubscribe from this group and all its topics, send an email

>>> to stan-users+...@googlegroups.com <javascript:>.

>>> For more options, visit https://groups.google.com/d/optout

>>> <https://groups.google.com/d/optout>.

>>
>> --
>> You received this message because you are subscribed to a topic in
>> the Google Groups "stan users mailing list" group.
>> To unsubscribe from this topic, visit
>> https://groups.google.com/d/topic/stan-users/8ILE6JLcENU/unsubscribe.
>> To unsubscribe from this group and all its topics, send an email to
>> stan-users+...@googlegroups.com

>> <mailto:stan-users+...@googlegroups.com>.

[Marco Inacio]

Here's a minimal reproducible example:

require(rstan)
fit <- stan(model_code=
'data {
real shape;
real mu;
}
parameters {
real any; //just a dummy to make this compile
}
model {
any ~ normal(0, 1);
}
generated quantities {
real y_correct[1000];
real y_incorrect[1000];

for (i in 1:1000) {
y_correct [i] <- poisson_rng (gamma_rng (shape, shape/mu));
}

for (i in 1:1000) {
y_incorrect [i] <- neg_binomial_rng(shape, shape/mu);
}
}', data=list(shape=4, mu=2))

as.matrix(fit) -> mfit
correct_sample <- mfit[, grep("y_correct", colnames(mfit))]
incorrect_sample <- mfit[, grep("y_incorrect", colnames(mfit))]

r_test <- rnbinom(1000, mu=2, size=4)

mean(correct_sample)
mean(incorrect_sample)
mean(r_test)
sd(correct_sample)
sd(incorrect_sample)
sd(r_test)

[Marco Inacio]

hope you don't mind I have pushed a branch feature/fix-issue591-nb_rng
that might fix the problem.

[Bob Carpenter]

That's great. We love other people fixing the problems.

Our goal is to have both the RNG and the unmarked distribution
("neg_binomial") be the one from Gelman et al.'s Bayesian Data
Analysis. Is that what the fix entails?

- Bob

[Marco Inacio]

On 14-03-08 04:24 PM, Bob Carpenter wrote:
> That's great. We love other people fixing the problems.
>
> Our goal is to have both the RNG and the unmarked distribution
> ("neg_binomial") be the one from Gelman et al.'s Bayesian Data
> Analysis. Is that what the fix entails?
>
> - Bob
>

yes, it is. I created a pull request explain this.

[Avraham Adler]

So it was a problem of using the reciprocal of beta?

[Marco Inacio]

I think it is.

[Avraham Adler]

Great; I'm glad that has worked out. The negative binomial causes much confusion with its alternate definitions in literature (similar to the generalized Pareto, although there, at least, the two different versions have different numbers of parameters so it's easy to see which is being used). I may take a shot at adding some text to the manual to clarify a few things one of these days. Thanks!

[Linas Mockus]

Below is output.  The rstan I have is 2.1.0.  It looks like that

        y_incorrect [i] <- neg_binomial_rng(shape, shape/mu);

generates incorrect samples.

Linas