1

I am trying to plot the ksdensity of an Inverse gamma(0.001, 0.001) but the plot has only one point. The commands I used are

 alpha1 = 0.001;
    beta1 = 0.001;
    n = 1e+5;
    r=1./gamrnd(alpha1,1/beta1,n,1);  
   [f,xi] = ksdensity(r);
    plot(xi,f,'--m');

The first term of f is a real number and all the others are NaN. The first term of xi is a real number and all the others are Inf.

Could you please me help me with this.

Thank you very much.

David van Driessche
  • 6,602
  • 2
  • 28
  • 41
F.F.
  • 89
  • 7

1 Answers1

1

Did you look at the generated data in r? I get almost 49 % of the time the value Inf. It appears that many values of the Gamma distribution with the parameters you chose are 0, or actually, so small that they cannot be represented in the standard double precision numeric format of Matlab (the smallest possible value is 2.225 · 10^-308).

If you look at the Wikipedia page for the inverse Gamma distribution, you will see that the mean is not defined for alpha <= 1 and the variance is not defined for alpha <= 2, etc. According to Wolfram Alpha, the mode of your distribution is at 0.000999001, while the median is at 1.90687 · 10^298 (only a few orders of magnitude below the largest possible double value, 1.797 · 10^308!).

The density around the mode

x = 0.0001:0.0001:0.1;
plot(x, beta1 ^ alpha1 / gamma(alpha1) .* x .^ -(alpha1 + 1) .* exp(- beta1 ./ x))

looks like this

but this covers only a tiny fraction of the total distribution (the 1st percentile is at 41.2211).

So the problem here is not the kernel density estimation, but that the distribution you are looking at has extreme properties which make it hard to plot the density from the analytical formula, let alone estimate it from simulated random numbers.

A. Donda
  • 8,381
  • 2
  • 20
  • 49