Matlab之检验假设

Matlab 之检验假设 专业 : 天体物理姓名 : 聂俊丹学号 :712162 在统计中常见的是 : 需要多大的样本? 这是我们很关心的一个问题 在 matlab 统计工具箱中有一个函数 :sampsizepwr 可以用来计算样本大小 这篇论文的目的就是阐述如何来使用这个函数 文章中通过特殊的例子来实现具体的计算过程 同时 sampsizepwr 这个函数还有其它的功能 : 可以用来计算功效 在本文中也具体介绍了如何用 sampsizepwr 来计算功效函数值 除此之外, 我们还列举了一些其它的例子 当 sampsizepwr 函数不能使用的情况下如何来确定样本大小 1. sampsizepwr 函数计算样本数及 power 值 Sampsizepwr 函数可以用来计算双边检验的样本大小和 power 值 但 sampsizepwr 函数不是在任何情况下都可以使用的, 它只能用在假设检验中 假设检验有两种情况 : 一种是单边检验, 一种是双边检验 Sampsizepwr 在双边检验中用得比较多 当不知道标准偏差的情况下进行均值检验, 可以采用双边检验 所谓双边检验是 : 在原假设不成立的情况下进行备择检验, 不管样本均值是偏大还是偏小 即 : H 1 : u = u, H : u u. 其中代表原假设, H 代表备择假设 在这种检验中, 统计量是 t 统计量, H 1 它服从 : u u t ~ δ x 在原假设下,t 服从学生式 t 分布, 具有 N 1个自由度 ; 而在备择检验的情 况下它是一个有偏的统计量, 而且这个有偏的参数的值为真实值与检验均值的标 准差 顺便提及下单边检验, 它的具体形式是 : H : u = u, H 1 : u > u或 u < u 1

进行双边检验时, 假设原假设错误的机率是 5%( 显著水平 ) 如果原假设的统计量属于拒绝域, 就拒绝原假设, 在备择假设下进行双边检验 下面的这个程序是进行双边检验的具体实现步骤 : N = 16; df = N-1; alpha =.5; conf = 1-alpha; cutoff1 = tinv(alpha/2,df); cutoff2 = tinv(1-alpha/2,df); x = [linspace(-5,cutoff1), linspace(cutoff1,cutoff2),linspace(cutoff2,5)]; y = tpdf(x,df); h1 = plot(x,y); xlo = [x(x<=cutoff1),cutoff1]; ylo = [y(x<=cutoff1),]; xhi = [cutoff2,x(x>=cutoff2)]; yhi = [, y(x>=cutoff2)]; patch(xlo,ylo,'b'); patch(xhi,yhi,'b'); title('distribution of t statistic, N=16'); xlabel('t'); ylabel('density'); text(2.5,.5,sprintf('reject if t>%.4g\nprob =.25',cutoff2),'Color','b'); text(-4.5,.5,sprintf('reject if t<%.4g\nprob =.25',cutoff1),'Color','b'); 程序说明 : 自由度是 N = 16, 显著水平是.5,cutoff1 和 cutoff2 是拒绝域的临界值 通过运行程序, 得到 cutoff1= -2.131,cutoff2=2.131 如果统计量 t 属于拒绝域 : t< cutoff1 或 t>cutoff2, 就拒绝原假设, 而进行备择假设 如何进行备择假设? 一般是计算备择假设的功效函数 功效函数 (power 2

function) 的定义 : 备择假设成立的情况下拒绝原假设的概率 它的值取决于备择假设的 u 值, 以及样本的大小 计算功效函数的目的也就是想知道备择假设成立的概率 一般不考虑原假设及 u 值下的功效函数, 而是把功效函数看成是一个 u 关于的函数 u 越远离原假设的 u, 功效函数的值 (power) 的值越大 power 值可以通过 sampsizepwr 函数计算得到 Sampsizepwr 函数的形式有以下几种 : n = sampsizepwr(testtype,p,p1) N=sampsizepwr(testtype,p,p1,power) power = sampsizepwr(testtype,p,p1,[ ],n) p1 = sampsizepwr(testtype,p,[ ],power,n) 函数说明 :testtype 是指统计量的类型, 比如 :t 统计量,z 统计量,p 统计量等 p 是 原假设下参数的值 ( 如 ),p1 是备择假设下的参数值 ( 如, 或 u ) u u1 2 n = sampsizepwr(testtype,p,p1) 返回值 n 为样本大小, 默认的功效值为.9, 置信水平为.5 N=sampsizepwr(testtype,p,p1,power) 当已知 power 值时返回样本数 N power = sampsizepwr(testtype,p,p1,[ ],n) 当已知样本数 n 时返回 power 值 p1 = sampsizepwr(testtype,p,[ ],power,n) 当已知 power 和 n 时返回备择假设的参数值 同样也可以求 p,p2 下面介绍在程序中怎么来使用 sampsizepwr 函数来计算 power: 假设标准偏差大约为 2, 样本数为 16,u 取值范围 :9 127 N = 16; x = linspace(9,127); power = sampsizepwr('t',[1 2],x,[],N); plot(x,power); xlabel('true mean') ylabel('power') title('power function for N=16') 3

函数说明 : 从图中可以看到在 =1 的两边, 越远离 u, u 对应的 power 的值越大 u 当 N=16, u =11 时 power 的值为 5% 左右 假如想要得到 power=8%,u =11, 显然 N=16 是不够的 那么要多大的样本才符合要求呢? 从 sampersizepwr 函数的几种形式中看到, 可以使用 N=sampsizepwr(testtype,p,p1,power) 来计算样本大小 : N=sampsizepwr( t,[1 2],11,.8), 得到 N=34 这样的话就需要在原样本数上每边加上 9 的样本数 假如想得到不同 power 下的 N 值, 可以 : Nvec = 2:4; power = sampsizepwr('t',[1 2],11,[],Nvec); plot(nvec,power,'bo-',[ 4],[DesiredPower DesiredPower], 'k:'); xlabel('n = sample size') ylabel('power') title('power function for the alternative: \mu = 11') 4

怎样才知道我们做得对不对, 我们要进行进一步的验证 : nsamples = 4; samplenum = 1:nsamples; N = 34; h = zeros(1,nsamples); h1 = h; Z = normrnd(mu,sig,n,1); h(j) = ttest(z,mu,alpha); Z1 = normrnd(mu1,sig,n,1); h1(j) = ttest(z1,mu,alpha); p = cumsum(h)./ samplenum; p1 = cumsum(h1)./ samplenum; plot(samplenum,p,'b-',samplenum,p1,'r-') xlabel('sample number'); ylabel('proportion significant') title('verification of power computation') leg('null hypothesis','alternative hypothesis','location','east') 5

程序说明 :mu 是原假设的值,mu1 是备择假设的一个值 ttest 是 t 检验, 它的返回值是 : 当置信度为 5% 时, 拒绝原假设则返回值是 h=1, 不拒绝原假设则 h= P 时原假设的置信水平,p= 累计值 / 样本数 由图可以看出, 原假设的置信水平大约为.5,power 的值大约为.8 这正是我们原来程序所计算出来的 通过这个程序我们验证了备择假设是对的 如果样本数为 N=5, 那么在 u =11 时,power 肯定大于 8%( 因为原先认为样本数为 34 时,power 正好是 8%) 那么要达到 N=5,power=8%, 备择参数值 mu1 应为多少? 可以通过下面的语句实现 : mu1 = sampsizepwr('t',[1 2],[ ],.8,5) 得到 mu1=18.837 2. 检验相关性 Sampsizepwr 函数不是在任何情况下都可以用来计算样本数, 当不是双边检验时怎样来计算样本大小? 不同的检验类型有不同的方法 以下要介绍的方法是通过不断试探来确定样本数 这里以检验相关系数为例 假设可以从 X,Y 变量中抽取样本, 而且又想知道在备择假设的情况下, 相关系数为.4 时, 这两个抽取的样本应取多大 ( 因为备择下相关系数是.4, 基本认为是相关的 如果备择假设正确而拒绝原假设 (power), 就可以知道原假设是不相关的, 所以才拒绝 ) 现在的一个问题就是要多大的样本, 才能使原假设两样本不相关 显然 power 越大, 他们越不相关 以下是实现的步骤 : 6

当给定一个样本时, 可以用 Monte carol 方法来确定相关检验的一个边界值 cotoff 然后重复运行多次, 可以得到准确的 cutoff 假设样本是 25: nsamples = 1; N = 25; alpha =.5; conf = 1-alpha; r = zeros(1,nsamples); xy = normrnd(,1,n,2); r(j) = corr(xy(:,1),xy(:,2)); cutoff = quantile(r,conf) 运行后得 : Cutoff=.3382 这样的话, 找到了准确的 cutoff 然后可以在备择下产生样本, 并且计算 power 值 : nsamples = 1; mu = [; ]; sig = [1.4;.4 1]; r = zeros(1,nsamples); xy = mvnrnd(mu,sig,n); r(j) = corr(xy(:,1),xy(:,2)); [power,powerci] = binofit(sum(r>cutoff),nsamples) 程序说明 :phat=binofit(x,n) 函数 二项式参数估计, 返回值是 : 在 N 次独立试验中, 成功地次数为 x 时返回成功概率的极大似然估计 如果 x 为一个向量的话 (x(1)-x(k)), 那么返回的是 k 的极大似然估计值 每次试验都是相互独立的 [phat, pi]=binofit(x,n) 返回的是概率估计 phat, 以及 95% 的置信区间 pci 和其它分布拟合函数不同的是,x 是指从相互独立样本中取的数, 假如把 x 看成是单个样本, 而且产生一个参数估计的话, 就可以用 binofit(sum(x),n) 来实现 运行后 :power =.649 powerci =.6185.6786 可以近似地认为 power 是 65%, 有 95% 的把握认为真实值的 power 在 62%~ 68% 之间 所以为了得到 8% 的, 需要更大的样本 可以尝试样本为 5 然后重复以上的程序, 看运行后的 power 为多大 nsamples = 1; N = 5; alpha =.5; conf = 1-alpha; 7

r = zeros(1,nsamples); xy = normrnd(,1,n,2); r(j) = corr(xy(:,1),xy(:,2)); cutoff = quantile(r,conf) nsamples = 1; mu = [; ]; sig = [1.4;.4 1]; r = zeros(1,nsamples); xy = mvnrnd(mu,sig,n); r(j) = corr(xy(:,1),xy(:,2)); [power,powerci] = binofit(sum(r>cutoff),nsamples) 运行后 :cutoff =.2387 power =.899 powerci =.8786.917 power 值偏大, 那么需要减少样本, 以达到更接近 8% 不断的试探, 直到得理想的结果 这种方法的特点就是要不断地试探结果, 最终得到理想的样本数 3. 结论 计算双边检验的样本大小可以用 sampsizepwr 函数来实现,sampsizepwr 也可 以用来计算备择假设的功效值 当不是双边检验时, 可以采用其它方法计算样本 数, 比如检验相关性的时候可以通过不断的试探得到满意的样本数 8