function out = randn_bogotify(N, M, p)

% Produce an NxM matrix of samples from a standard normal 
%  distribution, except: with probability p, these values are
%  set to NaN

if(nargin < 1)
  N = 1;
end;

if(nargin < 2)
  M = 1;
end;

if(nargin < 3)
  p = 0.5;
end;

% Sample of size NxM
out = randn(N, M);

censored_elems = find(rand(N,M) < p);
out(censored_elems) = NaN;

