diff options
Diffstat (limited to 'pkg/utilities/pfregres.f')
| -rw-r--r-- | pkg/utilities/pfregres.f | 183 |
1 files changed, 183 insertions, 0 deletions
diff --git a/pkg/utilities/pfregres.f b/pkg/utilities/pfregres.f new file mode 100644 index 00000000..2edb1884 --- /dev/null +++ b/pkg/utilities/pfregres.f @@ -0,0 +1,183 @@ +c subroutine regres.f +c +c source +c Bevington, pages 172-175. +c +c purpose +c make a mulitple linear regression fit to data with a specified +c function which is linear in coefficients +c +c usage +c call regres (x, y, sigmay, npts, nterms, m, mode, yfit, +c a0, a, sigma0, sigmaa, r, rmul, chisqr, ftest) +c +c description of parameters +c x - array of points for independent variable +c y - array of points for dependent variable +c sigmay - array of standard deviations for y data points +c npts - number of pairs of data points +c nterms - number of coefficients +c m - array of inclusion/rejection criteria for fctn +c mode - determines method of weighting least-squares fit +c +1 (instrumental) weight(i) = 1./sigmay(i)**2 +c 0 (no weighting) weight(i) = 1. +c -1 (statistical) weight(i) = 1./y(i) +c yfit - array of calculated values of y +c a0 - constant term +c a - array of coefficients +c sigma0 - standard deviation of a0 +c sigmaa - array of standard deviations for coefficients +c r - array of linear correlation coefficients +c rmul - multiple linear correlation coefficient +c chisqr - reduced chi square for fit +c ftest - value of f for test of fit +c +c subroutines and function subprograms required +c fctn (x, i, j, m) +c evaluates the function for the jth term and the ith data point +c using the array m to specify terms in the function +c matinv (array, nterms, det) +c inverts a symmetric two-dimensional matrix of degree nterms +c and calculates its determinant +c +c comments +c (dim npts changed 100->1000 21-may-84 dct) +c dimension statement valid for npts up to 100 and nterms up to 10 +c sigmaag changed to sigmaa in statement following statement 132 +c + subroutine pfregs (x,y,sigmay,npts,nterms,m,mode,yfit, + *a0,a,sigma0,sigmaa,r,rmul,chisqr,ftest,fctn) + double precision array,sum,ymean,sigma,chisq,xmean,sigmax + dimension x(1),y(1),sigmay(1),m(1),yfit(1),a(1),sigmaa(1), + *r(1) + dimension weight(1000),xmean(10),sigmax(10),array(10,10) + external fctn +c +c initialize sums and arrays +c +11 sum=0. + ymean=0. + sigma=0. + chisq=0. + rmul=0. + do 17 i=1,npts +17 yfit(i)=0. +21 do 28 j=1,nterms + xmean(j)=0. + sigmax(j)=0. + r(j)=0. + a(j)=0. + sigmaa(j)=0. + do 28 k=1,nterms +28 array(j,k)=0. +c +c accumulate weighted sums +c +30 do 50 i=1,npts +31 if (mode) 32,37,39 +32 if (y(i)) 35,37,33 +33 weight(i)=1./y(i) + goto 41 +35 weight(i)=1./(-y(i)) + goto 41 +37 weight(i)=1. + goto 41 +39 weight(i)=1./sigmay(i)**2 +41 sum=sum+weight(i) + ymean=ymean+weight(i)*y(i) + do 44 j=1,nterms +44 xmean(j)=xmean(j)+weight(i)*fctn(x,i,j,m) +50 continue +51 ymean=ymean/sum + do 53 j=1,nterms +53 xmean(j)=xmean(j)/sum + fnpts=npts + wmean=sum/fnpts + do 57 i=1,npts +57 weight(i)=weight(i)/wmean + +c +c accumulate matrices r and array +c +61 do 67 i=1,npts + sigma=sigma+weight(i)*(y(i)-ymean)**2 + do 67 j=1,nterms + sigmax(j)=sigmax(j)+weight(i)*(fctn(x,i,j,m)-xmean(j))**2 + r(j)=r(j)+weight(i)*(fctn(x,i,j,m)-xmean(j))*(y(i)-ymean) + do 67 k=1,j +67 array(j,k)=array(j,k)+weight(i)*(fctn(x,i,j,m)-xmean(j))* + *(fctn(x,i,k,m)-xmean(k)) +71 free1=npts-1 +72 sigma=dsqrt(sigma/free1) + do 78 j=1,nterms +74 sigmax(j)=dsqrt(sigmax(j)/free1) + r(j)=r(j)/(free1*sigmax(j)*sigma) + do 78 k=1,j + array(j,k)=array(j,k)/(free1*sigmax(j)*sigmax(k)) +78 array(k,j)=array(j,k) +c +c invert symmetric matrix +c +81 call matinv (array,nterms,det) + if (det) 101,91,101 +91 a0=0. + sigma0=0. + rmul=0. + chisqr=0. + ftest=0. + goto 150 +c +c calculate coefficients, fit, and chi square +c +101 a0=ymean +102 do 108 j=1,nterms + do 104 k=1,nterms +104 a(j)=a(j)+r(k)*array(j,k) +105 a(j)=a(j)*sigma/sigmax(j) +106 a0=a0-a(j)*xmean(j) +107 do 108 i=1,npts +108 yfit(i)=yfit(i)+a(j)*fctn(x,i,j,m) +111 do 113 i=1,npts + yfit(i)=yfit(i)+a0 +113 chisq=chisq+weight(i)*(y(i)-yfit(i))**2 + freen=npts-nterms-1 +115 chisqr=chisq*wmean/freen +c +c calculate uncertainties +c +121 if (mode) 122,124,122 +122 varnce=1./wmean + goto 131 +124 varnce=chisqr +131 do 133 j=1,nterms +132 sigmaa(j)=array(j,j)*varnce/(free1*sigmax(j)**2) + if (sigmaa(j)) 835, 835, 836 +835 sigmaa(j) = 0.0 + goto 133 +836 sigmaa(j)=sqrt(sigmaa(j)) +133 rmul=rmul+a(j)*r(j)*sigmax(j)/sigma + freej=nterms +c +noao: When rmul = 1, the following division (stmt 135) would blow up. +c It has been changed so ftest is set to -99999. in this case. + if (rmul) 935, 136, 136 +935 ftest = -99999. + rmul = -99999. + goto 141 +c -noao +136 if (1.0 - abs(rmul)) 1035, 1036, 1037 +1035 rmul=-99999. + ftest = -99999. + goto 141 +1036 ftest = -99999. + rmul = 1.0 + goto 141 +1037 ftest=(rmul/freej)/((1.-rmul)/freen) + rmul=sqrt(rmul) +141 sigma0=varnce/fnpts + do 145 j=1,nterms + do 145 k=1,nterms +145 sigma0=sigma0+varnce*xmean(j)*xmean(k)*array(j,k)/ + *(free1*sigmax(j)*sigmax(k)) +146 sigma0=sqrt(sigma0) +150 return + end |