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(α, θ) 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