111 lines
3.2 KiB
FortranFixed
111 lines
3.2 KiB
FortranFixed
subroutine ftest(data1,n1,data2,n2,f,prob)
|
|
implicit none
|
|
integer n1,n2
|
|
double precision f,prob,data1(n1),data2(n2)
|
|
!given the arrays data1 and data2, this routine returns the value of f
|
|
!and the significance as prob. Small values of prob indicate that the
|
|
!two arrays have significantly different variances
|
|
double precision ave1,ave2,df1,df2,var1,var2,betai,
|
|
&xmin,xmax
|
|
call stdmaxmeanmin(n1,data1,var1,ave1,xmin,xmax)
|
|
var1=var1*var1
|
|
call stdmaxmeanmin(n2,data2,var2,ave2,xmin,xmax)
|
|
var2=var2*var2
|
|
if(var1.gt.var2)then
|
|
f=var1/var2
|
|
df1=dble(n1-1)
|
|
df2=dble(n2-1)
|
|
else
|
|
f=var2/var1
|
|
df1=dble(n2-1)
|
|
df2=dble(n1-1)
|
|
endif
|
|
prob=2.0d0*betai(0.5d0*df2,0.5d0*df1,df2/(df2+df1*f))
|
|
if(prob.gt.1.0d0)prob=2.0d0-prob
|
|
return
|
|
end
|
|
|
|
c####################################################################
|
|
subroutine ftestrsq_rms(ymeas,ypred,n0,nparams,rsq,rms,agrind,
|
|
&rmse_norm,rmse_perc,aic,aicc,f,prob)
|
|
implicit double precision (a-h,l,o-z)
|
|
dimension ymeas(n0),ypred(n0),y1(n0),y2(n0)
|
|
fn9999=-9999.0d0
|
|
tiny=1.0d-7
|
|
n=0
|
|
do i=1,n0
|
|
if(dabs(ymeas(i)-fn9999).gt.tiny.and.
|
|
&dabs(ypred(i)-fn9999).gt.tiny)then
|
|
n=n+1
|
|
y1(n)=ymeas(i)
|
|
y2(n)=ypred(i)
|
|
endif
|
|
enddo
|
|
ymin=y1(1)
|
|
ymax=y1(1)
|
|
do i=2,n
|
|
if(y1(i).lt.ymin)ymin=y1(i)
|
|
if(y1(i).gt.ymax)ymax=y1(i)
|
|
enddo
|
|
sum=0.0d0
|
|
rmse_perc=0.0d0
|
|
do 10 i=1,n
|
|
sum=sum+(y1(i)-y2(i))*(y1(i)-y2(i))
|
|
rmse_perc=rmse_perc+(y1(i)-y2(i))*(y1(i)-y2(i))/(y2(i)*y2(i))
|
|
10 continue
|
|
rms=dsqrt(sum/dble(n))
|
|
if(nparams.gt.0)then
|
|
aic=dble(n)*dlog(rms*rms)+2.0d0*dble(nparams)
|
|
aicc=aic+2.0d0*dble(nparams*(nparams+1))/dble(n-nparams-1)
|
|
else
|
|
aic=-9999.0d0
|
|
aicc=-9999.0d0
|
|
endif
|
|
rmse_norm=rms/(ymax-ymin)
|
|
rmse_perc=100.0d0*dsqrt(rmse_perc/dble(n))
|
|
ymean1=0.0d0
|
|
ymean2=0.0d0
|
|
do 20 i=1,n
|
|
ymean1=ymean1+y1(i)
|
|
ymean2=ymean2+y2(i)
|
|
20 continue
|
|
ymean1=ymean1/dble(n)
|
|
ymean2=ymean2/dble(n)
|
|
sum1=0.0d0
|
|
sum2=0.0d0
|
|
sum3=0.0d0
|
|
sum4=0.0d0
|
|
sum5=0.0d0
|
|
do 30 i=1,n
|
|
sum1=(y1(i)-ymean1)*(y2(i)-ymean2)+sum1
|
|
sum2=(y1(i)-ymean1)*(y1(i)-ymean1)+sum2
|
|
sum3=(y2(i)-ymean2)*(y2(i)-ymean2)+sum3
|
|
sum4=(y1(i)-y2(i))*(y1(i)-y2(i))+sum4
|
|
sum5=(dabs(y2(i)-ymean1)+dabs(y1(i)-ymean1))*
|
|
&(dabs(y2(i)-ymean1)+dabs(y1(i)-ymean1))+sum5
|
|
sum6=sum6+(y2(i)-ymean1)*(y2(i)-ymean1)
|
|
30 continue
|
|
if((sum2*sum3).eq.0.0d0)then
|
|
rsq=-9999.0d0
|
|
else
|
|
rsq=sum1/dsqrt(sum2*sum3)
|
|
rsq=rsq*rsq
|
|
endif
|
|
if(sum5.eq.0.0d0)then
|
|
agrind=-9999.0d0
|
|
else
|
|
agrind=1.0d0-sum4/sum5
|
|
endif
|
|
if(rsq.gt.0.0d0)then
|
|
df1=dble(nparams-1)
|
|
df2=dble(n-nparams)
|
|
f=(sum6/df1)/(sum4/df2)
|
|
prob=betai(0.5d0*df2,0.5d0*df1,df2/(df2+df1*f))
|
|
else
|
|
f=-9999.0d0
|
|
prob=-9999.0d0
|
|
endif
|
|
return
|
|
end
|
|
!####################################################################
|