Sampling from a Polya Gamma distribution?

Hello everyone.
Anybody tried to implement the Polya-Gamma distribution?
Or is there any plan for supporting the Polya-Gamma distribution?


This paper claims that sampling from PG distribution is beneficial for logistic models.
Wonder if anybody tried this on Stan.

This is quite useful for using a Gibbs sampler to estimate the coefficients of logit models using data augmentation, and in fact, it can be proven that this particular Gibbs sampler is uniformally ergodic, which is rare. But the approach isn’t that useful for Stan, which uses No U-turn Sampling and joint proposals rather than full-conditional proposals and usually samples from the posterior distribution of the coefficients of a logit model (without data augmentation) quite efficiently as well.

The Polya Gamma distribution might be useful for other purposes, although I haven’t seen any and I don’t know of anyone who has tried to implement it in Stan in part because that infinite sum with alternating signs would be challenging to do in a numerically stable way.

1 Like

Thanks for the answer. I was trying to implement [1, 2] which they used PG for modeling the likelihood of stick-breaking processes. Maybe it’s not beneficial in this case too?

[1] Poddar, Lahari, Wynne Hsu, and Mong Li Lee. “Quantifying aspect bias in ordinal ratings using a bayesian approach.” arXiv preprint arXiv:1705.05098 (2017).

[2] Chen, Jianfei, et al. “Scalable inference for logistic-normal topic models.” Advances in Neural Information Processing Systems . 2013.

Continuing the discussion from Sampling from a Polya Gamma distribution?:

Hello,I’m also interested in this problem. And I wonder if you have solved the problem? Did you implement the code in the articles

As far as I know there arent any publicly available implementations for Stan.

Not sure if this helps and apologize if you already aware of this but the BayesLogit R packages has the implementation of PolyaGamma distribution. Its done in C so its faster than if you would reimplement PG in R but its really not super efficient.

3+ years ago me and Erik rewrote their implementation in C and then also made a GPU variant of it that is roughly a 100x faster altogether. Its also in a CRAN R package tho that is not being actively maintained and not sure if I will ever come back to that (https://cran.r-project.org/web/packages/bayesCL/index.html).

But if anyone wants to rewrite my C code in C++ for Stan Math I am happy to provide more info or help on that.

1 Like

@ergou_chen Yes I implemented a similar model on Stan.
Instead of using the Polya-gamma sampling scheme,
I simply typed in the model and ran the default NUTS sampler and it just worked.