wolffd@0: %DEMGAUSS Demonstrate sampling from Gaussian distributions.
wolffd@0: %
wolffd@0: %	Description
wolffd@0: %
wolffd@0: %	DEMGAUSS provides a simple illustration of the generation of data
wolffd@0: %	from Gaussian distributions. It first samples from a one-dimensional
wolffd@0: %	distribution using RANDN, and then plots a normalized histogram
wolffd@0: %	estimate of the distribution using HISTP together with the true
wolffd@0: %	density calculated using GAUSS.
wolffd@0: %
wolffd@0: %	DEMGAUSS then demonstrates sampling from a Gaussian distribution in
wolffd@0: %	two dimensions. It creates a mean vector and a covariance matrix, and
wolffd@0: %	then plots contours of constant density using the function GAUSS. A
wolffd@0: %	sample of points drawn from this distribution, obtained using the
wolffd@0: %	function GSAMP, is then superimposed on the contours.
wolffd@0: %
wolffd@0: %	See also
wolffd@0: %	GAUSS, GSAMP, HISTP
wolffd@0: %
wolffd@0: 
wolffd@0: %	Copyright (c) Ian T Nabney (1996-2001)
wolffd@0: 
wolffd@0: clc
wolffd@0: mean = 2; var = 5; nsamp = 3000;
wolffd@0: xmin = -10; xmax = 10; nbins = 30;
wolffd@0: disp('Demonstration of sampling from a uni-variate Gaussian with mean')
wolffd@0: dstring = [num2str(mean), ' and variance ', num2str(var), '.  ', ...
wolffd@0:     num2str(nsamp), ' samples are taken.'];
wolffd@0: disp(dstring);
wolffd@0: x = mean + sqrt(var)*randn(nsamp, 1);
wolffd@0: fh1 = figure;
wolffd@0: histp(x, xmin, xmax, nbins);
wolffd@0: hold on;
wolffd@0: axis([xmin xmax 0 0.2]);
wolffd@0: plotvals = linspace(xmin, xmax, 200)';
wolffd@0: probs = gauss(mean, var, plotvals);
wolffd@0: plot(plotvals, probs, '-r');
wolffd@0: xlabel('X')
wolffd@0: ylabel('Density')
wolffd@0: 
wolffd@0: disp(' ')
wolffd@0: disp('Press any key to continue')
wolffd@0: pause; 
wolffd@0: mu = [3 2];
wolffd@0: lam1 = 0.5;
wolffd@0: lam2 = 5.0;
wolffd@0: Sigma = lam1*[1,1]'*[1,1] + lam2*[1,-1]'*[1,-1];
wolffd@0: disp(' ')
wolffd@0: disp('Demonstration of sampling from a bi-variate Gaussian.  The mean is')
wolffd@0: dstring = ['[', num2str(mu(1)), ', ', num2str(mu(2)), ...
wolffd@0:       '] and the covariance matrix is'];
wolffd@0: disp(dstring)
wolffd@0: disp(Sigma);
wolffd@0: ngrid = 40;
wolffd@0: cmin = -5; cmax = 10; 
wolffd@0: cvals = linspace(cmin, cmax, ngrid);
wolffd@0: [X1, X2] = meshgrid(cvals, cvals);
wolffd@0: XX = [X1(:), X2(:)];
wolffd@0: probs = gauss(mu, Sigma, XX);
wolffd@0: probs = reshape(probs, ngrid, ngrid);
wolffd@0: 
wolffd@0: fh2 = figure;
wolffd@0: contour(X1, X2, probs, 'b');
wolffd@0: hold on
wolffd@0: 
wolffd@0: nsamp = 300;
wolffd@0: dstring = [num2str(nsamp), ' samples are generated.'];
wolffd@0: disp('The plot shows the sampled data points with a contour plot of their density.')
wolffd@0: samples = gsamp(mu, Sigma, nsamp);
wolffd@0: plot(samples(:,1), samples(:,2), 'or');
wolffd@0: xlabel('X1')
wolffd@0: ylabel('X2')
wolffd@0: grid off;
wolffd@0: 
wolffd@0: disp(' ')
wolffd@0: disp('Press any key to end')
wolffd@0: pause; 
wolffd@0: close(fh1);
wolffd@0: close(fh2);
wolffd@0: clear all;