Hi George,
The inverse gamma 2 density appears naturally in Bayesian statistics when considering conjugate priors in linear Gaussian models. The inverse gamma densities are quite common in Bayesian statistics (the references you are giving are not pertinent with respect to this literature).
Furthermore the inverse gamma 2 is just a re-parametrization of the probability textbook (or common if you want) inverse gamma density. This is shown in the attached pdf. So I don't think we should provide the inverse gamma prior you are suggesting in Dynare.
Best, Stéphane.
On 2 June 2010 10:35, George Perendia george@perendia.orangehome.co.ukwrote:
Hi Michel
Thanks for the pdf document .
- My apologies - Having had a chance to see the full representation of
both inv gamma, 1 and 2 pdf in the document you sent I do not anymore think there is a bug in the implementation - the implementation seem to conform closely to the given representations.
- I can not tell if the notes are incorrect but they do seem mutually
consistent, so that e.g.
lpdfig1(x,s,n) = lpdfig2(x*x,s,n) = lpdfgam(1/(x*x),n/2,2/s) - 2*log(x*x) = 0.5673 for my test
On the other hand, both inv gamma 1 and 2 representations differ from what I find the most common definition of the inverse gamma from various sources
invGammaPDF(x,a,b) = b^a / gamma(a) * 1/x ^(a+1) * exp( -b/x) which, in turn = GammaPDF(1/x,a,1/b) / x^2
(though, alternative representations of gammaPDF with inverse scale beta = 1/b would give = GammaPDF(1/x,a,beta) / x^2 )
(e.g. see the enclosed document or internet reference such as http://www.chemistrydaily.com/chemistry/Inverse-gamma_distribution or even wikipedia http://en.wikipedia.org/wiki/Inverse-gamma_distribution )
We may need to add another inv Gamma for that common representation.
Best George
On 01/06/2010 20:37, Michel Juillard wrote:
Hi George,
inverse gamma distributions are tricky because different authors use different notations. I attach the notes that I had made in 2004. Could you please tell me
- if you think that these notes are erroneous. In this case, please,
indicate the sources that you use. 2) if the code currently implemented in Dynare Matlab conforms with these notes.
I'm leaving tomorrow for the Dynare Conference. I will be back this week end and will try to look into the issue next week.
All the best,
Michel
On 06/01/2010 09:00 PM, George Perendia wrote:
PS: On the other hand, using the commonly used inv gamma formula
pdgInvGamma (x,n,s) = 1/(gamma(n) * s^n) * x ^(-n-1) exp (1/ (x*s) )
it appear that
pdgInvGamma (x,n,s) = pdfGamma (1/x, n, s) / x^2
or, in logs,
lpdfinvg (x,n,s) = -gammaln(n) - n*log(s) + (n+1)*log(1/x) - (1/x)/s = [ -gammaln(n) - n*log(s) + (n-1)*log(1/x) - (1/x)/s ] - 2*log(x) = [ lpdfgam(1/x,n,s) ] - 2*log(x)
Best regards
George
On 01/06/2010 16:24, George Perendia wrote:
Dear Michel
thanks
I have 2 issues:
- According to few sources,
If x~InvGamma(a,b) , then 1/x ~Gamma(a,1/b) distribution however, in our dynare inv gamma 2 (lpdfig2.m) we have different definition: X ~ IG2(s,nu) if X = inv(Z) where Z ~ G(nu/2,2/s) (Gamma distribution) (which is equivalent to If x~InvGamma(a,b) , then 1/x ~Gamma(b/2,2/a) distribution)
- I suspect there are a bugs in implementation for both invg2 and invg1
even as per our specifications of invg2 aind invg1, i.e. for invg2: we have
ldens(idx) = -gammaln(.5*nu) - (.5*nu)*(log(2)-log(s)) -.5*(nu+2)*log(x(idx)) -.5*s./x(idx);
and for e.g. s=2, n=5 and x= 0.5 as a result we have (incorrectly?):
lpdfig2(x,s,n) = -gammaln(.5*n) - (.5*n)*(log(2)-log(s)) - .5*(n+2)*log(x) -.5*s./x
ans =
0.1413whilst
lpdfgam(1/x,n/2,2/s)
ans =
-1.2450
and modified (corrected?) invgam2 with ...(nu - 2).which meets our spec.. instead of (nu + 2)...
K>> -gammaln(.5*n) - (.5*n)*(log(2)-log(s)) - .5*(n-2)*log(x) -.5*s./x
ans =
-1.2450
thus, gives same result as lpdfgam(1/x,n/2,2/s).
Similarly, invg1
ldens(idx) = log(2) - gammaln(.5*nu) - .5*nu*(log(2)-log(s)) -(nu+1)*log(x(idx)) - .5*s./(x(idx).*x(idx)) ;
lpdfig1(x,s,n)= 0.5673 = lpdfig2(x*x,s,n)= 0.5673
whilst lpdfgam(1/(x*x),n/2,2/s)
ans =
-2.2052
and modified invg1 without leading log(2) and with nu-2 instead nu+1
- gammaln(.5*nu) - .5*nu*(log(2)-log(s)) - (nu-2)*log(x) - .5*s./(x.*x)
ans =
-2.2052
i.e. = lpdfgam(1/(x*x),n/2,2/s) as per definition
(... though I may have got something wrong)
If they are bugs, shall I log them and correct them?
Best george
On 19/05/2010 20:29, Michel Juillard wrote:
Dear George,
thank you for the message. I agree that we should go with BOOST.
- In the estimation DLL, we only need a uniform and a normal random
generator and the pdfs for the normal, beta, gamma, inverse gamma 1, inverse gamma 2 and the uniform. You shouldn't try to implement more than that at this stage.
- for the inverse gamma 1 and 2 pdfs, you can directly translate the
current matlab function in ./matlab/distributions. We should add the source of the algorithm if we can remember.
- You need to look at the Matlab implementation, because we need the same
kind of generalizations:
- gamma is defined as greater than p3 rather than greater than 0
- beta is defined over [p3; p4] instead of over [0,1]
The relevant Matlab functions are set_prior.m and priordens.m
All the best,
Michel
On 05/19/2010 06:42 PM, George Perendia wrote:
Hi Sebastien
I started looking into priordens and pdfs as well as the relvant random number generators (rng-s) for estimated parameters (*) and investigated two options: C gsl and the C++ header BOOST as initially suggested by Stephane.
Both provide for the most of the requirements but there are seeral insufficiencies.
gsl provides the necessary pair of pdf and random number generators functions based on given uniform random number but not the uniform random number generator - rng (or the pdf but that is trivial to add).
On the other hand whilst BOOST provides a large variety of functions including the pdfs and cdfs in a template form as well as a variety of uniform and other distributions rngs and ways to combine them. Though it does not seem to provide a rng for beta, the beta rng may be easily derived using gamma rng and, e.g. the relation:
"If *X* and *Y* are independently distributed Gammahttp://wiki/Gamma_distribution (*α*, *θ*) and Gamma(*β*, *θ*) respectively, then *X* / (*X* + *Y*) is distributed Beta(*α*, *β*), "
...which is how gsl beta rng is implemented anyway:
Neither gsl nor Boost however seem to provide pdfs or rngs for inverse_gamma and that would need to be derived for either of them too.
Whilst the gsl is simpler it is C whlst the BOOST is C++. Also, gsl would require a library compilation before it can be used and it would also be an additional package requirement for building Dynare (i.e. whilst we are already using BOOST header library anyway).
It seems that BOOST pdfs and rngs are better documented, provide a larger variety of options for later extensions and fit better into our design strategy whilst, as part of already used BOOST, they do not require any extension on the overall maintenance and build requirements of Dynare.
I wold be inclined therefore to suggest us to use BOOST but please do let me know if you have any objections or other suggestions.
If the performance is the issue, we can run a comparative test e.g. for throwing 100000 beta rng-s and deriving their pdf-s.using gsl and BOOSt and see the difference but that may a be bit excessive - just use of a different uniform rng for generating specific distribution rng from a variety of those available in BOOST can change performance dramatically.
*) assuming the following distributions are required beta_pdf | gamma_pdf | normal_pdf | uniform_pdf | inv_gamma_pdf | inv_gamma1_pdf | inv_gamma2_pdf
Best regards George
On 18/05/2010 17:26, Sébastien Villemot wrote:
Hi George,
On Wed, May 12, 2010 at 09:36:26AM +0100, George Perendia wrote:
It was good to work together last week. Aftert yesterday's ACES meeting it looks like I may not be engaged full time on ACES for now, so, I could contiue on Estimation DLL, 1-2 days a week for time being, e.g. on prior dens, optmiiser integration, MCMC or something else. Let me know what wold be best for you.
Sorry for my late answer.
I am also happy with the progress that we made when you came to Paris.
Unfortunately, since I have been working on other projects, I have not been able to finalize the likelihood DLL yet, but I am quite close to it.
I am leaving for vacations until May 30, and I hope to finish the DLL by the first days of June. Note that I will still read and answer to emails during my vacations.
Please feel free to advance on the implementation of other parts of the design. In particular, we need to decide how to implement the prior densities. We should rely on existing code rather than reimplementing all the probability density functions. My feeling is that using the GNU Scientific Library for that purpose is the easiest way, since all the PDFs are there and it is licensed under the GPL. Do you have any other option in mind ?
Best,
--
--
--
--
Best regards George 07951415480
Hi Stephane
Thanks for your very detailed notes. I am afraid that although they are consistent within themselves and your comment that the inv gamma2 is just reparametrised inv gamma, they however appear to be inconsistent with both Michel's priors.pdf and the current implementation where the inv gamma 2 (and 1) modified parameters appear to be permuted when applied on gamma in comparison to your notes. I.e.:
% X ~ IG2(s,nu) if X = inv(Z) where Z ~ G(nu/2,2/s) (Gamma distribution) % X ~ IG1(s,nu) if X = sqrt(Y) where Y ~ IG2(s,nu) and Y = inv(Z) with Z ~ G(nu/2,2/s) (Gamma distribution)
and, as result, e.g.
lpdfig1(x,s,n) = lpdfig2(x*x,s,n) = lpdfgam(1/(x*x),n/2,2/s) - 2*log(x*x) = 0.5673 for my test (see below)
whilst , permuting s and n as in your equation usually gives a different result, e.g.
lpdfgam(1/(x*x),s/2,2/n) - 2*log(x*x) = -6.3111
It thus appear to me that the Michel's priors.pdf and our current implementation either reflect a different notion of inv gamma 1 and 2, or, our inv gamma 1 and 2 gamma-shape and gamma-scale parameters are somehow (or somewhere) permuted ?
best george
On 02/06/2010 15:57, Stéphane Adjemian wrote:
Hi George,
The inverse gamma 2 density appears naturally in Bayesian statistics when considering conjugate priors in linear Gaussian models. The inverse gamma densities are quite common in Bayesian statistics (the references you are giving are not pertinent with respect to this literature).
Furthermore the inverse gamma 2 is just a re-parametrization of the probability textbook (or common if you want) inverse gamma density. This is shown in the attached pdf. So I don't think we should provide the inverse gamma prior you are suggesting in Dynare.
Best, Stéphane.
On 2 June 2010 10:35, George Perendia <george@perendia.orangehome.co.uk mailto:george@perendia.orangehome.co.uk> wrote:
Hi Michel Thanks for the pdf document . 2) My apologies - Having had a chance to see the full representation of both inv gamma, 1 and 2 pdf in the document you sent I do not anymore think there is a bug in the implementation - the implementation seem to conform closely to the given representations. 1) I can not tell if the notes are incorrect but they do seem mutually consistent, so that e.g. lpdfig1(x,s,n) = lpdfig2(x*x,s,n) = lpdfgam(1/(x*x),n/2,2/s) - 2*log(x*x) = 0.5673 for my test On the other hand, both inv gamma 1 and 2 representations differ from what I find the most common definition of the inverse gamma from various sources invGammaPDF(x,a,b) = b^a / gamma(a) * 1/x ^(a+1) * exp( -b/x) which, in turn = GammaPDF(1/x,a,1/b) / x^2 (though, alternative representations of gammaPDF with inverse scale beta = 1/b would give = GammaPDF(1/x,a,beta) / x^2 ) (e.g. see the enclosed document or internet reference such as http://www.chemistrydaily.com/chemistry/Inverse-gamma_distribution or even wikipedia http://en.wikipedia.org/wiki/Inverse-gamma_distribution ) We may need to add another inv Gamma for that common representation. Best George On 01/06/2010 20:37, Michel Juillard wrote:Hi George, inverse gamma distributions are tricky because different authors use different notations. I attach the notes that I had made in 2004. Could you please tell me 1) if you think that these notes are erroneous. In this case, please, indicate the sources that you use. 2) if the code currently implemented in Dynare Matlab conforms with these notes. I'm leaving tomorrow for the Dynare Conference. I will be back this week end and will try to look into the issue next week. All the best, Michel On 06/01/2010 09:00 PM, George Perendia wrote:PS: On the other hand, using the commonly used inv gamma formula pdgInvGamma (x,n,s) = 1/(gamma(n) * s^n) * x ^(-n-1) exp (1/ (x*s) ) it appear that pdgInvGamma (x,n,s) = pdfGamma (1/x, n, s) / x^2 or, in logs, lpdfinvg (x,n,s) = -gammaln(n) - n*log(s) + (n+1)*log(1/x) - (1/x)/s = [ -gammaln(n) - n*log(s) + (n-1)*log(1/x) - (1/x)/s ] - 2*log(x) = [ lpdfgam(1/x,n,s) ] - 2*log(x) Best regards George On 01/06/2010 16:24, George Perendia wrote:Dear Michel thanks I have 2 issues: 1) According to few sources, If x~InvGamma(a,b) , then 1/x ~Gamma(a,1/b) distribution however, in our dynare inv gamma 2 (lpdfig2.m) we have different definition: X ~ IG2(s,nu) if X = inv(Z) where Z ~ G(nu/2,2/s) (Gamma distribution) (which is equivalent to If x~InvGamma(a,b) , then 1/x ~Gamma(b/2,2/a) distribution) 2) I suspect there are a bugs in implementation for both invg2 and invg1 even as per our specifications of invg2 aind invg1, i.e. for invg2: we have ldens(idx) = -gammaln(.5*nu) - (.5*nu)*(log(2)-log(s)) - .5*(nu+2)*log(x(idx)) -.5*s./x(idx); and for e.g. s=2, n=5 and x= 0.5 as a result we have (incorrectly?): lpdfig2(x,s,n) = -gammaln(.5*n) - (.5*n)*(log(2)-log(s)) - .5*(n+2)*log(x) -.5*s./x ans = 0.1413 whilst lpdfgam(1/x,n/2,2/s) ans = -1.2450 and modified (corrected?) invgam2 with ...(nu - 2).which meets our spec.. instead of (nu + 2)... K>> -gammaln(.5*n) - (.5*n)*(log(2)-log(s)) - .5*(n-2)*log(x) -.5*s./x ans = -1.2450 thus, gives same result as lpdfgam(1/x,n/2,2/s). Similarly, invg1 ldens(idx) = log(2) - gammaln(.5*nu) - .5*nu*(log(2)-log(s)) - (nu+1)*log(x(idx)) - .5*s./(x(idx).*x(idx)) ; lpdfig1(x,s,n)= 0.5673 = lpdfig2(x*x,s,n)= 0.5673 whilst lpdfgam(1/(x*x),n/2,2/s) ans = -2.2052 and modified invg1 without leading log(2) and with nu-2 instead nu+1 - gammaln(.5*nu) - .5*nu*(log(2)-log(s)) - (nu-2)*log(x) - .5*s./(x.*x) ans = -2.2052 i.e. = lpdfgam(1/(x*x),n/2,2/s) as per definition (... though I may have got something wrong) If they are bugs, shall I log them and correct them? Best george On 19/05/2010 20:29, Michel Juillard wrote:Dear George, thank you for the message. I agree that we should go with BOOST. 1) In the estimation DLL, we only need a uniform and a normal random generator and the pdfs for the normal, beta, gamma, inverse gamma 1, inverse gamma 2 and the uniform. You shouldn't try to implement more than that at this stage. 2) for the inverse gamma 1 and 2 pdfs, you can directly translate the current matlab function in ./matlab/distributions. We should add the source of the algorithm if we can remember. 3) You need to look at the Matlab implementation, because we need the same kind of generalizations: - gamma is defined as greater than p3 rather than greater than 0 - beta is defined over [p3; p4] instead of over [0,1] The relevant Matlab functions are set_prior.m and priordens.m All the best, Michel On 05/19/2010 06:42 PM, George Perendia wrote:Hi Sebastien I started looking into priordens and pdfs as well as the relvant random number generators (rng-s) for estimated parameters (*) and investigated two options: C gsl and the C++ header BOOST as initially suggested by Stephane. Both provide for the most of the requirements but there are seeral insufficiencies. gsl provides the necessary pair of pdf and random number generators functions based on given uniform random number but not the uniform random number generator - rng (or the pdf but that is trivial to add). On the other hand whilst BOOST provides a large variety of functions including the pdfs and cdfs in a template form as well as a variety of uniform and other distributions rngs and ways to combine them. Though it does not seem to provide a rng for beta, the beta rng may be easily derived using gamma rng and, e.g. the relation: "If /X/ and /Y/ are independently distributed Gamma <http://wiki/Gamma_distribution>(/α/, /θ/) and Gamma(/β/, /θ/) respectively, then /X/ / (/X/ + /Y/) is distributed Beta(/α/, /β/), " ...which is how gsl beta rng is implemented anyway: Neither gsl nor Boost however seem to provide pdfs or rngs for inverse_gamma and that would need to be derived for either of them too. Whilst the gsl is simpler it is C whlst the BOOST is C++. Also, gsl would require a library compilation before it can be used and it would also be an additional package requirement for building Dynare (i.e. whilst we are already using BOOST header library anyway). It seems that BOOST pdfs and rngs are better documented, provide a larger variety of options for later extensions and fit better into our design strategy whilst, as part of already used BOOST, they do not require any extension on the overall maintenance and build requirements of Dynare. I wold be inclined therefore to suggest us to use BOOST but please do let me know if you have any objections or other suggestions. If the performance is the issue, we can run a comparative test e.g. for throwing 100000 beta rng-s and deriving their pdf-s.using gsl and BOOSt and see the difference but that may a be bit excessive - just use of a different uniform rng for generating specific distribution rng from a variety of those available in BOOST can change performance dramatically. *) assuming the following distributions are required |beta_pdf| | |gamma_pdf| | |normal_pdf| | |uniform_pdf| | |inv_gamma_pdf| | |inv_gamma1_pdf| | |inv_gamma2_pdf| Best regards George On 18/05/2010 17:26, Sébastien Villemot wrote:> Hi George, > > On Wed, May 12, 2010 at 09:36:26AM +0100, George Perendia wrote: > >> It was good to work together last week. >> Aftert yesterday's ACES meeting it looks like I may not be engaged >> full time on ACES for now, so, I could contiue on Estimation DLL, >> 1-2 days a week for time being, e.g. on prior dens, optmiiser >> integration, MCMC or something else. >> Let me know what wold be best for you. >> > Sorry for my late answer. > > I am also happy with the progress that we made when you came to Paris. > > Unfortunately, since I have been working on other projects, I have not been able to finalize the likelihood DLL yet, but I am quite close to it. > > I am leaving for vacations until May 30, and I hope to finish the DLL by the first days of June. Note that I will still read and answer to emails during my vacations. > > Please feel free to advance on the implementation of other parts of the design. In particular, we need to decide how to implement the prior densities. We should rely on existing code rather than reimplementing all the probability density functions. My feeling is that using the GNU Scientific Library for that purpose is the easiest way, since all the PDFs are there and it is licensed under the GPL. Do you have any other option in mind ? > > Best, > >
-------- Best regards George 07951415480-- Stéphane Adjemian CEPREMAP & Université du Maine
Tel: (33)1-43-13-62-39
Hi George,
I am not sure to understand your problem... Comparing Michel notes and the codes with my notes I can't see any difference... Except the order of declaration of the parameters nu and s in IG_2 and IG_1. I In my pdf I write IG_2(nu,s) but in the codes we wrote IG_2(s,nu) (as in Michel's notes). I will change these notations in my pdf.
Best, Stéphane.
On 3 June 2010 11:49, George Perendia george@perendia.orangehome.co.ukwrote:
Hi Stephane
Thanks for your very detailed notes. I am afraid that although they are consistent within themselves and your comment that the inv gamma2 is just reparametrised inv gamma, they however appear to be inconsistent with both Michel's priors.pdf and the current implementation where the inv gamma 2 (and 1) modified parameters appear to be permuted when applied on gamma in comparison to your notes. I.e.:
% X ~ IG2(s,nu) if X = inv(Z) where Z ~ G(nu/2,2/s) (Gamma distribution) % X ~ IG1(s,nu) if X = sqrt(Y) where Y ~ IG2(s,nu) and Y = inv(Z) with Z ~ G(nu/2,2/s) (Gamma distribution)
and, as result, e.g.
lpdfig1(x,s,n) = lpdfig2(x*x,s,n) = lpdfgam(1/(x*x),n/2,2/s) - 2*log(x*x) = 0.5673 for my test (see below)
whilst , permuting s and n as in your equation usually gives a different result, e.g.
lpdfgam(1/(x*x),s/2,2/n) - 2*log(x*x) = -6.3111
It thus appear to me that the Michel's priors.pdf and our current implementation either reflect a different notion of inv gamma 1 and 2, or, our inv gamma 1 and 2 gamma-shape and gamma-scale parameters are somehow (or somewhere) permuted ?
best george
On 02/06/2010 15:57, Stéphane Adjemian wrote:
Hi George,
The inverse gamma 2 density appears naturally in Bayesian statistics when considering conjugate priors in linear Gaussian models. The inverse gamma densities are quite common in Bayesian statistics (the references you are giving are not pertinent with respect to this literature).
Furthermore the inverse gamma 2 is just a re-parametrization of the probability textbook (or common if you want) inverse gamma density. This is shown in the attached pdf. So I don't think we should provide the inverse gamma prior you are suggesting in Dynare.
Best, Stéphane.
On 2 June 2010 10:35, George Perendia george@perendia.orangehome.co.ukwrote:
Hi Michel
Thanks for the pdf document .
- My apologies - Having had a chance to see the full representation of
both inv gamma, 1 and 2 pdf in the document you sent I do not anymore think there is a bug in the implementation - the implementation seem to conform closely to the given representations.
- I can not tell if the notes are incorrect but they do seem mutually
consistent, so that e.g.
lpdfig1(x,s,n) = lpdfig2(x*x,s,n) = lpdfgam(1/(x*x),n/2,2/s) - 2*log(x*x) = 0.5673 for my test
On the other hand, both inv gamma 1 and 2 representations differ from what I find the most common definition of the inverse gamma from various sources
invGammaPDF(x,a,b) = b^a / gamma(a) * 1/x ^(a+1) * exp( -b/x) which, in turn = GammaPDF(1/x,a,1/b) / x^2
(though, alternative representations of gammaPDF with inverse scale beta = 1/b would give = GammaPDF(1/x,a,beta) / x^2 )
(e.g. see the enclosed document or internet reference such as http://www.chemistrydaily.com/chemistry/Inverse-gamma_distribution or even wikipedia http://en.wikipedia.org/wiki/Inverse-gamma_distribution )
We may need to add another inv Gamma for that common representation.
Best George
On 01/06/2010 20:37, Michel Juillard wrote:
Hi George,
inverse gamma distributions are tricky because different authors use different notations. I attach the notes that I had made in 2004. Could you please tell me
- if you think that these notes are erroneous. In this case, please,
indicate the sources that you use. 2) if the code currently implemented in Dynare Matlab conforms with these notes.
I'm leaving tomorrow for the Dynare Conference. I will be back this week end and will try to look into the issue next week.
All the best,
Michel
On 06/01/2010 09:00 PM, George Perendia wrote:
PS: On the other hand, using the commonly used inv gamma formula
pdgInvGamma (x,n,s) = 1/(gamma(n) * s^n) * x ^(-n-1) exp (1/ (x*s) )
it appear that
pdgInvGamma (x,n,s) = pdfGamma (1/x, n, s) / x^2
or, in logs,
lpdfinvg (x,n,s) = -gammaln(n) - n*log(s) + (n+1)*log(1/x) - (1/x)/s = [ -gammaln(n) - n*log(s) + (n-1)*log(1/x) - (1/x)/s ] - 2*log(x) = [ lpdfgam(1/x,n,s) ] - 2*log(x)
Best regards
George
On 01/06/2010 16:24, George Perendia wrote:
Dear Michel
thanks
I have 2 issues:
- According to few sources,
If x~InvGamma(a,b) , then 1/x ~Gamma(a,1/b) distribution however, in our dynare inv gamma 2 (lpdfig2.m) we have different definition: X ~ IG2(s,nu) if X = inv(Z) where Z ~ G(nu/2,2/s) (Gamma distribution) (which is equivalent to If x~InvGamma(a,b) , then 1/x ~Gamma(b/2,2/a) distribution)
- I suspect there are a bugs in implementation for both invg2 and invg1
even as per our specifications of invg2 aind invg1, i.e. for invg2: we have
ldens(idx) = -gammaln(.5*nu) - (.5*nu)*(log(2)-log(s)) -.5*(nu+2)*log(x(idx)) -.5*s./x(idx);
and for e.g. s=2, n=5 and x= 0.5 as a result we have (incorrectly?):
lpdfig2(x,s,n) = -gammaln(.5*n) - (.5*n)*(log(2)-log(s)) - .5*(n+2)*log(x) -.5*s./x
ans =
0.1413whilst
lpdfgam(1/x,n/2,2/s)
ans =
-1.2450
and modified (corrected?) invgam2 with ...(nu - 2).which meets our spec.. instead of (nu + 2)...
K>> -gammaln(.5*n) - (.5*n)*(log(2)-log(s)) - .5*(n-2)*log(x) -.5*s./x
ans =
-1.2450
thus, gives same result as lpdfgam(1/x,n/2,2/s).
Similarly, invg1
ldens(idx) = log(2) - gammaln(.5*nu) - .5*nu*(log(2)-log(s)) -(nu+1)*log(x(idx)) - .5*s./(x(idx).*x(idx)) ;
lpdfig1(x,s,n)= 0.5673 = lpdfig2(x*x,s,n)= 0.5673
whilst lpdfgam(1/(x*x),n/2,2/s)
ans =
-2.2052
and modified invg1 without leading log(2) and with nu-2 instead nu+1
- gammaln(.5*nu) - .5*nu*(log(2)-log(s)) - (nu-2)*log(x) - .5*s./(x.*x)
ans =
-2.2052
i.e. = lpdfgam(1/(x*x),n/2,2/s) as per definition
(... though I may have got something wrong)
If they are bugs, shall I log them and correct them?
Best george
On 19/05/2010 20:29, Michel Juillard wrote:
Dear George,
thank you for the message. I agree that we should go with BOOST.
- In the estimation DLL, we only need a uniform and a normal random
generator and the pdfs for the normal, beta, gamma, inverse gamma 1, inverse gamma 2 and the uniform. You shouldn't try to implement more than that at this stage.
- for the inverse gamma 1 and 2 pdfs, you can directly translate the
current matlab function in ./matlab/distributions. We should add the source of the algorithm if we can remember.
- You need to look at the Matlab implementation, because we need the same
kind of generalizations:
- gamma is defined as greater than p3 rather than greater than 0
- beta is defined over [p3; p4] instead of over [0,1]
The relevant Matlab functions are set_prior.m and priordens.m
All the best,
Michel
On 05/19/2010 06:42 PM, George Perendia wrote:
Hi Sebastien
I started looking into priordens and pdfs as well as the relvant random number generators (rng-s) for estimated parameters (*) and investigated two options: C gsl and the C++ header BOOST as initially suggested by Stephane.
Both provide for the most of the requirements but there are seeral insufficiencies.
gsl provides the necessary pair of pdf and random number generators functions based on given uniform random number but not the uniform random number generator - rng (or the pdf but that is trivial to add).
On the other hand whilst BOOST provides a large variety of functions including the pdfs and cdfs in a template form as well as a variety of uniform and other distributions rngs and ways to combine them. Though it does not seem to provide a rng for beta, the beta rng may be easily derived using gamma rng and, e.g. the relation:
"If *X* and *Y* are independently distributed Gammahttp://wiki/Gamma_distribution (*α*, *θ*) and Gamma(*β*, *θ*) respectively, then *X* / (*X* + *Y*) is distributed Beta(*α*, *β*), "
...which is how gsl beta rng is implemented anyway:
Neither gsl nor Boost however seem to provide pdfs or rngs for inverse_gamma and that would need to be derived for either of them too.
Whilst the gsl is simpler it is C whlst the BOOST is C++. Also, gsl would require a library compilation before it can be used and it would also be an additional package requirement for building Dynare (i.e. whilst we are already using BOOST header library anyway).
It seems that BOOST pdfs and rngs are better documented, provide a larger variety of options for later extensions and fit better into our design strategy whilst, as part of already used BOOST, they do not require any extension on the overall maintenance and build requirements of Dynare.
I wold be inclined therefore to suggest us to use BOOST but please do let me know if you have any objections or other suggestions.
If the performance is the issue, we can run a comparative test e.g. for throwing 100000 beta rng-s and deriving their pdf-s.using gsl and BOOSt and see the difference but that may a be bit excessive - just use of a different uniform rng for generating specific distribution rng from a variety of those available in BOOST can change performance dramatically.
*) assuming the following distributions are required beta_pdf | gamma_pdf | normal_pdf | uniform_pdf | inv_gamma_pdf | inv_gamma1_pdf | inv_gamma2_pdf
Best regards George
On 18/05/2010 17:26, Sébastien Villemot wrote:
Hi George,
On Wed, May 12, 2010 at 09:36:26AM +0100, George Perendia wrote:
It was good to work together last week. Aftert yesterday's ACES meeting it looks like I may not be engaged full time on ACES for now, so, I could contiue on Estimation DLL, 1-2 days a week for time being, e.g. on prior dens, optmiiser integration, MCMC or something else. Let me know what wold be best for you.
Sorry for my late answer.
I am also happy with the progress that we made when you came to Paris.
Unfortunately, since I have been working on other projects, I have not been able to finalize the likelihood DLL yet, but I am quite close to it.
I am leaving for vacations until May 30, and I hope to finish the DLL by the first days of June. Note that I will still read and answer to emails during my vacations.
Please feel free to advance on the implementation of other parts of the design. In particular, we need to decide how to implement the prior densities. We should rely on existing code rather than reimplementing all the probability density functions. My feeling is that using the GNU Scientific Library for that purpose is the easiest way, since all the PDFs are there and it is licensed under the GPL. Do you have any other option in mind ?
Best,
--
--
--
--
Best regards George 07951415480
-- Stéphane Adjemian CEPREMAP & Université du Maine
Tel: (33)1-43-13-62-39
--
Best regards George 07951415480
Hi Stephane
Thanks.
After further inspection of our code I believe that you should not change anything in your document and that our code works OK, according to expectation and both your and Michel's documents.
The main (minor) discrepancy is just the permutation of the input parameters for inv gamma 1 and 2 in our code compensated by the permutation of the outputs in inverse_gamma_specification.m (all inline with the prior.pdf document) .
Hence, just to avoid any possibly confusion, it may be better to change our code instead and permute the nu and s parameters in calling lpdfig1 and lpdfig2 (and the internal explanations) and in the output of the inverse_gamma_specification.m modules so that we can call the log pdfs as e.g.
lpdfig1(x,n,s) = lpdfig2(x*x,n,s) = lpdfgam(1/(x*x),n/2,2/s) - 2*log(x*x)
This change will (and should ! ) not change any of our numerical results, and will be aimed just to bring the code up to more commonly used notation and the consistent derivation in your documentation which, as Sebastien already suggested, should be included in the Dynare documentation as is
In meantime, however, poss. just an additional note in the documentation pointing to the permutation in our implementation will do the job.
best regards George
On 03/06/2010 12:32, Stéphane Adjemian wrote:
Hi George,
I am not sure to understand your problem... Comparing Michel notes and the codes with my notes I can't see any difference... Except the order of declaration of the parameters nu and s in IG_2 and IG_1. I In my pdf I write IG_2(nu,s) but in the codes we wrote IG_2(s,nu) (as in Michel's notes). I will change these notations in my pdf.
Best, Stéphane.
On 3 June 2010 11:49, George Perendia <george@perendia.orangehome.co.uk mailto:george@perendia.orangehome.co.uk> wrote:
Hi Stephane Thanks for your very detailed notes. I am afraid that although they are consistent within themselves and your comment that the inv gamma2 is just reparametrised inv gamma, they however appear to be inconsistent with both Michel's priors.pdf and the current implementation where the inv gamma 2 (and 1) modified parameters appear to be permuted when applied on gamma in comparison to your notes. I.e.: % X ~ IG2(s,nu) if X = inv(Z) where Z ~ G(nu/2,2/s) (Gamma distribution) % X ~ IG1(s,nu) if X = sqrt(Y) where Y ~ IG2(s,nu) and Y = inv(Z) with Z ~ G(nu/2,2/s) (Gamma distribution) and, as result, e.g. lpdfig1(x,s,n) = lpdfig2(x*x,s,n) = lpdfgam(1/(x*x),n/2,2/s) - 2*log(x*x) = 0.5673 for my test (see below) whilst , permuting s and n as in your equation usually gives a different result, e.g. lpdfgam(1/(x*x),s/2,2/n) - 2*log(x*x) = -6.3111 It thus appear to me that the Michel's priors.pdf and our current implementation either reflect a different notion of inv gamma 1 and 2, or, our inv gamma 1 and 2 gamma-shape and gamma-scale parameters are somehow (or somewhere) permuted ? best george On 02/06/2010 15:57, Stéphane Adjemian wrote:Hi George, The inverse gamma 2 density appears naturally in Bayesian statistics when considering conjugate priors in linear Gaussian models. The inverse gamma densities are quite common in Bayesian statistics (the references you are giving are not pertinent with respect to this literature). Furthermore the inverse gamma 2 is just a re-parametrization of the probability textbook (or common if you want) inverse gamma density. This is shown in the attached pdf. So I don't think we should provide the inverse gamma prior you are suggesting in Dynare. Best, Stéphane. On 2 June 2010 10:35, George Perendia <george@perendia.orangehome.co.uk <mailto:george@perendia.orangehome.co.uk>> wrote: Hi Michel Thanks for the pdf document . 2) My apologies - Having had a chance to see the full representation of both inv gamma, 1 and 2 pdf in the document you sent I do not anymore think there is a bug in the implementation - the implementation seem to conform closely to the given representations. 1) I can not tell if the notes are incorrect but they do seem mutually consistent, so that e.g. lpdfig1(x,s,n) = lpdfig2(x*x,s,n) = lpdfgam(1/(x*x),n/2,2/s) - 2*log(x*x) = 0.5673 for my test On the other hand, both inv gamma 1 and 2 representations differ from what I find the most common definition of the inverse gamma from various sources invGammaPDF(x,a,b) = b^a / gamma(a) * 1/x ^(a+1) * exp( -b/x) which, in turn = GammaPDF(1/x,a,1/b) / x^2 (though, alternative representations of gammaPDF with inverse scale beta = 1/b would give = GammaPDF(1/x,a,beta) / x^2 ) (e.g. see the enclosed document or internet reference such as http://www.chemistrydaily.com/chemistry/Inverse-gamma_distribution or even wikipedia http://en.wikipedia.org/wiki/Inverse-gamma_distribution ) We may need to add another inv Gamma for that common representation. Best George On 01/06/2010 20:37, Michel Juillard wrote:Hi George, inverse gamma distributions are tricky because different authors use different notations. I attach the notes that I had made in 2004. Could you please tell me 1) if you think that these notes are erroneous. In this case, please, indicate the sources that you use. 2) if the code currently implemented in Dynare Matlab conforms with these notes. I'm leaving tomorrow for the Dynare Conference. I will be back this week end and will try to look into the issue next week. All the best, Michel On 06/01/2010 09:00 PM, George Perendia wrote:PS: On the other hand, using the commonly used inv gamma formula pdgInvGamma (x,n,s) = 1/(gamma(n) * s^n) * x ^(-n-1) exp (1/ (x*s) ) it appear that pdgInvGamma (x,n,s) = pdfGamma (1/x, n, s) / x^2 or, in logs, lpdfinvg (x,n,s) = -gammaln(n) - n*log(s) + (n+1)*log(1/x) - (1/x)/s = [ -gammaln(n) - n*log(s) + (n-1)*log(1/x) - (1/x)/s ] - 2*log(x) = [ lpdfgam(1/x,n,s) ] - 2*log(x) Best regards George On 01/06/2010 16:24, George Perendia wrote:Dear Michel thanks I have 2 issues: 1) According to few sources, If x~InvGamma(a,b) , then 1/x ~Gamma(a,1/b) distribution however, in our dynare inv gamma 2 (lpdfig2.m) we have different definition: X ~ IG2(s,nu) if X = inv(Z) where Z ~ G(nu/2,2/s) (Gamma distribution) (which is equivalent to If x~InvGamma(a,b) , then 1/x ~Gamma(b/2,2/a) distribution) 2) I suspect there are a bugs in implementation for both invg2 and invg1 even as per our specifications of invg2 aind invg1, i.e. for invg2: we have ldens(idx) = -gammaln(.5*nu) - (.5*nu)*(log(2)-log(s)) - .5*(nu+2)*log(x(idx)) -.5*s./x(idx); and for e.g. s=2, n=5 and x= 0.5 as a result we have (incorrectly?): lpdfig2(x,s,n) = -gammaln(.5*n) - (.5*n)*(log(2)-log(s)) - .5*(n+2)*log(x) -.5*s./x ans = 0.1413 whilst lpdfgam(1/x,n/2,2/s) ans = -1.2450 and modified (corrected?) invgam2 with ...(nu - 2).which meets our spec.. instead of (nu + 2)... K>> -gammaln(.5*n) - (.5*n)*(log(2)-log(s)) - .5*(n-2)*log(x) -.5*s./x ans = -1.2450 thus, gives same result as lpdfgam(1/x,n/2,2/s). Similarly, invg1 ldens(idx) = log(2) - gammaln(.5*nu) - .5*nu*(log(2)-log(s)) - (nu+1)*log(x(idx)) - .5*s./(x(idx).*x(idx)) ; lpdfig1(x,s,n)= 0.5673 = lpdfig2(x*x,s,n)= 0.5673 whilst lpdfgam(1/(x*x),n/2,2/s) ans = -2.2052 and modified invg1 without leading log(2) and with nu-2 instead nu+1 - gammaln(.5*nu) - .5*nu*(log(2)-log(s)) - (nu-2)*log(x) - .5*s./(x.*x) ans = -2.2052 i.e. = lpdfgam(1/(x*x),n/2,2/s) as per definition (... though I may have got something wrong) If they are bugs, shall I log them and correct them? Best george On 19/05/2010 20:29, Michel Juillard wrote:Dear George, thank you for the message. I agree that we should go with BOOST. 1) In the estimation DLL, we only need a uniform and a normal random generator and the pdfs for the normal, beta, gamma, inverse gamma 1, inverse gamma 2 and the uniform. You shouldn't try to implement more than that at this stage. 2) for the inverse gamma 1 and 2 pdfs, you can directly translate the current matlab function in ./matlab/distributions. We should add the source of the algorithm if we can remember. 3) You need to look at the Matlab implementation, because we need the same kind of generalizations: - gamma is defined as greater than p3 rather than greater than 0 - beta is defined over [p3; p4] instead of over [0,1] The relevant Matlab functions are set_prior.m and priordens.m All the best, Michel On 05/19/2010 06:42 PM, George Perendia wrote:> Hi Sebastien > > I started looking into priordens and pdfs as well as the > relvant random number generators (rng-s) for estimated > parameters (*) and investigated two options: C gsl and > the C++ header BOOST as initially suggested by Stephane. > > Both provide for the most of the requirements but there > are seeral insufficiencies. > > gsl provides the necessary pair of pdf and random number > generators functions based on given uniform random > number but not the uniform random number generator - rng > (or the pdf but that is trivial to add). > > On the other hand whilst BOOST provides a large variety > of functions including the pdfs and cdfs in a template > form as well as a variety of uniform and other > distributions rngs and ways to combine them. Though it > does not seem to provide a rng for beta, the beta rng > may be easily derived using gamma rng and, e.g. the > relation: > > "If /X/ and /Y/ are independently distributed Gamma > http://wiki/Gamma_distribution(/α/, /θ/) and > Gamma(/β/, /θ/) respectively, then /X/ / (/X/ + /Y/) is > distributed Beta(/α/, /β/), " > > ...which is how gsl beta rng is implemented anyway: > > Neither gsl nor Boost however seem to provide pdfs or > rngs for inverse_gamma and that would need to be derived > for either of them too. > > Whilst the gsl is simpler it is C whlst the BOOST is > C++. Also, gsl would require a library compilation > before it can be used and it would also be an additional > package requirement for building Dynare (i.e. whilst we > are already using BOOST header library anyway). > > It seems that BOOST pdfs and rngs are better documented, > provide a larger variety of options for later extensions > and fit better into our design strategy whilst, as part > of already used BOOST, they do not require any extension > on the overall maintenance and build requirements of Dynare. > > I wold be inclined therefore to suggest us to use BOOST > but please do let me know if you have any objections or > other suggestions. > > If the performance is the issue, we can run a > comparative test e.g. for throwing 100000 beta rng-s and > deriving their pdf-s.using gsl and BOOSt and see the > difference but that may a be bit excessive - just use of > a different uniform rng for generating specific > distribution rng from a variety of those available in > BOOST can change performance dramatically. > > *) assuming the following distributions are required > |beta_pdf| | |gamma_pdf| | |normal_pdf| | |uniform_pdf| > | |inv_gamma_pdf| | |inv_gamma1_pdf| | |inv_gamma2_pdf| > > Best regards > George > > On 18/05/2010 17:26, Sébastien Villemot wrote: >> Hi George, >> >> On Wed, May 12, 2010 at 09:36:26AM +0100, George Perendia wrote: >> >>> It was good to work together last week. >>> Aftert yesterday's ACES meeting it looks like I may not be engaged >>> full time on ACES for now, so, I could contiue on Estimation DLL, >>> 1-2 days a week for time being, e.g. on prior dens, optmiiser >>> integration, MCMC or something else. >>> Let me know what wold be best for you. >>> >> Sorry for my late answer. >> >> I am also happy with the progress that we made when you came to Paris. >> >> Unfortunately, since I have been working on other projects, I have not been able to finalize the likelihood DLL yet, but I am quite close to it. >> >> I am leaving for vacations until May 30, and I hope to finish the DLL by the first days of June. Note that I will still read and answer to emails during my vacations. >> >> Please feel free to advance on the implementation of other parts of the design. In particular, we need to decide how to implement the prior densities. We should rely on existing code rather than reimplementing all the probability density functions. My feeling is that using the GNU Scientific Library for that purpose is the easiest way, since all the PDFs are there and it is licensed under the GPL. Do you have any other option in mind ? >> >> Best, >> >> > > > -- > > > > >
------ Best regards George -- Stéphane Adjemian CEPREMAP & Université du Maine Tel: (33)1-43-13-62-39-- Best regards George-- Stéphane Adjemian CEPREMAP & Université du Maine
Tel: (33)1-43-13-62-39
Dev mailing list Dev@dynare.org https://www.dynare.org/cgi-bin/mailman/listinfo/dev
-
George Perendia -
Stéphane Adjemian