Prev: How to save an axes inside a GUI?
Next: HELP: Object Oriented Programming
From: Misha Koshelev on 1 Apr 2010 17:19
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)
From: Nick on 1 Apr 2010 16:34
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.
From: Misha Koshelev on 1 Apr 2010 21:19
Nick <nickchng(a)gmail.com> wrote in message <902172621.487835.1270168484560.JavaMail.root(a)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
 | 
Pages: 1
Prev: How to save an axes inside a GUI?
Next: HELP: Object Oriented Programming