MATLAB代码
%plot confidence ellipse
%run in MATLAB2018b
t = [0.9 1.1 4.8 3.2 7.8 6.8 11.6 1.6 1.2 7.2 3.2]';
x1 = [67.3 111.3 173.0 80.8 199.7 139.4 368.2 95.7 109.6 196.2 102.2]';
% t = rand(100,1);
% x1 = rand(100,1);
A = [t,x1];
p = 0.95;
errorEllipse(A,p)
function errorEllipse(datamatrix,p)
%print 2-demension confidence ellipse
%In:n×2 matrix, confidence probability p
data = datamatrix;
covariance = cov(data);
[eigenvec,eigenval] = eig(covariance);
[sortEigenval,index] = sort(diag(eigenval),'descend');
sortEigenvec = eigenvec(:,index);
largestEigenval = sortEigenval(1);
smallestEigenval = sortEigenval(end); %find the minimum eigenvalue
largestEigenvec = sortEigenvec(:,1); %find the maximum eigenvector
angle = atan2(largestEigenvec(2), largestEigenvec(1)); %calculate the angle between x-axis and the maximum eigenvector, [-pi,pi]
if(angle < 0)
angle = angle + 2*pi;
end
avg = mean(data); %calculate the mean of two columns of data
%configure the parameters of the confidence ellipse
chisquareVal = sqrt(chi2inv(p,2)); %chi-square value
thetaGrid = linspace(0,2*pi);
phi = angle; %rotation angle
X0=avg(1);
Y0=avg(2);
a=chisquareVal*sqrt(largestEigenval); %the wheelbase length
b=chisquareVal*sqrt(smallestEigenval);
ellipseXR = a*cos( thetaGrid ); %onto rectangular axis
ellipseYR = b*sin( thetaGrid );
R = [ cos(phi) sin(phi); -sin(phi) cos(phi) ]; %rotation matrix
rEllipse = [ellipseXR;ellipseYR]' * R; %rotation
plot(rEllipse(:,1) + X0,rEllipse(:,2) + Y0,'-',data(:,1), data(:,2), '.') %print
axis square
end
运行结果:
![9f9ca5bbd03e6aeec861d871afe1c8d0.png](https://i-blog.csdnimg.cn/blog_migrate/f655bb4eda7c054ab55d402f7f99e281.png)