Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Please help - computing inverse gamma PDFs

Subject: Please help - computing inverse gamma PDFs

From: Misha Koshelev

Date: 1 Apr, 2010 21:19:04

Message: 1 of 3

Dear All:

I am drawing numbers from an inverse gamma with parameters a and b using three different methods. All methods give same results per histogram.

However, when computing the PDFs I am unable to use any MATLAB functions to do this. Is there a MATLAB function I can use? I have tried quite a few permutations below but do not seem to get correct result.

Thank you!
Misha

%
% histograms are the same
%

a=2;
b=2;
N=10000;

v=(2*b)./chi2rnd(2*a,1,N);
fig=figure('visible','off');
hist(log(v),100);
print(fig,'-dpdf','1.pdf');

v=1./gamrnd(a,1/b,1,N);
fig=figure('visible','off');
hist(log(v),100);
print(fig,'-dpdf','2.pdf');

v=b./randgamma(a*ones(1,N));
fig=figure('visible','off');
hist(log(v),100);
print(fig,'-dpdf','3.pdf');

%
% why aren't the PDFs???
%

%
% START WITH GAMMA(a,b)
%
s=a*b;

% these are same - matlab parameterization
1/(b^a*gamma(a))*s^(a-1)*exp(-s/b)
gampdf(s,a,b)

% logs - matlab parameterization
-(a*log(b)+gammaln(a))+(a-1)*log(s)-s/b
log(gampdf(s,a,b))

% alternate parameterization
s=a/b;
a*log(b)-gammaln(a)+(a-1)*log(s)-s*b
log(gampdf(s,a,1/b))

% now inverse gamma - no longer works - why???
s=b/(a-1);
a*log(b)-gammaln(a)-(a+1)*log(s)-b/s
log(gampdf(s^-1,a,1/b))

% others that don't work
chi2pdf((2*b)./s,2*a)

Subject: Please help - computing inverse gamma PDFs

From: Nick

Date: 2 Apr, 2010 00:34:14

Message: 2 of 3

x = 0.01:0.01:7;
invgam = @(x,a,b) b^a/gamma(a).*(1./x).^(a+1).*exp(-b./x);
plot(x, invgam(x,2,2));

An inverse gamma pdf.

Subject: Please help - computing inverse gamma PDFs

From: Misha Koshelev

Date: 2 Apr, 2010 01:19:04

Message: 3 of 3

Nick <nickchng@gmail.com> wrote in message <902172621.487835.1270168484560.JavaMail.root@gallium.mathforum.org>...
> x = 0.01:0.01:7;
> invgam = @(x,a,b) b^a/gamma(a).*(1./x).^(a+1).*exp(-b./x);
> plot(x, invgam(x,2,2));
>
> An inverse gamma pdf.

Ah. Thank you!

Misha

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us