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 |