<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
