应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

CyRhmU.jpeg
小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » Fortran » 请问，我想用fortran计算统计中的p值以及95%信度空间
51/1返回列表
查看: 1113  |  回复: 4
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
本帖产生 1 个 程序强帖 ，点击这里进行查看

ian_zhangty

木虫 (著名写手)

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

请问有没有现成的程序，谢谢
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼2011-08-05 16:56:12
已阅   回复此楼   关注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 上找找

其实，这种事情，找个统计软件做了就结了……
赞一下(1人) 回复此楼
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编的，如果有直接的程序就非常省事了
赞一下(1人) 回复此楼
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
赞一下(2人) 回复此楼
5楼2011-08-06 10:41:54
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 瘦竹云淡风轻 的主题更新
51/1返回列表
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭