55 lines
1.4 KiB
FortranFixed
55 lines
1.4 KiB
FortranFixed
subroutine linregres(npoints0,x0,y0,Sign_Level,a,b,
|
|
&r,bradius,isslopezero)
|
|
implicit none
|
|
!fit for y=a+bx
|
|
integer npoints0,isslopezero
|
|
!isslopezero=-1, slope differs from zero
|
|
!isslopezero=1, slope does not differ from zero
|
|
double precision x0(npoints0),y0(npoints0),r,a,b,
|
|
&bradius,Sign_Level
|
|
integer i,npoints
|
|
double precision xmean,ymean,lxx,lyy,lxy,seb,t0,
|
|
&tstat,tstatr,student_t,fn9999,tiny,x(npoints0),y(npoints0)
|
|
parameter(fn9999=-9999.0d0,tiny=1.0d-7)
|
|
|
|
npoints=0
|
|
do i=1,npoints0
|
|
if(dabs(x0(i)-fn9999).gt.tiny.and.
|
|
&dabs(y0(i)-fn9999).gt.tiny)then
|
|
npoints=npoints+1
|
|
x(npoints)=x0(i)
|
|
y(npoints)=y0(i)
|
|
endif
|
|
enddo
|
|
|
|
xmean=0.0d0
|
|
ymean=0.0d0
|
|
do i=1,npoints
|
|
xmean=xmean+x(i)
|
|
ymean=ymean+y(i)
|
|
enddo
|
|
xmean=xmean/dble(npoints)
|
|
ymean=ymean/dble(npoints)
|
|
lxx=0.0d0
|
|
lyy=0.0d0
|
|
lxy=0.0d0
|
|
do i=1,npoints
|
|
lxx=lxx+(x(i)-xmean)**2
|
|
lyy=lyy+(y(i)-ymean)**2
|
|
lxy=lxy+(x(i)-xmean)*(y(i)-ymean)
|
|
enddo
|
|
b=lxy/lxx
|
|
a=ymean-b*xmean
|
|
r=lxy/dsqrt(lxx*lyy)
|
|
seb=dsqrt((lyy-b*lxy)/(lxx*dble(npoints-2)))
|
|
t0=student_t(npoints-2,Sign_Level)
|
|
bradius=seb*t0
|
|
tstat=dabs(b/seb)
|
|
if(tstat.gt.t0)then
|
|
isslopezero=-1
|
|
else
|
|
isslopezero=1
|
|
endif
|
|
return
|
|
end
|