matlab 误差椭圆,求3倍标准差误差椭圆分析的程序

根据《白话空间统计之九：方向分布(标准差椭圆)修正版》(有些地方没有理解清楚)，写了下面的程序。但是好像结果不对

Z=mvnrnd([0.5 1.5], [0.025 0.03 ; 0.03 0.16], 50);

X=Z(:,1);  Y=Z(:,2);

mean_X=nanmean(X);  mean_Y=nanmean(Y);   %椭圆圆心

%确定长短半轴

SDE_X=nanstd(X,1);

SDE_Y=nanstd(Y,1);

%确定方向

N=size(X,1)*size(X,2)-size(find(isnan(X)),1);  %非NaN元素个数

X2=nanvar(X(:),1)*N;

Y2=nanvar(Y(:),1)*N;

A1=nancov(X,Y,1);  % 是2*2矩阵，含4个元素，元素1是S(X(:))，及X的方差；元素4是Y的方差；元素2与3相等，是X与Y的协方差；

X1=A1(1)*N; % (X-mean_X)平方求和

Y1=A1(4)*N; % (Y-mean_Y) 平方求和

XY=A1(2)*N; %(X-mean_X)(Y-mean_Y) 求和

A=X1-Y1;

B=sqrt(A^2+4*XY^2);

C=2*XY;

theta=atan((A+B)/C);  % 椭圆方向，狐度，以X轴为准，正北方(12点方向)为0度，顺时针旋转    *180/pi=角度

%确定x、y轴的标准差

XStd=sqrt((cos(theta)^2*X1+sin(theta)^2*Y1-sin(2*theta)*XY)/N)*sqrt(2);  %X轴标准差，不知道是否有*sqrt(2)

YStd=sqrt((sin(theta)^2*X1+cos(theta)^2*Y1+sin(2*theta)*XY)/N)*sqrt(2);  %Y轴标准差，不知道是否有*sqrt(2)

% XStd=sqrt((cos(theta)^2*X1+sin(theta)^2*Y1-sin(2*theta)*XY)/N);  %X轴标准差，不知道是否有*sqrt(2)

% YStd