How to implement sum of N pareto distributions where N is poisson distributed and its approximation

[Albert Lee]

Hello! 

Thanks for accepting me into the group and also creating a software like stan! 

I am collaborating with a mathematician to fit some complicated biological process and got stuck in the middle of the design.

I read the manual, but couldn't find the answer I wanted, so here I ask. My questions are three fold.

1. I am trying to construct a model involving a distribution characterized as a sum of N-many iid Pareto random variables where N itself is a Poisson distributed. Is such thing possible in stan?

//pseudo code

N ~ poisson(lambda)

for (i in 1:N) { //

    M += pareto(y_min, alpha)

}

That is, can I use for loop based on a parameter? ( I saw Michael's lecture and he sort of mentions something alone the line (http://www.youtube.com/watch?v=xWQpEAyI5s8&t=16m34s) but wanted to ask in the context of my problem). If so how?

2. I learned that there is a approximation for this sort of compound Poisson when lambda is >=100, which is Landau distribution

X ~ Landau  // pdf has the form : 1/pi * integral from 0 to infinity (e ^-(-t*log(t)-xt) * sin(pi * t)) dt)

M = y_min * lambda (X + ln(lambda))

But then again, can Landau distribution be implemented in stan when it involves integral?

3. If #1 and #2 do not work, can I then pass the lookup table of log Landau pdf (generated from mathematica) to stan since its parameter-less?  As far as I know, the only interface to giving a custom distribution to stan is via increment_log_prob, but I was curious if "the passing the lookup table"  is possible.

Any help/advice is appreciated! thank you very much!!

Best,

Albert

[Michael Betancourt]

>
> 1. I am trying to construct a model involving a distribution characterized as a sum of N-many iid Pareto random variables where N itself is a Poisson distributed. Is such thing possible in stan?
>
> //pseudo code
> N ~ poisson(lambda)
> for (i in 1:N) { //
> M += pareto(y_min, alpha)
> }
>
> That is, can I use for loop based on a parameter? ( I saw Michael's lecture and he sort of mentions something alone the line (http://www.youtube.com/watch?v=xWQpEAyI5s8&t=16m34s) but wanted to ask in the context of my problem). If so how?

The problem here is that you want to infer the discrete parameter N, but Stan can’t fit discrete parameters.
Your best bet is to assume a super population — i.e. a large static population, each with a probability of
being active (in which case they add a pareto distribution to the likelihood) or inactive (in which case they
add nothing to the likelihood). For an ecological example of this approach see, for example,
https://github.com/stan-dev/example-models/blob/master/knitr/dorazio-royle-occupancy/dorazio-et-al-knitr.Rmd

> 2. I learned that there is a approximation for this sort of compound Poisson when lambda is >=100, which is Landau distribution
>
> X ~ Landau // pdf has the form : 1/pi * integral from 0 to infinity (e ^-(-t*log(t)-xt) * sin(pi * t)) dt)
> M = y_min * lambda (X + ln(lambda))
>
> But then again, can Landau distribution be implemented in stan when it involves integral?

Not unless you can find some analytic form that approximates it well enough.

> 3. If #1 and #2 do not work, can I then pass the lookup table of log Landau pdf (generated from mathematica) to stan since its parameter-less? As far as I know, the only interface to giving a custom distribution to stan is via increment_log_prob, but I was curious if "the passing the lookup table" is possible.

A lookup table won’t work because it doesn’t provide the associated gradient which is needed for HMC.

[Daniel Rosenbloom]

Hi Michael,

Thanks for your response. I'm working with Albert on this project, and I'm interested in following up on your answer to #3 about the lookup table.

Say that we have a one-dimensional lookup table for this Landau distribution, and for values between the provided values in the table, we are happy with linear interpolation. If the precomputed points are a0, a1, a2...aN then the distribution function could look schematically like this:

if x < a0:

   return (simple function describing left-tail behavior)

if a0 < x < a1:

   return (linear interpolation between f(a0) and f(a1))

if a1 < x < a2:

   return (linear interpolation between f(a1) and f(a2))

...

if x > aN:

   return (simple function describing right-tail behavior)

We have a table with 4000 precomputed points and could automatically generate a very long if statement. Would this not play well with HMC? Of course the derivatives can be computed anywhere between the precomputed points -- is this fact something that Stan effectively "knows," or is there some way to tell it that?

Re: Answer to #1 -- Thanks for the idea! We may end up trying something like this for the case where lambda is less than 100 or so. But lambda can be far larger than that in many cases (say, ~10^9), meaning that the Landau distribution alternative would be required.

Thanks again,

Daniel

[Bob Carpenter]

Stan assumes that the model plus transforms define
a continuously differentiable density. For now, we only
use first-order gradient information.

If there are kinks at the boundaries of the piecewise
function, Stan could have trouble sampling across boundaries
unless step sizes are small enough that the discontinuity in
derivatives doesn't affect the slice rejection step in NUTS.

- 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.

[Daniel Rosenbloom]

That's perfect. Exactly what I needed to know -- so it sounds like if we use a smooth interpolation strategy, we will be safe.

[Joseph Catanzarite]

Hi Daniel

I would like to follow up with you on this post. Were you able to successfully sample from a 1-dimensional lookup table in Stan?
If so, could you please explain how it was implemented?

I also want to be able to sample from posteriors that can be specified in lookup tables (and that are presumably smooth).

Thanks!

Joe Catanzarite

[Daniel Rosenbloom]

Hi Joe,

Yes, this worked -- I posted the answer in response to the new thread you made.

Daniel

On Tue, Sep 13, 2016 at 9:02 PM, Joseph Catanzarite <jcat...@gmail.com> wrote:

Hi Daniel

I wanted to follow up with you on this post. Were you able to successfully sample from a 1-dimensional lookup table in Stan?

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