RNG for non-standard discrete distribution?

Hello, is there an example on how to code up a custom RNG for a non-standard discrete distribution?

The examples in the manual are only for continuous distribution?

I am interested in the Discrete-Weibull distribution.

But any example is welcome and I will try to adapt it for my need. Thank you.

Iā€™m assuming you are interested in discrete user-defined probability functions (and not in the _rng functions)?

Here is a model that fits a discrete weibull distribution (parameterized by \alpha and \beta as in the wikipedia article) to some data y:

functions {
    real discrete_weibull_lpmf(int x, real alpha, real beta) {

        return log_diff_exp(-pow(x * 1.0 / alpha, beta), -pow(((x + 1.0) / alpha), beta));
    }
}

data {
    int N;
    int<lower=0> y[N];
}

parameters {
    real<lower=0> alpha;
    real<lower=0> beta;
}

model {
    alpha ~ normal(0, 5);
    beta ~ normal(0, 5);

    for (i in 1:N) {
        y[i] ~ discrete_weibull(alpha, beta);
    }
}

log_diff_exp(a, b) is the same as log(exp(a) - exp(b)), just a bit more numerically stable.

In case you are actually interested in a discrete_weibull_rng function, a common trick is to draw from a standard uniform distribution and transform these by the inverse cdf.

Thanks. I actually AM interested in a discrete_weibull_rng function!

Once I solve for x using the inverse cdf, do I use the floor function to ensure it is an integer?

That would be my first attempt. You can always produce samples using your scheme and then simulate a large number of random draws and compare the frequencies with the true probabilities.

Note that all rounding functions return a real, so your random draws are always represented as 2.0, 3.0 etc. There is no good workaround to convert them to integers right now. Depending on your use case the easiest solution is probably to do the sampling outside of Stan.

2 Likes