Generate Random Numbers from a 2D Discrete Distribution

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