| 查看: 233 | 回复: 0 | |||
[交流]
【求助】求一统计方面的MATLAB高手帮我解决，谢谢，急
|
|
求一统计方面的MATLAB高手帮我解决，谢谢，急 function fdr=rand01(); %向量a，b，为随机数构成 的向量 %向量 c，为a和b的组合， %mark为一向量，标记c中的元素是a中的元素还是b中的元素 lvv=[]; lm=[]; t=[]; lvvb=[]; lmb=[]; tb=[]; newfdr=[]; BHfdr=[]; newpower=[]; BHpower=[]; sk=0; N=500 for n1=100:100:N q=0.05; m1=100; m2=100; n2=N-n1; for l=1:1000 a=mean(normrnd(0,1,m1,n1)); p1=[]; for j=1:n1 p1(j) =1-normcdf(a(1,j)*(sqrt(m1)),0,1); end %抽取真实原假设的样本，并计算P值 if n2 > 0 b=mean(normrnd(1,1,m2,n2)); p2=[]; for k=1:n2 p2(k) =1-normcdf(b(1,k)*(sqrt(m2)),0,1); end %抽取错误原假设的样本，并计算P值 c=[p1,p2]; mark=[zeros(1,n1),ones(1,n2)]; for i=1n1+n2-1) for j=i+1n1+n2) if c(i)>c(j) t1=c(i); c(i)=c(j); c(j)=t1; t2=mark(i); mark(i)=mark(j);%排序 mark(j)=t2; end end end r1=[]; for m=1n1+n2) r1(m)=m*q/(n1+n2); end s=n1+n2; j1=0; jj1=0; while s>0 && c(s)> r1(s) if mark(s)==0 j1=j1+1; end if mark(s)==1 jj1=jj1+1; end s=s-1; end totalb=s ; %拒绝的总个数 ttb=n1-j1 ; %错误拒绝的个数 total2b=n2-jj1 ; %正确拒绝的个数 tb(l)=total2b/n2 ; %功效 if totalb == 0 lvvb(l)=0; else lvvb(l)=ttb/totalb ; %V/R end end %第一个while语句结束 if n2==0 r1=[]; for m=1:n1 r1(m)=m*q/(n1+n2); end d=[p1]; c=sort(d); ig=n1; while ig>0 && c(ig)>r1(ig) ig=ig-1; end jg=0; if ig>0 jg=jg+1; end lvvb(l)=jg ; tb(l)=0; end %第二个if语句结束 end sk=sk+1 newpower(sk)=mean(tb); newfdr(sk)=mean(lvvb); end newpower newfdr x=[100:100:N]; subplot(2,1,1) plot(x,newfdr,'-k'); axis ([0,N,0,0.1]); legend('newfdr'); ylabel('FDR');xlabel('真实原假设个数'); grid on subplot(2,1,2) function fdr=rand01(); %向量a，b，为随机数构成 的向量 %向量 c，为a和b的组合， %mark为一向量，标记c中的元素是a中的元素还是b中的元素 lvv=[]; lm=[]; t=[]; lvvb=[]; lmb=[]; tb=[]; newfdr=[]; BHfdr=[]; newpower=[]; BHpower=[]; sk=0; N=500 for n1=100:100:N q=0.05; m1=100; m2=100; n2=N-n1; for l=1:1000 a=mean(normrnd(0,1,m1,n1)); p1=[]; for j=1:n1 p1(j) =1-normcdf(a(1,j)*(sqrt(m1)),0,1); end %抽取真实原假设的样本，并计算P值 if n2 > 0 b=mean(normrnd(1,1,m2,n2)); p2=[]; for k=1:n2 p2(k) =1-normcdf(b(1,k)*(sqrt(m2)),0,1); end %抽取错误原假设的样本，并计算P值 c=[p1,p2]; mark=[zeros(1,n1),ones(1,n2)]; for i=1n1+n2-1) for j=i+1n1+n2) if c(i)>c(j) t1=c(i); c(i)=c(j); c(j)=t1; t2=mark(i); mark(i)=mark(j);%排序 mark(j)=t2; end end end r1=[]; for m=1n1+n2) r1(m)=m*q/(n1+n2); end s=n1+n2; j1=0; jj1=0; while s>0 && c(s)> r1(s) if mark(s)==0 j1=j1+1; end if mark(s)==1 jj1=jj1+1; end s=s-1; end totalb=s ; %拒绝的总个数 ttb=n1-j1 ; %错误拒绝的个数 total2b=n2-jj1 ; %正确拒绝的个数 tb(l)=total2b/n2 ; %功效 if totalb == 0 lvvb(l)=0; else lvvb(l)=ttb/totalb ; %V/R end end %第一个while语句结束 if n2==0 r1=[]; for m=1:n1 r1(m)=m*q/(n1+n2); end d=[p1]; c=sort(d); ig=n1; while ig>0 && c(ig)>r1(ig) ig=ig-1; end jg=0; if ig>0 jg=jg+1; end lvvb(l)=jg ; tb(l)=0; end %第二个if语句结束 end sk=sk+1 newpower(sk)=mean(tb); newfdr(sk)=mean(lvvb); end newpower newfdr x=[100:100:N]; subplot(2,1,1) plot(x,newfdr,'-k'); axis ([0,N,0,0.1]); legend('newfdr'); ylabel('FDR');xlabel('真实原假设个数'); grid on subplot(2,1,2) plot(x,newpower,'-k'); axis ([0,N,0,1]); legend('newpower'); ylabel('功效');xlabel('真实原假设个数'); grid on 上面这个程序是多重假设检验的内容，希望有高手帮我看看这个程序有没有问题，我总觉得不对劲 多重假设检验中，真实原假设服从N(0,1)，错误原假设服从N（1,1），检验的总个数是N=500，所以随机抽取了500组数，其中有n1组来自N(0,1)，n2组来自N（1,1），取每组的均值计算的P值，样本量是100，谢谢帮我看看， |
回复此楼» 猜你喜欢
假如你的研究生提出不合理要求
已经有6人回复
-
对氯苯硼酸纯化
已经有3人回复
-
求助:我三月中下旬出站，青基依托单位怎么办？
已经有12人回复
-
不自信的我
已经有12人回复
-
所感
已经有4人回复
-
论文终于录用啦！满足毕业条件了
已经有28人回复
-
要不要辞职读博？
已经有7人回复
-
北核录用
已经有3人回复
-
实验室接单子
已经有3人回复
-
磺酰氟产物，毕不了业了！
已经有8人回复