48 lines
1.5 KiB
FortranFixed
48 lines
1.5 KiB
FortranFixed
subroutine linuncertainty(nsamp,a1_mean,a1_sigma,b1_mean,
|
|
&b1_sigma,y1_mean,y1_sigma,a2_mean,a2_sigma,b2_mean,b2_sigma,
|
|
&y2_mean,y2_sigma,x1,x2,x1_mean,x1_sigma,x2_mean,x2_sigma)
|
|
implicit none
|
|
!given:
|
|
!y1=a1*x1+b1*x2
|
|
!y2=a2*x1+b1*x2
|
|
!y1 ~ (y1_mean, y1_sigma)
|
|
!y2 ~ (y2_mean, y2_sigma)
|
|
!a1 ~ (a1_mean, a1_sigma)
|
|
!b1 ~ (b1_mean, b1_sigma)
|
|
!a2 ~ (a2_mean, a2_sigma)
|
|
!b2 ~ (b2_mean, b2_sigma)
|
|
!find:
|
|
!x1 ~ (x1_mean, x1_sigma)
|
|
!x2 ~ (x2_mean, x2_sigma)
|
|
integer nsamp,i
|
|
double precision a1_mean,a1_sigma,b1_mean,b1_sigma,y1_mean,
|
|
&y1_sigma,a2_mean,a2_sigma,b2_mean,b2_sigma,y2_mean,y2_sigma,
|
|
&x1(nsamp),x2(nsamp),x1_mean,x1_sigma,x2_mean,x2_sigma,gasdev2,
|
|
&a1,b1,a2,b2,y1,y2
|
|
x1_mean=0.0d0
|
|
x2_mean=0.0d0
|
|
do i=1,nsamp
|
|
y1=gasdev2()*y1_sigma+y1_mean
|
|
y2=gasdev2()*y2_sigma+y2_mean
|
|
a1=gasdev2()*a1_sigma+a1_mean
|
|
b1=gasdev2()*b1_sigma+b1_mean
|
|
a2=gasdev2()*a2_sigma+a2_mean
|
|
b2=gasdev2()*b2_sigma+b2_mean
|
|
x1(i)=(b2*y1-b1*y2)/(a1*b2-a2*b1)
|
|
x2(i)=(a1*y2-a2*y1)/(a1*b2-a2*b1)
|
|
x1_mean=x1_mean+x1(i)
|
|
x2_mean=x2_mean+x2(i)
|
|
enddo
|
|
x1_mean=x1_mean/dble(nsamp)
|
|
x2_mean=x2_mean/dble(nsamp)
|
|
x1_sigma=0.0d0
|
|
x2_sigma=0.0d0
|
|
do i=1,nsamp
|
|
x1_sigma=x1_sigma+(x1(i)-x1_mean)*(x1(i)-x1_mean)
|
|
x2_sigma=x2_sigma+(x2(i)-x2_mean)*(x2(i)-x2_mean)
|
|
enddo
|
|
x1_sigma=x1_sigma/dble(nsamp)
|
|
x2_sigma=x2_sigma/dble(nsamp)
|
|
return
|
|
end
|