few desired features for Rstan

[Linas Mockus]

Hi,

I would like to add some simple features to some future release of Rstan (there are workarounds of course).

1. I have been struggling with neg_binomial_2_log_rng for a long time.  While the chains converge during warmup the error is thrown up because the parameter for the underlying poisson_rng is too large.
Ben suggested the following excellent workaround:

  if (is_inf(theta))
    gamma_temp <- nu;
  else
    gamma_temp <- gamma_rng(theta, theta / nu);
  if (gamma_temp < poisson_max)
    neg_bin_rng <- poisson_rng(gamma_temp);
  else
    neg_bin_rng <- normal_rng(gamma_temp, sqrt(gamma_temp));
Another alternative is to allow generated quantities block to be executed only during sampling.

2. It is nice to see the progress of parallel chains not only in Internet Explorer but also as text file. Below is the workaround I came with by using previous postings. The output of all chains is redirected to one outfile.

par.stan <- function (file, code, data, warmup, n.draws, chains, pars,
    init, outfile)
{
    # compile C++ code
    sm <- stan (file = file, model_code = code, data = data, chains = 1,
        iter = 0)
    CL <- makeCluster (chains, outfile = outfile)
    clusterExport (cl = CL, c ("data", "pars", "warmup", "n.draws"))
    sflist <- parLapply (CL, 1 : chains, fun = function (cid) {
        require (rstan)
        stan (fit = sm, data = data, chains = 1, warmup = warmup,
            iter = n.draws + warmup, pars = pars, init = init,
            chain_id = cid)
    })
    stopCluster (CL)
    if (chains == 1) {
        return (sflist [[1]])
    } else {
        return (sflist2stanfit (sflist))
    }
}
stanmod <- function (chain, file, sm, data, warmup, n.draws, pars, init) {
    if (length (init) == 1) {
        if (init == "random") {
            init <- "random"
        } else {
            init <- list (init [[chain]])
        }
    } else {
        init <- list (init [[chain]])
    }
    sampling (sm, data = data, pars = pars, chains = 1,
        iter = warmup + n.draws, warmup = warmup,
        init = init, 
#        diagnostic_file = diagnostic_file,
#        sample_file = sample_file,
        verbose = TRUE, refresh = max ((warmup + n.draws) / 10, 1),
        chain_id = chain)
}

Thanks,
Linas

[Bob Carpenter]

The thing to do is create feature-request issues in
GitHub. Project (1) is for Stan, project (2) for RStan.
They have different GitHub repos and issue trackers.

1. I'm not sure what the feature is here. A new function?
A rewrite of the existing neg_bin_rng?

Assuming that neg_bin_rng is the result (I wouldn't call it that,
since I take an RNG to be a random number generator, not
a random draw), in one case it's an integer and in the
other it's real. Is there some kind of truncation?

Stan's data type only goes up to int for integers, so
that maxes out at 2^31. We put our bound where we did
so as to be able to handle the whole of range of
poisson_rng.

To avoid that, you can do the test yourself for values
that are too large. Or you can sample outside of Stan.
If the mean's too large, then you wind up expecting draws
that are too large to fit in an int anyway.

Do you want to assign a real value to the result in a
different function?

Extend Stan's int to long? This would be great, but it's
a big project.


2. I'll leave that to the folks who know about R and/or Windows.

- Bob

> --
> 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.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Linas Mockus]

Thanks Bob.

My guess is that rewriting neg_binomial_rng and neg_binomial_2_rng
according to Ben's workaround will make those functions more usable.

Thank you,
Linas

--
Thank you,
Linas

[Bob Carpenter]

Ben's function returns real values (normal_rng), whereas neg_binomial_2 and
poisson return integers. In Stan, those can only go up to 2^31
on the positive side, so anything much outside the bounds we
have set now can't generate values we can represent in Stan's precision.

Once we're up in the 2^30 (1 billion) range, the integer rounding's probably
not going to matter much, but we can't handle the overflow.

It would totally make sense to extend Stan's integer type to long.
Then we could just extend the precision on our discrete RNGs.

- Bob

[Linas Mockus]

Thanks Bob. Long would be ideal.  Otherwise it really doesn't matter since overflow in generated block occurs during warm up and is not used for predictive posterior.  However it stops the further processing for a chain.  My workaround is to use neg_binomial_2_log_rng as is but run 4 chains. Typically overflow occurs in the first, sometimes in both the first and the second chains.  I am still arriving to a good converged result (as shown by trace).  I don't want to change the code in hopes that neg_binomial_2_log_rng will be eventually fixed.

Thanks, Linas

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/rgpHgBtFS3k/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Bob Carpenter]

I would suggest using a conditional in the generated
quantities:

if (mean > threshold)
sample[n] <- -1;
else
draw sample

That way, you know where it fails in the posterior, and
if it's just during warmup, you won't get -1 values.

- Bob

[Linas Mockus]

Bob,

In order to prevent work-arounds my suggestion is to rewrite
neg_binomial_2_rng(mu,psi) as suggested by Ben:
real gamma_temp;
if (is_inf(psi))
gamma_temp <- mu;
else
gamma_temp <- gamma_rng(psi, psi / mu);
if (gamma_temp < pow(2.0, 30.0))
result <- poisson_rng(gamma_temp);
else
result <- normal_rng(gamma_temp, sqrt(gamma_temp));

Thank you,
Linas
On Sat, 03 Oct 2015 16:16:04 -0400, Bob Carpenter <ca...@alias-i.com> wrote:

--
Thank you,
Linas

[Bob Carpenter]

It won't work. poisson_rng returns an integer
and normal_rng returns a real value. You could
try casting this back to an int, but that won't work
because the whole problem originates because the
result won't fit into 32-bit integers.

The workable solution is to change Stan's integer type
from 32-bit to 64-bit.

- Bob

[Linas Mockus]

Bob,

Actually my problem went away as long as I declared result as real type.  So this is not a theoretical neg_binomial but some reasonable approximation.  I also prefer integers to become 64-bit so work arounds are not needed.

Thanks,

Linas

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/rgpHgBtFS3k/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Bob Carpenter]

None of these changes are R specific. You could
create an issue in stan-dev/math for the function, then
it could be run through Stan itself in stan-dev/stan.
There's nothing to change in RStan itself.

But really we need to get Stan to 64-bits. Then we won't
have trouble until that mean gets near 2^60 (roughly 10^18).

- Bob