plotting contours according to distribution variance

Question

I'm running an experiment and I get as output a set of data in the form of 2D points

# read csv file
samples = csvread('results.csv');

assuming they are normally distributed I can compute their mean and covariance and evaluate the pdf

mu = mean(samples);
sigma = cov(samples);
y = mvnpdf(X,mu,sigma);

where X is a 2D grid computed with meshgrid.

What I would like to do now is to plot this distribution with iso-contours as here.

And, I would like to choose contours at specific values that are related to the covariance as in this picture:

As an example, I should know the value of the distribution in the point = (X: mu(1)+sigma(1,1), Y: mu(1)+sigma(1,1)). My guess is that I could use the normpdf function in this way

value = normpdf(point,mu,sigma)

but I get this as output:

value =

   116.39   297.63
   297.63   409.88

and I don't really know how to interpret it.

Can please someone show me how to do it?

Thanks.

Answer 1

I've lost you at the point where you started talking about a normpdf. You have a multivariate normal and that's what you need to work on. Not sure what you expect to get from the univariate normal here.

If you want to get the probability value of the multivariate normal at a mahalanobis distance of 1 (i.e. one standard deviation away from the mean of the distribution), you can do that by evaluating the mvnpdf at the value [m+s,0] where, m = mu(1) and s = sqrt( sigma(1,1) ).

Then you can use that with an appropriate contour function to get your desired contour.

EDIT: Here's an example.

  pkg load statistics

% create 1000 samples from a known multivariate normal
  Observations     = mvnrnd( [0,0], [4, 1; 1 ,2], 1000 );

% Get mean and covariance estimates from sample
  SampleMean       = mean( Observations, 1 ) % average along rows
  SampleCovariance = cov( Observations )

% Get mean and standard deviation in first dimension (i.e. "x-axis")
  Mean_firstDimension   = SampleMean(1)
  StdDev_firstDimension = sqrt( SampleCovariance(1,1) )

% Get gaussian values at malanobis distance of 1, 2, and 3
  MVN_value_at_mahalanobis_distance_of_one   = mvnpdf( [ Mean_firstDimension +     StdDev_firstDimension, 0], SampleMean, SampleCovariance )
  MVN_value_at_mahalanobis_distance_of_two   = mvnpdf( [ Mean_firstDimension + 2 * StdDev_firstDimension, 0], SampleMean, SampleCovariance )    
  MVN_value_at_mahalanobis_distance_of_three = mvnpdf( [ Mean_firstDimension + 3 * StdDev_firstDimension, 0], SampleMean, SampleCovariance )    

% Define grid:
  [GridX, GridY] = ndgrid( -8:0.1:8, -8:0.1:8 );
  GridValues = mvnpdf( [GridX(:), GridY(:)], SampleMean, SampleCovariance );
  GridValues = reshape( GridValues, size( GridX ) );

% Plot Observations
  plot( Observations(:,1), Observations(:,2), 'o', 'markerfacecolor', 'g', 'markeredgecolor', [0,0.5,0], 'linewidth', 1.5 );
  hold on;

% Plot contours over grid
  contour( GridX, GridY, GridValues,   ...
           [ MVN_value_at_mahalanobis_distance_of_three, ...
             MVN_value_at_mahalanobis_distance_of_two, ...
             MVN_value_at_mahalanobis_distance_of_one ...
           ],
           'linewidth', 2
  )
  hold off;

% Set nice limits and colours for background
  axis([-8, +8, -8, +8]); axis equal square;
  set(gca, 'color', 'k');
  set(gcf, 'color', [0.75, 0.75, 0.75]);