% This function draws N-dimensional points from a probability
% distribution that is a sum of a narrow gaussian (sigma=1) and
% a wide gaussian. The latter term of the PDF generates occasional
% outliers.
function [r, Nflyers, flyers] = pdf_flyers(Npoints, N, Pflyer, sigma_flyer)
% Npoints = the number of N-dimensional points 
% N = number of dimensions
% Pflyer = probability that a given point will be an outlier (a "flyer").
% sigma_flyer = standard deviation of the Gaussian random term
%     added to flyers.

r = randn(Npoints, N)/sqrt(N);
flyers = rand(Npoints,1) < Pflyer; % logical indicating which points are flyers.
Nflyers = sum(flyers); % number of flyers.
r(flyers,:) = r(flyers,:) + sigma_flyer * randn(Nflyers, N)/sqrt(N);