24小时热门版块排行榜    

查看: 1369  |  回复: 4
本帖产生 1 个 程序强帖 ,点击这里进行查看

ian_zhangty

木虫 (著名写手)

[求助] 请问,我想用fortran计算统计中的p值以及95%信度空间

请问有没有现成的程序,谢谢
回复此楼
你打你的,我打我的
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

snoopyzhao

至尊木虫 (职业作家)

【答案】应助回帖


jjdg(金币+1): 感谢您的关注 2011-08-06 01:01:29
ian_zhangty(金币+1): 2011-08-06 09:10:51
计算 p 值可不是简单的事情,呵呵……
去 statlib 和 netlib 上找找

其实,这种事情,找个统计软件做了就结了……
2楼2011-08-05 23:00:50
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ian_zhangty

木虫 (著名写手)


dubo(金币+1): 欢迎常来程序语言版讨论 2011-08-21 14:00:44
引用回帖:
2楼: Originally posted by snoopyzhao at 2011-08-05 23:00:50:
计算 p 值可不是简单的事情,呵呵……
去 statlib 和 netlib 上找找

其实,这种事情,找个统计软件做了就结了……

我也知道的。但是由于我的整个程序都是用fortran编的,如果有直接的程序就非常省事了
你打你的,我打我的
3楼2011-08-06 09:10:28
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

snoopyzhao

至尊木虫 (职业作家)

【答案】应助回帖

ian_zhangty(金币+4): 2011-08-06 10:36:46
那就上我说的那两个上面找找,也可以去看看 R (http://www.r-project.org) 或者 dataplot (http://www.itl.nist.gov/div898/software/dataplot/) 的源码中关于 p 的计算。后者主要是 fortran 编的,前者有一部分是 C,有一部分是 fortran,还有大部分是 R....
4楼2011-08-06 09:32:22
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ian_zhangty

木虫 (著名写手)

【答案】应助回帖

★ ★ ★ ★ ★ ★
jjdg(金币+1): 感谢分享你的解决方案 2011-08-06 15:28:45
dubo(金币+5, 程序强帖+1): good 2011-08-21 14:01:09
已经解决,方案如下,谢谢

FUNCTION betai(a,b,x)
REAL betai,a,b,x
!USES betacf,gammln
REAL bt,betacf,gammln
if(x<0..or.x>1.) pause 'bad argument x in betai'
if(x==0..or.x==1.) then
  bt=0.
else
  bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)&
     +b*log(1.-x))
endif
if(x<(a+1.)/(a+b+2.)) then
  betai=bt*betacf(a,b,x)/a
  return
else
  betai=1.-bt*betacf(b,a,1.-x)/b
  return
endif
END FUNCTION betai

FUNCTION betacf(a,b,x)
INTEGER maxit
REAL betacf,a,b,x,EPS,fpmin
PARAMETER (maxit=100,EPS=3.e-7,fpmin=1.e-30)
INTEGER m,m2
REAL aa,c,d,del,h,qab,qam,qap
qab=a+b
qap=a+1.
qam=a-1.
c=1.
d=1.-qab*x/qap
if(abs(d) d=1./d
h=d
do m=1,maxit
  m2=2*m
  aa=m*(b-m)*x/((qam+m2)*(a+m2))
  d=1.+aa*d
  if(abs(d)   c=1.+aa/c
  if(abs(c)   d=1./d
  h=h*d*c
  aa=-(a+m)*(qab+m)*x/((a+m2)*(qap+m2))
  d=1.+aa*d
  if(abs(d)   c=1.+aa/c
  if(abs(c)   d=1./d
  del=d*c
  h=h*del
  if(abs(del-1.)     betacf=h
    return
  end if
end do
pause 'a or b too big, or maxit too small in betacf'
END FUNCTION betacf

FUNCTION gammln(xx)
REAL gammln,xx
INTEGER j
DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
SAVE cof,stp
DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,&
     24.01409824083091d0,-1.231739572450155d0,&
         .1208650973866179d-2,-.5395239384953d-5,&
         2.5066282746310005d0/
x=xx
y=x
tmp=x+5.5d0
tmp=(x+0.5d0)*log(tmp)-tmp
ser=1.000000000190015d0
do j=1,6
  y=y+1.d0
  ser=ser+cof(j)/y
end do
gammln=tmp+log(stp*ser/x)
END FUNCTION gammln



program pvalue
        REAL df, r, t, prob
        df=17-2   ! 17 is sample size

        r=0.188446645   ! pearson coefficient
        t=(abs(r)*sqrt(df))/sqrt(1-r**2)
        prob=betai(0.5*df,0.5,df/(df+t**2))
        write (*,*) prob





end program pvalue
你打你的,我打我的
5楼2011-08-06 10:41:54
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 瘦竹云淡风轻 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 时间戳变了,能看出什么问题? +7 基诺咪客 2026-08-17 10/500 2026-08-18 16:24 by Lullabygxb
[基金申请] 重要消息,中午系统在维护 +7 yuleib84 2026-08-18 8/400 2026-08-18 16:05 by jnhyjjm
[基金申请] 欢迎发来filecode的Mz6后的代码验证其规律 +47 医学老男孩 2026-08-13 105/5250 2026-08-18 16:05 by 医学老男孩
[基金申请] 明天放榜? +4 Shxjjxjkx 2026-08-18 4/200 2026-08-18 15:59 by abcabc1133
[论文投稿] 投稿咨询 +4 wwm09 2026-08-17 6/300 2026-08-18 15:36 by wwm09
[教师之家] 两姐妹,一个辛辛苦苦30岁工作当大学教师,另一个23岁就工作当开车教练 +3 瞬息宇宙 2026-08-11 5/250 2026-08-18 14:55 by LNP@mRNA
[基金申请] 时间戳又变了8-15 +14 archvillain 2026-08-15 26/1300 2026-08-18 13:37 by phantomgost
[基金申请] 哪位老哥知道今年的国自然具体哪一天放榜? +13 Ldrop2023 2026-08-13 16/800 2026-08-18 12:25 by 淀粉搬运工
[基金申请] 2027广东省杰青 +3 奶牛小黑 2026-08-15 8/400 2026-08-18 11:47 by gltch
[基金申请] 快农历七夕节了,轻松一下,男人悄悄话,女施主请不要进来。 +4 Tide man 2026-08-14 5/250 2026-08-18 11:35 by Tide man
[基金申请] 93BebMhtakh前后11位开头都是大写 +4 且听虎啸 2026-08-17 5/250 2026-08-18 00:49 by 蔡棒棒菂
[基金申请] 感觉是下周放榜了 +6 angus9576 2026-08-17 11/550 2026-08-17 23:57 by angus9576
[基金申请] filecode=后面第一个是大写字母 +8 wangze12014 2026-08-14 10/500 2026-08-17 17:05 by xter9665
[基金申请] 今天系统多次维护,明天很可能放榜! +8 zju2000 2026-08-16 9/450 2026-08-17 12:20 by lmz0216
[基金申请] 咱们一起用铁证分析2026国家社科基金中标与否 +7 启萌科技 2026-08-12 26/1300 2026-08-16 12:35 by 启萌科技
[精细化工] 招聘 金属平磨液,抛光液研发工程师 +3 小天0311 2026-08-14 3/150 2026-08-16 07:31 by H9PLUS
[基金申请] 各位道友,我要去昆明玩几天,回来见。 +7 Tide man 2026-08-14 8/400 2026-08-15 01:11 by arzu_hma
[基金申请] 是这周出结果还是下周出结果? +4 yuleib84 2026-08-11 4/200 2026-08-14 23:05 by lfy8008
[基金申请] 不应该看fileCode +7 且听虎啸 2026-08-12 9/450 2026-08-13 14:27 by flydreamws
[基金申请] 结合人工智能,周易传统文化,filecode打分制来了,3分以上希望很大。 +3 Tide man 2026-08-12 4/200 2026-08-13 08:35 by ZJTJZ
信息提示
请填处理意见