Matlab之检验假设

Size: px

Start display at page:

Download "Matlab之检验假设"

畏懿邮成
7 years ago
Views:

1 Matlab 之检验假设专业 : 天体物理姓名 : 聂俊丹学号 : 在统计中常见的是 : 需要多大的样本? 这是我们很关心的一个问题在 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

2 进行双边检验时, 假设原假设错误的机率是 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= ,cutoff2=2.131 如果统计量 t 属于拒绝域 : t< cutoff1 或 t>cutoff2, 就拒绝原假设, 而进行备择假设如何进行备择假设? 一般是计算备择假设的功效函数功效函数 (power 2

3 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

4 函数说明 : 从图中可以看到在 =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

5 怎样才知道我们做得对不对, 我们要进行进一步的验证 : 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

6 程序说明 :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= 检验相关性 Sampsizepwr 函数不是在任何情况下都可以用来计算样本数, 当不是双边检验时怎样来计算样本大小? 不同的检验类型有不同的方法以下要介绍的方法是通过不断试探来确定样本数这里以检验相关系数为例假设可以从 X,Y 变量中抽取样本, 而且又想知道在备择假设的情况下, 相关系数为.4 时, 这两个抽取的样本应取多大 ( 因为备择下相关系数是.4, 基本认为是相关的如果备择假设正确而拒绝原假设 (power), 就可以知道原假设是不相关的, 所以才拒绝 ) 现在的一个问题就是要多大的样本, 才能使原假设两样本不相关显然 power 越大, 他们越不相关以下是实现的步骤 : 6

7 当给定一个样本时, 可以用 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 = 可以近似地认为 power 是 65%, 有 95% 的把握认为真实值的 power 在 62%~ 68% 之间所以为了得到 8% 的, 需要更大的样本可以尝试样本为 5 然后重复以上的程序, 看运行后的 power 为多大 nsamples = 1; N = 5; alpha =.5; conf = 1-alpha; 7

8 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 = power 值偏大, 那么需要减少样本, 以达到更接近 8% 不断的试探, 直到得理想的结果这种方法的特点就是要不断地试探结果, 最终得到理想的样本数 3. 结论计算双边检验的样本大小可以用 sampsizepwr 函数来实现,sampsizepwr 也可以用来计算备择假设的功效值当不是双边检验时, 可以采用其它方法计算样本数, 比如检验相关性的时候可以通过不断的试探得到满意的样本数 8

p-2.indd

p-2.indd 2 第一篇睽違 1342 年, 陳元光終於回到老家開漳聖王陳元光回老家活動, 是開漳聖王聖像自台灣啟程回固始老家, 讓 1342 年自固始出發前往閩南漳洲, 其後英靈飄洋渡海來到台灣的開漳聖王, 再回到原始出生地固始, 形成一個圓滿的迴