randpdf2(x, y, p, m)

function z=randpdf2(x, y, p, m)
% RANDPDF2
% Random numbers from a user defined 2-dimensional distribution
%
% SYNTAX:
% z = randpdf2(x, y, p, n)
% randpdf2(x, y, p, n)
%
% INPUT:
% x - x random variable - first dimension (1 x k) vector
% y - y random variable - second dimension (1 x k) vector
% p - probability density (1 x k) vector
% m - length of the output vector (scalar)
%
% OUTPUT:
% z - A (2 x n) matrix of random pairs (x,y).
% Run function without output for some plots.
%
% DESCRIPTION:
% z=randpdf2(x, y, p, n) returns the matrix of random numbers from
% probability density distribution defined by (x,y) and p. p is a 1 x k
% vector which is the probability density (pdf).
% x, y, and p must be of the same length.
% p must be non-negative (as a valid pdf). Its integral does not have to
% be 1, as it is normalized by randsample().
% m defines the output matrix length. Hence, the output is a 2 x n matrix
%
%
%
% EXAMPLE-1:
% Generation of random numbers from the distribution: f(x,y)=x+y for 0<=x,y<=1
% When one computes the correlation between x and y, one gets -1/11
%
% x=0:1e-2:1;
% y=x;
% [X Y]=meshgrid(x,y);
% X=X(:); % column vectors
% Y=Y(:);
% p=X+Y; % pdf
% X=X'; % row vectors
% Y=Y';
% p=p';
% Z=randpdf2(X,Y,p,100000);
% corr(Z')
% plot a 3D histogram
% hist3(Z')
% xlabel('X')
% ylabel('Y')
%
% EXAMPLE-2:
% A discrete distribution. x=-1 or 1; y=-1, 0, or 1 with the following
% probabilities:
%
% | -1 | 1 |
% -----------------
% -1| 1/8 | 1/8 |
% -----------------
% 0| 1/8 | 1/8 |
% -----------------
% 1| 2/8 | 2/8 |
%
% x=[-1 1];
% y=[-1 0 1];
% [X Y]=meshgrid(x,y);
% X=X(:);
% Y=Y(:);
% p=[1/8 1/8 2/8 1/8 1/8 2/8];
% X=X';
% Y=Y';
% Z=randpdf2(X,Y,p,100000);
%
% plot a 3D histogram
% ctrs={[-1 1]; [-1 0 1]};
% hist3(Z',ctrs)
% xlabel('X')
% ylabel('Y')