Re: [Boost-users] [distributions]: Inverse Gamma

Subject: Re: [Boost-users] [distributions]: Inverse Gamma
From: Paul A. Bristow (pbristow_at_[hidden])
Date: 2010-08-04 13:24:19

> -----Original Message-----
> From: boost-users-bounces_at_[hidden] [mailto:boost-users-bounces_at_[hidden]] On Behalf Of Thomas Mang
> Sent: Wednesday, August 04, 2010 4:02 PM
> To: boost-users_at_[hidden]
> Subject: Re: [Boost-users] [distributions]: Inverse Gamma
> 
> On 04.08.2010 10:37, Paul A. Bristow wrote:
> >
> >
> >> -----Original Message-----
> >> From: boost-users-bounces_at_[hidden] [mailto:boost-users-bounces_at_[hidden]] On Behalf Of Thomas Mang
> >> Sent: Tuesday, August 03, 2010 3:07 PM
> >> To: boost-users_at_[hidden]
> >> Subject: Re: [Boost-users] [distributions]: Inverse Gamma
> >>
> >   >>
> >>>> any plans of implementing the inverse gamma distribution as part of the distributions library ?
> >>>
> >>> This looks possible - but I'm curious about applications - you obviously have one, but Wikipedia doesn't mention
any
> >>>
> >>> http://en.wikipedia.org/wiki/Inverse-gamma_distribution
> >>>
> >>> But you obviously have one ;-)
> >>
> >> Yes I truly have one ;) The inverse gamma distribution and its special
> >> case, the scaled inverse chi-square distribution, is the conjugate prior
> >> to the normal distribution variance parameter in Bayesian statistics.
> >> Pretty much as uncommon and unheard of as it is outside Bayes world [to
> >> the best of my knowledge], it's very much central to Bayesian stats and
> >> appears in every textbook right after the introduction chapter ;)
> >>
> >> http://en.wikipedia.org/wiki/Scaled_inverse_chi-square_distribution
> >> http://en.wikipedia.org/wiki/Conjugate_prior
> >>
> >> Hence I wonder it has not been requested so far - but being a Bayesian
> >> C++ / booster I definitely want / need it :).
> >>
> >> @John: Yes it is a transformation deviate of the gamma, and an easy so.
> >> And it should be fairly easy to implement IMHO.
> >
> > The pdf and pdf etc looks fairly straightforward (only uses exp, pow and gamma?) so I might be persuaded to do
these.
> > But the inverses (qhantiles) may prove more troublesome if have to be done by brute force numerically
> >
> > R library does it numerically http://rss.acs.unt.edu/Rdoc/library/pscl/html/igamma.html
> >
> > - or are there analytic expressions for these?
> >
> > Or can it use the inverse of the gamma distribution?
> 
> 
> Yes pdf just uses exp, pow and gamma function, for the cdf /
> complementary cdf also the incomplete gamma function is needed - all in
> the library available.
> For quantiles yes you can use quantiles of the gamma distribution.
> Indeed all of pdf, cdf , complementary cdf and quantiles can be derived
> using calls to the corresponding functions of the gamma distribution:
> 
> If we define both the gamma and inverse gamma distribution with the
> shape and scale parameter (as is done for the boost implementation of
> the gamma) and in that constructor order, then the relations are [should
> be] as follows:
> 
> typedef gamma_distribution<...> gamma_dist;
> typedef inverse_gamma_distribution<...> inv_gamma_dist;
> 
> pdf(inv_gamma_dist(a, b), x) == pdf(gamma_dist(a, 1/b), 1/x) * 1/(x*x);
> 
> cdf(inv_gamma_dist(a, b), x) == 1 - cdf(gamma_dist(a, 1/b), 1/x);
> 
> quantile(inv_gamma(a, b), p) == 1 / quantile(gamma_dist(a, 1/b), 1/x);
> 
> 
> where for the pdf the extra term 1/(x*x) is the density transformation
> Jacobian term while for the probabilities no equivalent term is needed.
> Note that I just wrote "1" as nominator to save space and avoid a
> newsgroup-linebreak; it should be a floating point type in real code of
> course. And my "[should be]" refers to that I just worked that out and
> tested it in R, but not in C++. But it will be double-checked anyway
> after it's implemented. Why R evaluates the quantile numerically is
> beyond my understanding.
> Whether you implement the pdf / cdf directly or go for above relations
> is probably an implementation issue, in terms of performance it
> shouldn't make a big difference.
> 
> If you go for implementation that would be great and very much
> appreciated with lots of thanks. Could you then also implement the
> inverse chi-square and scaled inverse chi-square distribution as well
> (they are just special inverse gamma distributions, but their names
> might ring more bells and I'd find it appealing to offer them straight
> away, similar to Chi-square and Gamma distribution cases) ?

OK that seems to be what I would need to know.  I'll give this a try soon and tell you when a draft is available.

> > (And there is the question of testing - some parameters and value combinations (preferably exact) are needed for
sanity
> > and accuracy checks).

Any *exact* pdf and cdf values from *exact* parameters would be helpful - if you know any.  But I'll look more closely
at the equations to try to find some too.

Do you have easy access to R that you could produce some test values (to double precision only of course)?

Paul

---
Paul A. Bristow
Prizet Farmhouse
Kendal, UK   LA8 8AB
[phone], mobile [phone]
pbristow_at_[hidden]