利用MATLAB绘制置信区域

<MATLAB小技巧之二十四：利用MATLAB绘制置信区域>

***************************************

统计中经常会遇到求置信区间、置信区域（如置信椭圆、置信椭球）等，有时候需要把置信区域画出来，这样看起看更为直观，下面结合具体案例介绍调用自编函数ConfidenceRegion绘制置信区域。

【例1】绘制置信区间。产生一元正态分布随机数向量，绘制样本数据的95%置信区间。

x = normrnd(10,4,100,1);

subplot(1,2,1)

ConfidenceRegion(x)

subplot(1,2,2)

ConfidenceRegion(x,\'exp\')

效果图如下图所示：

【例2】绘制置信椭圆。产生二元正态分布随机数矩阵，绘制样本数据的95%置信椭圆区域。

x = mvnrnd([1;2],[1 4;4 25],100);

subplot(1,2,1)

ConfidenceRegion(x)

subplot(1,2,2)

ConfidenceRegion(x,\'exp\')

果图如下图所示：

【例3】绘制置信椭球。产生三元正态分布随机数矩阵，绘制样本数据的95%置信椭球区域。

x = mvnrnd([1;2;3],[3 0 0;0 5 -1;0 -1 1],100);

subplot(1,2,1)

ConfidenceRegion(x)

subplot(1,2,2)

ConfidenceRegion(x,\'exp\')

效果图如下图所示：

自编函数ConfidenceRegion程序代码如下：

function ConfidenceRegion(xdat,alpha,distribution)

% 绘制置信区域（区间、椭圆区域或椭球区域）

% ConfidenceRegion(xdat,alpha,distribution)

% xdat：样本观测值矩阵,p*N 或 N*p的矩阵，p = 1,2 或 3

% alpha：显著性水平，标量，取值介于[0,1]，默认值为0.05

% distribution：字符串（\'norm\'或\'experience\'），用来指明求置信区域用到的分布类型，

% distribution的取值只能为字符串\'norm\'或\'experience\'，分别对应正态分布和经验分布

% CopyRight：xiezhh（谢中华）

% date：2011.4.14

%

% example1：x = normrnd(10,4,100,1);

% ConfidenceRegion(x)

% ConfidenceRegion(x,\'exp\')

% example2：x = mvnrnd([1;2],[1 4;4 25],100);

% ConfidenceRegion(x)

% ConfidenceRegion(x,\'exp\')

% example3：x = mvnrnd([1;2;3],[3 0 0;0 5 -1;0 -1 1],100);

% ConfidenceRegion(x)

% ConfidenceRegion(x,\'exp\')

% 设定参数默认值

if nargin == 1

distribution = \'norm\';

alpha = 0.05;

elseif nargin == 2

if ischar(alpha)

distribution = alpha;

alpha = 0.05;

else

distribution = \'norm\';

end

end

% 判断参数取值是否合适

if ~isscalar(alpha) || alpha>=1 || alpha<=0

error(\'alpha的取值介于0和1之间\')

end

if ~strncmpi(distribution,\'norm\',3) && ~strncmpi(distribution,\'experience\',3)

error(\'分布类型只能是正态分布（\'\'norm\'\'）或经验分布（\'\'experience\'\'）\')

end

% 检查数据维数是否正确

[m,n] = size(xdat);

p = min(m,n); % 维数

if ~ismember(p,[1,2,3])

error(\'应输入一维、二维或三维样本数据,并且样本容量应大于3\')

end

% 把样本观测值矩阵转置，使得行对应观测，列对应变量

if m < n

xdat = xdat\';

end

xm = mean(xdat); % 均值

n = max(m,n); % 观测组数

% 分情况绘制置信区域

switch p

case 1 % 一维情形（置信区间）

xstd = std(xdat); % 标准差

if strncmpi(distribution,\'norm\',3)

lo = xm - xstd*norminv(1-alpha/2,0,1); % 正态分布情形置信下限

up = xm + xstd*norminv(1-alpha/2,0,1); % 正态分布情形置信上限

%lo = xm - xstd*tinv(1-alpha/2,n-1); % 正态分布情形置信下限

%up = xm + xstd*tinv(1-alpha/2,n-1); % 正态分布情形置信上限

TitleText = \'置信区间（基于正态分布）\';

else

lo = prctile(xdat,100*alpha/2); % 经验分布情形置信下限

up = prctile(xdat,100*(1-alpha/2)); % 经验分布情形置信上限

TitleText = \'置信区间（基于经验分布）\';

end

% 对落入区间内外不同点用不同颜色和符号绘图

xin = xdat(xdat>=lo & xdat<=up);

xid = find(xdat>=lo & xdat<=up);

plot(xid,xin,\'.\')

hold on

xout = xdat(xdat<lo | xdat>up);

xid = find(xdat<lo | xdat>up);

plot(xid,xout,\'r+\')

h = refline([0,lo]);

set(h,\'color\',\'k\',\'linewidth\',2)

h = refline([0,up]);

set(h,\'color\',\'k\',\'linewidth\',2)

xr = xlim;

yr = ylim;

text(xr(1)+range(xr)/20,lo-range(yr)/20,\'置信下限\',...

\'color\',\'g\',\'FontSize\',15,\'FontWeight\',\'bold\')

text(xr(1)+range(xr)/20,up+range(yr)/20,\'置信上限\',...

\'color\',\'g\',\'FontSize\',15,\'FontWeight\',\'bold\')

xlabel(\'观测序号\')

ylabel(\'观测值\')

title(TitleText)

hold off

case 2 % 二维情形（置信椭圆）

x = xdat(:,1);

y = xdat(:,2);

s = inv(cov(xdat)); % 协方差矩阵

xd = xdat-repmat(xm,[n,1]);

rd = sum(xd*s.*xd,2);

if strncmpi(distribution,\'norm\',3)

r = chi2inv(1-alpha,p);

%r = p*(n-1)*finv(1-alpha,p,n-p)/(n-p)/n;

TitleText = \'置信椭圆（基于正态分布）\';

else

r = prctile(rd,100*(1-alpha));

TitleText = \'置信椭圆（基于经验分布）\';

end

plot(x(rd<=r),y(rd<=r),\'.\',\'MarkerSize\',16) % 画样本散点

hold on

plot(x(rd>r),y(rd>r),\'r+\',\'MarkerSize\',10) % 画样本散点

plot(xm(1),xm(2),\'k+\'); % 椭圆中心

h = ellipsefig(xm,s,r,1); % 绘制置信椭圆

xlabel(\'X\')

ylabel(\'Y\')

title(TitleText)

hold off;

case 3 % 三维情形（置信椭球）

x = xdat(:,1);

y = xdat(:,2);

z = xdat(:,3);

s = inv(cov(xdat)); % 协方差矩阵

xd = xdat-repmat(xm,[n,1]);

rd = sum(xd*s.*xd,2);

if strncmpi(distribution,\'norm\',3)

r = chi2inv(1-alpha,p);

%r = p*(n-1)*finv(1-alpha,p,n-p)/(n-p)/n;

TitleText = \'置信椭球（基于正态分布）\';

else

r = prctile(rd,100*(1-alpha));

TitleText = \'置信椭球（基于经验分布）\';

end

plot3(x(rd<=r),y(rd<=r),z(rd<=r),\'.\',\'MarkerSize\',16) % 画样本散点

hold on

plot3(x(rd>r),y(rd>r),z(rd>r),\'r+\',\'MarkerSize\',10) % 画样本散点

plot3(xm(1),xm(2),xm(3),\'k+\'); % 椭球中心

h = ellipsefig(xm,s,r,2); % 绘制置信椭球

xlabel(\'X\')

ylabel(\'Y\')

zlabel(\'Z\')

title(TitleText)

hidden off;

hold off;

end

%--------------------------------------------------

% 子函数：用来绘制置信区域（椭圆或椭球）

%--------------------------------------------------

function h = ellipsefig(xc,P,r,tag)

% 画一般椭圆或椭球：(x-xc)\'*P*(x-xc) = r

[V, D] = eig(P);

if tag == 1

aa = sqrt(r/D(1));

bb = sqrt(r/D(4));

t = linspace(0, 2*pi, 60);

xy = V*[aa*cos(t);bb*sin(t)]; % 坐标旋转

h = plot(xy(1,:)+xc(1),xy(2,:)+xc(2), \'k\', \'linewidth\', 2);

else

aa = sqrt(r/D(1,1));

bb = sqrt(r/D(2,2));

cc = sqrt(r/D(3,3));

[u,v] = meshgrid(linspace(-pi,pi,30),linspace(0,2*pi,30));

x = aa*cos(u).*cos(v);

y = bb*cos(u).*sin(v);

z = cc*sin(u);

xyz = V*[x(:)\';y(:)\';z(:)\']; % 坐标旋转

x = reshape(xyz(1,:),size(x))+xc(1);

y = reshape(xyz(2,:),size(y))+xc(2);

z = reshape(xyz(3,:),size(z))+xc(3);

h = mesh(x,y,z); % 绘制椭球面网格图

end