prob_check.m

matlab从入门到精通第6章书后源码对初学者很有帮助

% prob_check.m
% 假设检验演示程序
% 单侧假设检验 %
% 设置检验参数
mu0=100;sig=20;N=16;
% 设置假设检验的置信水平
alpha=0.05;conf=1-alpha;
% 设置正态分布的截断点
cutoff=norminv(conf,mu0,sig/sqrt(N));
% 产生数据点
x=[linspace(90,cutoff),linspace(cutoff,127)];
y=normpdf(x,mu0,sig/sqrt(N));
% 绘制正态分布图形
h1=plot(x,y);
xhi=[cutoff,x(x>=cutoff)];
yhi=[0,y(x>=cutoff)];
patch(xhi,yhi,'b');
% 添加注释
title('distribution hypothesis, N=16');
xlabel('sample mean');ylabel('density');
text(96,0.01,sprintf('rejected if mean>%.4g\nProb=0.05',cutoff),'color','b');

% 修改检验假设的数据 %
mu1=110;
y2=normpdf(x,mu1,sig/sqrt(N));
% 绘制图形
h2=line(x,y2,'color','r');
% 绘制假设检验的面积图
yhi=[0,y2(x>=cutoff)];
patch(xhi,yhi,'r','FaceAlpha',0.25);
% 添加图形注释
P=1-normcdf(cutoff,mu1,sig/sqrt(N));
text(115,0.06,sprintf('rejected if T>%.4g\nProb=%.2g',cutoff,P),'color',[1 0 0]);
legend([h1 h2],'null hypothesis','another hypothesis');axis tight;

% 计算累计密度函数图形
ynull=normcdf(x,mu0,sig/sqrt(N));
yalt=normcdf(x,mu1,sig/sqrt(N));
% 绘制累计密度函数图形
figure;
plot(x,ynull,'b-',x,yalt,'r-');
% 计算置信条件水平下的反正态分布数值
zval=norminv(conf);
cutoff=mu0+zval*sig/sqrt(N);
% 绘制图形
line([90,cutoff,cutoff],[conf,conf,0],'linestyle',':');
msg=sprintf('cutoff=\\mu_0+%.2g\\sigma/\\surd{n}',zval);
text(cutoff,0.15,msg,'color','b');
text(min(x),conf,sprintf('%g%%test',100*alpha),'color','b',...
    'verticalalignment','top');
palt=normcdf(cutoff,mu1,sig/sqrt(N));
line([90,cutoff],[palt,palt],'color','r','linestyle',':');
text(91,palt+0.02,sprintf('power is 1-%.2g=%.2g',palt,1-palt),'color',[1 0 0]);

% 定义power需要的参数数值
desiredPower=0.80;
Nvec=1:30;
cutoff=mu0+norminv(conf)*sig./sqrt(Nvec);
% 计算假设检验的power数值
power=1-normcdf(cutoff,mu1,sig./sqrt(Nvec));
% 绘制图形
figure;
plot(Nvec,power,'bo-',[0 30],[desiredPower desiredPower],'k:');
xlabel('N=sample size');ylabel('Power')
title('power for the alternate hypothesis: \mu=110');

% 使用Monte Carlo方法检验power函数的结果
% 定义Monte Carlo模拟的参数
nsamples=400;
samplenum=1:nsamples;
N=25;
% 创建零值矩阵
hh0=zeros(1,nsamples);
hh1=hh0;
% 进行右侧已知方差条件下均值假设检验
for j=1:nsamples
    Z0=normrnd(mu0,sig,N,1);
    hh0(j)=ztest(Z0,mu0,sig,alpha,'right');
    Z1=normrnd(mu1,sig,N,1);
    hh1(j)=ztest(Z1,mu0,sig,alpha,'right');
end
p0=cumsum(hh0)./samplenum;
p1=cumsum(hh1)./samplenum;
% 绘制对应的图形
figure;
plot(samplenum,p0,'b-',samplenum,p1,'r-');
xlabel('sample number');
ylabel('proportion significant');
title('verification of power computation');
legend('null hypothesis','another hypothesis');