File Exchange

## Generate Random Numbers from a 2D Discrete Distribution

version 1.3.0.0 (6.9 KB) by Tristan Ursell

### Tristan Ursell (view profile)

Random numbers from any 2D discrete probability distribution, at any resolution.

Updated 10 Feb 2016

Tristan Ursell
2D Random Number Generator for a Given Discrete Distribution
March 2012
[x0,y0]=pinky(Xin,Yin,dist_in,varargin);
'Xin' is a vector specifying the equally spaced values along the x-axis.
'Yin' is a vector specifying the equally spaced values along the y-axis.
'dist_in' (dist_in > 0) is a matrix with dimensions length(Yin) x length(Xin), whose values specify a 2D discrete probability distribution. The distribution does not need to be normalized.

'res' (res > 1) is a multiplicative increase in the resolution of chosen random number, using cubic spline interpolation of the values in 'dist_in'. Using the 'res' option can significantly slow down the code, due to the computational costs of interpolation, but allows one to generate more continuous values from the distribution.

[x0,y0] is the output doublet of random numbers consistent with dist_in.

The 'gendist' function required by this script, is included in this m-file.

Example with an anisotropic Gaussian at native resolution:

Xin=-10:0.1:10;
Yin=-5:0.1:5;

Xmat = ones(length(Yin),1)*Xin;
Ymat = Yin'*ones(1,length(Xin));

D1=1;
D2=3;
Dist=exp(-Ymat.^2/(2*D1^2)-Xmat.^2/(2*D2^2));

N=10000;

vals=zeros(2,N);
for i=1:N
[vals(1,i),vals(2,i)]=pinky(Xin,Yin,Dist);
end

figure;
hold on
imagesc(Xin,Yin,Dist)
colormap(gray)
plot(vals(1,:),vals(2,:),'r.')
xlabel('Xin')
ylabel('Yin')
axis equal tight
box on

Example with multiple peaks at 10X resolution:

Xin=-10:0.1:10;
Yin=-5:0.1:5;

Xmat = ones(length(Yin),1)*Xin;
Ymat = Yin'*ones(1,length(Xin));

D1=0.5;
D2=1;
Dist=exp(-(Ymat-3).^2/(4*D2^2)-(Xmat+6).^2/(2*D1^2))+exp(-(Ymat+2).^2/(2*D1^2)-(Xmat+1).^2/(2*D2^2))+exp(-(Ymat-1).^2/(2*D1^2)-(Xmat-2).^2/(2*D2^2));

N=10000;
res=10;
vals=zeros(2,N);
for i=1:N
[vals(1,i),vals(2,i)]=pinky(Xin,Yin,Dist,res);
end

figure;
hold on
imagesc(Xin,Yin,Dist)
colormap(gray)
plot(vals(1,:),vals(2,:),'r.')
xlabel('Xin')
ylabel('Yin')
axis equal tight
box on

Derrick Gao

H H

### H H (view profile)

For those who experiences a flipped plot when using imagesc (it may happen if you define Dist in a different way), you can use this:
imagesc(rot90(flipud(Xin),-1), rot90(flipud(Yin),-1) ,Dist)
set(gca,'YDir','normal')

Tristan Ursell

### Tristan Ursell (view profile)

Hi HH -- thanks for the comment and I'm glad it's working for you. The code has *no* for-loops, I think what you're referring to is the example code I include with the help text ... that implements this code so that you can see it in use, the m-file has no for-loops. Hope that clarifies.

Ahinoam Pollack

### Ahinoam Pollack (view profile)

Super nice.
Thanks!

Alexander Kharlan

### Alexander Kharlan (view profile)

Good function indeed, only I have got an issue calling this one thousands of times in a FOR loop. Really slows stuff down. I am struggling to get it to generate multiple x-y couples from one distribution with a single function call - as, for example, randsample does. But I'm not having much success just yet. Does anyone have a clue on how do I generate loads of random numbers at once?

Steven White

### Steven White (view profile)

Very useful - thanks!

Tristan Ursell

### Tristan Ursell (view profile)

Good catch -- I just put a check in the code to make sure the interpolated values are always positive.

M S

### M S (view profile)

Excellent code, and amazingly, the only FEX submission for random sampling from 2D distributions.

I did run into a little bug. When using interpolation (res>1), I occasionally get "Error: All elements of first argument, P, must be positive." thrown by gendist(). This looks like it comes from interpolation giving negative numbers. The only hint at the culprit is that the joint PDFs I'm sending to your function have very very small numbers, (e.g., 1.4822e-323). When more reasonable numbers are used, the interpolation is fine.

Filip Trönnberg

### Filip Trönnberg (view profile)

Wonderful function!!