CODEWITHSUNDEEP
matlabrandomdistributiongamma-distribution
Kelsey
Kelsey

Reputation: 21

How to draw random numbers from a gamma distribution without the Statistics Toolbox?

I am varying the signal strength for synthetic images. I need the signal to vary between 0 and 0.1, but I need to do this with a gamma distribution so that more of them fall around the .01/.02 range. The problem is that I am using the 2010 version of Matlab without the Statistics Toolbox that doesn't have the gamrnd function a part of its library.

Any and all help is greatly appreciated.

Upvotes: 2

Views: 2261

Answers (1)

EBH
EBH

Reputation: 10440

You can use the Inverse transform sampling method to convert a uniform distribution to any other distribution:

P = rand(1000);
X = gaminv(P(:),2,2); % with k = 2 and theta = 2

Here is a litle demonstration:

for k = [1 3 9]
    for theta = [0.5 1 2]
        X = gaminv(P(:),k,theta);
        histogram(X,50)
        hold on
    end
end

Which gives:

gamma dist

Edit:

Without the statistics toolbox, you can use the Marsaglia's simple transformation-rejection method to generate random numbers from gamma distribution with rand and randn:

N = 10000; % no. of tries
% distribution parameters:
a = 0.5;
b = 0.1;
% Marsaglia's simple transformation-rejection:
d = a - 1/3;
x = randn(N,1);
U = rand(N,1);
v = (1+x./sqrt(9*d)).^3;
accept = log(U)<(0.5*x.^2+d-d*v+d*log(v));
Y = d*(v(accept)).*b;

Now Y is distributed like gamma(a,b). We can test the result using the gamrnd function:

n = size(Y,1);
X = gamrnd(a,b,n,1);

And the histograms of Y, and X are:

gamma dist 2

However, keep in mind that gamma distribution might not fit your needs because it has no specific upper bound (i.e. goes to infinity). So you may want to use another (bounded) distribution, like beta divided by 10.

Upvotes: 5

Related Questions