error_ellipse

二维空间协方差矩阵可视化为一个误差椭圆的matlab和C++代码，C++代码应用到了opencv的库函数，所以如果需要运行这个代码需要配置opencv环境。

clear all;
close all;

% Create some random data
s = [2 2];
x = randn（3341）;
y1 = normrnd（s（1）.*x1）;
y2 = normrnd（s（2）.*x1）;
data = [y1 y2];

% Calculate the eigenvectors and eigenvalues
covariance = cov（data）;
[eigenvec eigenval ] = eig（covariance）;

% Get the index of the largest eigenvector
[largest_eigenvec_ind_c r] = find（eigenval == max（max（eigenval）））;
largest_eigenvec = eigenvec（: largest_eigenvec_ind_c）;

% Get the largest eigenvalue
largest_eigenval = max（max（eigenval））;

% Get the smallest eigenvector and eigenvalue
if（largest_eigenvec_ind_c == 1）
    smallest_eigenval = max（eigenval（:2））
    smallest_eigenvec = eigenvec（:2）;
else
    smallest_eigenval = max（eigenval（:1））
    smallest_eigenvec = eigenvec（1:）;
end

% Calculate the angle between the x-axis and the largest eigenvector
angle = atan2（largest_eigenvec（2） largest_eigenvec（1））;

% This angle is between -pi and pi.
% Let‘s shift it such that the angle is between 0 and 2pi
if（angle < 0）
    angle = angle + 2*pi;
end

% Get the coordinates of the data mean
avg = mean（data）;

% Get the 95% confidence interval error ellipse
chisquare_val = 2.4477;
theta_grid = linspace（02*pi）;
phi = angle;
X0=avg（1）;
Y0=avg（2）;
a=chisquare_val*sqrt（largest_eigenval）;
b=chisquare_val*sqrt（smallest_eigenval）;

% the ellipse in x and y coordinates 
ellipse_x_r  = a*cos（ theta_grid ）;
ellipse_y_r  = b*sin（ theta_grid ）;

%Define a rotation matrix
R = [ cos（phi） sin（phi）; -sin（phi） cos（phi） ];

%let‘s rotate the ellipse to some angle phi
r_ellipse = [ellipse_x_r;ellipse_y_r]‘ * R;

% Draw the error ellipse
plot（r_ellipse（:1） + X0r_ellipse（:2） + Y0‘-‘）
hold on;

% Plot the original data
plot（data（:1） data（:2） ‘.‘）;
mindata = min（min（data））;
maxdata = max（max（data））;
Xlim（[mindata-3 maxdata+3]）;
Ylim（[mindata-3 maxdata+3]）;
hold on;

% Plot the eigenvectors
quiver（X0 Y0 largest_eigenvec（1）*sqrt（largest_eigenval） largest_eigenvec（2）*sqrt（largest_eigenval） ‘-m‘ ‘LineWidth‘2）;
quiver（X0 Y0 smallest_eigenvec（1）*sqrt（smallest_eigenval） smallest_eigenvec（2）*sqrt（smallest_eigenval） ‘-g‘ ‘LineWidth‘2）;
hold on;

% Set the axis labels
hXLabel = xlabel（‘x‘）;
hYLabel = ylabel（‘y‘）;

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        2265  2015-05-23 01:00  error_ellipse.m
     文件        1785  2015-05-23 01:02  error_ellipse.txt