home *** CD-ROM | disk | FTP | other *** search
Text File | 1992-09-03 | 17.9 KB | 480 lines | [ TEXT/NCII]
################ # Fitting Functions to Arrays # # # program FitLine(xData,yData, A,B, r,ErrorEstimate)#Finds A,B for y=A x + b # # program FitPoly(xData,yData,Sigma, J,A, Prob) #Find A to fit x,yData[1…N] with G(x)=A[1]+A[2]x+…+A[J] x^(J-1). # function Poly(A,x) # Returns A[1] + A[2] x + … + A[J] x^(J-1). # # program FitLegendre(xData,yData,Sigma, J,A, Prob) #Find A[1…J] such that G(x) = A[j] • P_j(x) fits yData[1…N]. # function CalcLegendre(J, x) # Returns P[l,1…size(x)] = P_l(x) for l=1…J. # function LegendreModel(A,x) # ReturnsG(x) = sum over i=1…J { A[i] P_i(x) } # # program FitUser(xData,yData,Sigma, J,A, Prob, BasisUser) # Fits yData to a sum of user supplied basis functions. # program SampleBasisUser(Jp1, zIO) # Sample BasisUser. See program FitUser. # # # This text file explains and gives examples of # some of the routines in the file All Library Routines, which # should be Opened before trying any of these examples. # ################# ############# # Here are the xCOD's : xFIT_FUNCS # xCOD xFIT_FUNCS(prog BasisFuncs(Jp1:num; zIO:array); n:num; x,y,dy:array; j:num; A,Avar, w1,w2,w3:array; chisq,err:num), fits y=sum over j A(j)*Fj(x) to x[1…n], y[1…n]±dy[1…n]. BasisFuncs calculates J=Jp1-1 functions zIO[1…j]=Fj(x) for single x=zIO[Jp1]. Work arrays: w1[1…n,1…j],w2[1…j,1…j],w3[1…j]. Output: A[1…j]=coeffs, Avar[1…j,1…j]=covariance o fA, ChiSquare of fit. Err: -1=> didn't converge; an external program. xBASIS_POLY # xCOD xBASIS_POLY(Jp1:num; zIO:array), calculates J=Jp1-1 polynomial basis functions for a single value of x. Input: x=zIO(jp1). Output: zIO(1…j)=x^j. Used with xFIT_FUNCS; an external program. xBASIS_LEGENDRE # xCOD xBASIS_LEGENDRE(Jp1:num; zIO:array), calculates J=Jp1-1 Legendre basis functions for a single value of x. Input: x=zIO(jP1). Output: zIO(1…j)=P_j(x). Legendre definition: P1(x)=1; P2(x)=x; (n+1) P[n+1](x)=(2n+1) x Pn(x) - n P[n-1](x). Used with xFIT_FUNCS; an external program. ############# # All the routines in this file are based on xFIT_FUNCS, which finds the # coefficients A[1…j] such that G(x)= sum over j of { A[j] Fj(x) } best # fits xData[1…N], yData[1…N], sigma[1…n] = uncertainties in yData, # where Fj(x) are basis functions defined by either another xCOD # or an NCII user program. Examples of both are given below. # # Besides A[1…j], the ChiSqaure (c2) "goodness" of the fit is returned, # which is ChiSq = sum( yData - G(xData))^2, # as well as a co-variance matrix Avar[1…j,1…j] describing # the uncertainties in the fit. The square root of the diagonal elements # of this matrix are the standard deviations (sigmas) of the expected # errors in the coefficients A. # # Note that if we assume the model curve fits the data, then # an estimate of the uncertainties inherent in yData is # SigmaGuess = sqrt( ChiSq/ size(yData) ). # ############### # Fit a line through some data points. # program FitLine(xData,yData, A,B, r,ErrorEstimate)#Finds A,B for y=A x + b . var N,J,jj,nn,Sigma,C,Cvar,w1,w2,w3,ChiSq, err . # Inputs: . # xData[1…n], yData[1…n] . # Outputs: . # A = slope, B = intercept of best fit line y = A x + b . # r = correlation coefficient, -1<=r<=1. . # (r near zero implies that the data is not correlated, i.e. not linear.) . # ErrorEstimate = approximate errors in yData assuming fit to line is OK. . begin . J = 2; # number of coeffients to fine . jj=1…J; . N = size(xData); # number of data points. . nn = 1…N; . if size(yData)<>N then . Print(" ERROR: xData,yData not same size") . else . begin . C[jj] = 0; # space for coefficients . Cvar[jj,jj]=0; # and co-variance array. . w1[nn,jj] = 0; # work space . w2[jj,jj] = 0; . w3[jj] = 0; . ChiSq = 0; err=0; . Sigma[nn] = 1; # assume uniform y uncertainties . xFIT_FUNCS(xBASIS_POLY, N,xData,yData,Sigma, J,C,Cvar, w1,w2,w3, ChiSq, err) . if err<>0 then Print(" ERROR in xFIT_FUNCS = "+err) . A = C(2); B=C(1); # return coefficients . r = -Cvar(1,2)/sqrt( Cvar(1,1)*Cvar(2,2) ) # and correlation coeff r. . ErrorEstimate = sqrt(ChiSq/(N-J)) . end # of else . end . # Fit a general polynomial of degree J-1 to xData, yData, Sigma. # This essentially just calls xFIT_FUNCS with xBASIS_POLY. # The covariance of the A polynomial coefficients isn't used in # this interface routine; modify it if you want this information. # # Prob = the probability that a fit to the data this poor could have # been caused by random chance alone, assuming that the errors # are distributed normally with the given Sigma. A tiny value # of Prob means that this degree polynomial does not fit the data within # the given uncertainties. # program FitPoly(xData,yData,Sigma, J,A, Prob) #Find A to fit x,yData[1…N] with G(x)=A[1]+A[2]x+…+A[J] x^(J-1). . var N,nn,jj,Avar,w1,w2,w3,ChiSq, err . # Inputs: . # xData, yData, Sigma = yData errors, . # J = number of parameters A[j] to fit, . # ( with model G(x) = A[1] + A[2] x + a[3] x^2 + … + a[J] x^(J-1) ) . # Outputs: . # A[1…J] = best fit coefficients, . # Prob = probability that a fit this bad could be normal errors with . # stnd. devs. Sigma. . begin . jj=1…J; . N = size(xData); # number of data points. . nn = 1…N . if size(Sigma)=1 then . Sigma[nn]=Sigma # make sure sigma is an array . if size(Sigma)<>N then . Print("ERROR: xData, Sigma not same size") . else if size(yData)<>N then . Print(" ERROR: xData,yData not same size") . else . begin . A[jj] = 0; . Avar[jj,jj]=0; # and co-variance array. . w1[nn,jj] = 0; # work space . w2[jj,jj] = 0; . w3[jj] = 0; . ChiSq = 0; err=0; . xFIT_FUNCS(xBASIS_POLY, N,xData,yData,Sigma, J,A,Avar, w1,w2,w3, ChiSq, err) . if err<>0 then Print(" ERROR in xFIT_FUNCS = "+err) . Prob = ChiSqProbability( ChiSq, J) # see below. . end # of else . end . # Evaluate the Polynomial given by the A[1…J] coefficients. function Poly(A,x) # Returns A[1] + A[2] x + … + A[J] x^(J-1). . var jj, J, ans . # Inputs: A[1…J] coeffs., . # x = number or array . # Output : . # Poly = A[1] + A[2] x + A[3] x^2 + … + A[J] x^(J-1) . begin . J = size(A) . ans = A[J]; . for jj= 1 … J-1 do ans= A[J-jj] + x*ans . Poly = ans . end . # # ( IncompleteGamma and ChiSqProbability are in the Special Functions file. ) # # Fit a sum of J Legendre polynomials to xData, yData, Sigma. # This essentially just calls xFIT_FUNCS with xBASIS_LEGENDRE. # The covariance of the A polynomial coefficients isn't used in # this interface routine; modify it if you want this information. # # Prob = the probability that a fit to the data this poor could have # been caused by random chance alone, assuming that the errors # are distributed normally with the given Sigma. A tiny value # of Prob means that this degree polynomial does not fit the data within # the given uncertainties. # program FitLegendre(xData,yData,Sigma, J,A, Prob) #Find A[1…J] such that G(x) = A[j] • P_j(x) fits yData[1…N]. . var N,nn,jj,Avar,w1,w2,w3,ChiSq, err . # Inputs: . # xData[1…N], yData[1…N] data points, . # Sigma (real or array) = yData errors, . # J = number of parameters A[j] to fit, . # ( with g(x) = A[1] P_1(x) + A[2] P_2(x) + … + A[J] P_j(x) ) . # Outputs: . # A[1…J] = best fit coefficients, . # Prob = probability of this fit matching data with given Sigma. . # . # Note that here the Legendre polys. are P_j with j=1,2,3,…, . # NOT j=0,1,2,… as is the convention sometimes. . # . begin . jj=1…J; . N = size(xData); # number of data points. . nn = 1…N . if size(Sigma)=1 then . Sigma[nn]=Sigma # make sure sigma is an array . if size(Sigma)<>N then . Print("ERROR: xData, Sigma not same size") . else if size(yData)<>N then . Print(" ERROR: xData,yData not same size") . else . begin . Avar[jj,jj]=0; # and co-variance array. . w1[nn,jj] = 0; # work space . w2[jj,jj] = 0; . w3[jj] = 0; . ChiSq = 0; err=0; . A[jj] = 0; . xFIT_FUNCS(xBASIS_LEGENDRE, N,xData,yData,Sigma, J,A,Avar, w1,w2,w3, ChiSq, err) . if err<>0 then Print(" ERROR in xFIT_FUNCS = "+err) . Prob = ChiSqProbability( ChiSq, J) # see below. . end # of else . end . # Find J Legendre polynomials, P[1…J, 1…size(x)]. # This could also be done with xBASIS_LEGENDRE, but # that would mean looping over x, rather than j, which is probably smaller. # function CalcLegendre(J, x) # Returns P[l,1…size(x)] = P_l(x) for l=1…J. . var Pz, n . # Input: x = array of points to calulate Legendre polynomials on. . # J = positive integer > 0 . # Output: CalcLegendre = P[1…J, 1…size(x)] with . # P[j,] = one of the polynomials = P_j(x) . # NOTE that here the lowest possible J is 1, rather than 0 . # as is sometimes done. Thus P_j is a (j-1)'th order polynomial in x. . begin . Pz[1…J, 1…size(x)] = 0 # reserve space . Pz[1,] = 1; . if J>1 then . begin . Pz[2,] = x . for n= 2…(J-1) do Pz[n+1,] = 1/(n+1) { (2n+1)*x*Pz[n,] - n*Pz[n-1,] } . end . CalcLegendre = Pz . end . # Evaluate the sum of the Legendre polynomials P_l(x) times the A[l] coefficients # on the vector x. This gives the "best fit" curve that FitLegendre calculates. function LegendreModel(A,x) # ReturnsG(x) = sum over i=1…J { A[i] P_i(x) } . var P . # Input: A[j] coefficients . # x = array . # Output: . # LegendreModel = G(x) = sum over i=1…J A[i] P_i(x) . begin . P = CalcLegendre(size(A), x) . LegendreModel = A • P # dot product does sum j A[j]*P[j,] . end . # Fit a linear combination of user supplied functions # to xData, yData, just like the routines above. # # See below or the last example for a sample of BasisUser. # program FitUser(xData,yData,Sigma, J,A, Prob, BasisUser) # Fits yData to a sum of user supplied basis functions. . var N,nn,jj,Avar,w1,w2,w3,ChiSq, err . # Input: . # BasisUser = program to calculate the basis, F_j(x), . # like SampleBasisUser or xBASI_POLY. . # xData[1…N], yData[1…N] = data points . # Sigma (real or array) = uncertainties in yData . # J = number of A[1…J] parameters in basis . # Output: . # A[1…J] best-fit parameters, such that . # G(x) = sum over i=1…J { A[i] F_i(x) } fits yData. . # Prob = probability that the poorness of the fit is due to chance. . begin . jj=1…J; . N = size(xData); # number of data points. . nn = 1…N . if size(Sigma)=1 then . Sigma[nn]=Sigma # make sure sigma is an array . if size(Sigma)<>N then . Print("ERROR: xData, Sigma not same size") . else if size(yData)<>N then . Print(" ERROR: xData,yData not same size") . else . begin . Avar[jj,jj]=0; # and co-variance array. . w1[nn,jj] = 0; # work space . w2[jj,jj] = 0; . w3[jj] = 0; . ChiSq = 0; err=0; . A[jj] = 0; . xFIT_FUNCS(BasisUser, N,xData,yData,Sigma, J,A,Avar, w1,w2,w3, ChiSq, err) . if err<>0 then Print("ERROR in xFIT_FUNCS = "+err) . Prob = ChiSqProbability( ChiSq, J) # see below. . end # of else . end . # # This SampleBasisUser is here only to show the format # of a user routine that FitUser needs. # program SampleBasisUser(Jp1, zIO) # Sample BasisUser. See program FitUser. . var x, y, J . # Input: . # Jp1 = J+1 = size of zIO array . # zIO[Jp1] = x . # Output: . # zIO[1…J] = F_j(x) . # (This one a demo which does F_j(x) = x^(j-1), . # exactly the same as xCOD xBASIS_POLY does, only much slower.) . begin . x = zIO(Jp1) . J = Jp1 - 1 . zIO[1] = 1 . for jj=2…J do zIO[jj] = x*zIO[jj-1] . end . ####### # Examples: ###### ############## # Fitting a line to some data points: # xData=[1,2,3,4]; yData=[1.1, 1.9, 3.2, 3.9]; # Here are some points. # FitLine(xData, yData, a,b,r, ApproxError) # calulcate the best fit line. # This returns: a a = 0.97 # slope b b = 0.1 # intercept. r r = 0.91287093 # correlation coeff: yes, this data is almost a perfect line. ApproxError ApproxError = 0.17748239 # approximate errors in y. # yFit = a xData + b yFit[1…4] = [1.07, 2.04, 3.01, 3.98] # The line y(xData). # table(xData, yData, yFit) # xData yData yFit # - - - - - - - - - - - - - - - - - - - 1 1.1 1.07 2 1.9 2.04 3 3.2 3.01 4 3.9 3.98 # # To see a plot of this, plot : # [x,y]: [xData, yData] # as a point plot , and ( Item 1 ) # or # [xData, yData, 0, 2 ApproxError] # as points with error bars. # # Range: x = xMin…xMax@3xPix # as a curve. ( Item 2 ) # [x,y]: [x, a x + b] # ############## # A quadratic to these same points, with some fairly small error bars: # sigma = [0.2, 0.04, 0.04, 0.2] ; # errors on yData. # Now fit guess(x) = C(1) + C(2)*x + C(3)*x^2. FitPoly(xData,yData,Sigma, 3,C, prob) prob prob = 0.0110876 # a fit this bad would happen ~ 1% of the time by chance. c C[1…3] = [0.09826, 0.8713, 0.02522] quadFit = c[1] + c[2] xData + c[3] xData^2 quadFit[1…4] = [0.99478, 1.94174, 2.93914, 3.98698] # ( Or use quadFit = Poly(C,xData), which does the same thing.) # # Table(xData, yData, sigma, quadFit) # evaluate this without front "#" # xData , yData ± sigma , quadFit # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1 1.1 0.2 0.99478 2 1.9 0.04 1.94174 3 3.2 0.04 2.93914 4 3.9 0.2 3.98698 # So, with these errors, the data are NOT well fit by a quadratic. # # To see a plot of this, plot : # [x,y]: [xData, yData, 0, 2Sigma] # as a point plot with error bars, and # # Range: x = xMin…xMax@3xPix # as a curve. # [x,y]: [x, Poly(C,x)] # ############## # Some Legendre polynomials, which are meaningful for -1<=x<=1: # xData = -.8….8@.2 xData[1…9] = [-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8] yData = [0, .2,.2,.2, 0, -.2,-.2,-.2,0]; # a semi-sinusoid. Sigma = 0.05; # errors for yData FitLegendre(xData, yData, sigma, 4, Al, prob) # prob prob = 0.51972608 # so this fit is reasonable. # Al Al[1…4] = [ 1.7294e-21, -0.04377, -1.3042e-20, 0.48822] # Notice that A[2] and A[4], which are the coefficients fo the # odd polynomials, are the only significant ones, since yData is # also odd. # yLeg = LegendreModel(Al,xData) yLeg[1…9] = [0.05455, 0.28013, 0.30067, 0.18451, 1.0424e-20, -0.18451, -0.30067, -0.28013, -0.05455] # # table(xData, yData,sigma, yLeg) # xData , yData ± Sigma , yLeg # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.8 0 0.05 0.05455 -0.6 0.2 0.05 0.28013 -0.4 0.2 0.05 0.30067 -0.2 0.2 0.05 0.18451 0 0 0.05 1.0424e-20 0.2 -0.2 0.05 -0.18451 0.4 -0.2 0.05 -0.30067 0.6 -0.2 0.05 -0.28013 0.8 0 0.05 -0.05455 # # To see the Legendre polynomials themselves, define x=-1…1@.05; P = CalcLegendre(4,x); # and plot # [x,y]: [x,P[1,]]; [x,P[2,]] ; [x,P[3,]] ; [x,P[4,]] # from x=-1 to 1, which shows the first 4. # # Note that if dx= x[2] - x[1], then these arrays have # P[i,] • P[j,] * dx = 0 if i<>j (i.e. they are orthogonal). # = 2/(2j+3) if i=j. # # To plot the data and the fit, plot # [x,y]: [xData,yData,0,2Sigma] # as a point plot with error bars, # # and # Range: x=xMin…xMax@3xPix # [x,y]: [x, LegendreModel(Al, x)] # as a curve. # ############## # Finally, an example of calling xFIT_FUNCS with # an NC user program as a fitting function. # # Suppose we wish to fit a function of the form # G(x) = A[1] sin(x) + A[2] exp(x) # to the data points xData = [1,2,3,4]; yData=[1,2,3,4]; Sigma = 1; Jparam=2; # # # Then the routine to evalute G(x) should be program Basis_Strange(jP1, zIO) . var x . begin # n = number of basis functions to evaluate at zIO[n+1]. . x = zIO[jP1] # jP1=J+1 should = 3 here. . zIO[1] = sin(x) . zIO[2] = exp(x) . end . # # Then, just as for the polynomial and Legendre examples above, # Find Auser: FitUser(xData,yData,Sigma, Jparam,Auser, Prob, Basis_Strange) # prob # probability of this fit : Prob = 0.61348783 # This fit is OK (since the errors are large). Auser Auser[1…2] = [1.3807, 0.09774] # the coefficients yFit = Auser[1] sin(xData) + Auser[2] exp(xData) # the fit to the data yFit[1…4] = [1.4275, 1.97767, 2.158, 4.29149] # # Table(xData,yData,sigma, yFit) # xData yData Sigma yFit # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1 1 1 1.4275 2 2 1 1.97767 3 3 1 2.158 4 4 1 4.29149 
