Data Loading...
Appendix C: fit plot - Tools for Science Flipbook PDF
Appendix C: fit& plotExample This appendix provides a longer example using fitand plotto analyze and display data. T
112 Views
77 Downloads
FLIP PDF 102.2KB
Appendix C: fit & plot Example This appendix provides a longer example using fit and plot to analyze and display title=Resonance * border * dpoint Plot is to use exactly the same function as fit: * set f(x)=k1/sqrt(k2^2+(1/(2*pi*k3*x)-2*pi*k4*x)^2) and I can copy&paste2 from the fit screen the values fit found for those parameters: * set K1= 2.262721 K2= 190.7750 K3= 0.9930065E-08 * set K4= 0.8230284E-02 Then I can have plot display that function curve and print a hardcopy: * fcurve * pcopy 2
Warning repeated: Ctrl C kills not copies; highlight and middle-mouse for copy
182
Appendix C: fit & plot Example Resonance .012
RMS Current (A)
.010 .008 .006 .004 .002 10000
15000
20000 25000 Frequency (Hz)
30000
If I go back to the fit screen and type s in answer to the question about matrices: Display covariance/curvature matrices? No, Screen, File [N,S,F] s Results of a fit to the function f(x)=K1/SQRT(K2^2+(1/(2*PI*K3*X)-2*PI*K4*X)^2) at 15:15 on 13-MAY-2009 found a REDUCED chi-square of 3.790897 and a chi-squared of using the following parameter estimates K1= 2.262721 K2= 190.7750 K3= 0.9930065E-08 K4= 0.8230284E-02 and 100 data points at row 2 in the file: LRC.dat
363.9261
Below are displayed the covariance matrix = inverse curvature matrix and the curvature matrix = .5 * Hessian = .5 * matrix of 2nd partials of chi-square see Numerical Recipes 15.5-15.6 for a discussion of these matrices Matrix form:
K1 K2 K3 K4 K1 K2 K3 K4
COVARIANCE -.55E+01 -.46E+03 0.24E-07 -.20E-01
MATRIX: -.46E+03 -.39E+05 0.20E-05 -.17E+01
0.24E-07 0.20E-05 -.11E-15 0.87E-10
-.20E-01 -.17E+01 0.87E-10 -.72E-04
CURVATURE MATRIX: 0.37E+06 -.12E+04 -.12E+04 0.83E+01 0.34E+14 -.54E+11 -.34E+08 0.63E+05
0.34E+14 -.54E+11 0.53E+23 0.56E+17
-.34E+08 0.63E+05 0.56E+17 0.76E+11
Appendix C: fit & plot Example
183
The error results are a total disaster! The usual rule is that the error in a parameter is given by the square root of the diagonal element of the covariance matrix. Thus (supposedly): p δk1 = −.55 × 101 (C.2) p 5 −.39 × 10 (C.3) δk2 = p −.11 × 10−15 (C.4) δk3 = p δk4 = −.72 × 10−4 (C.5) Clearly this is some sort of nonsense as the errors cannot be imaginary. Further note that in every case |δki | ∼ ki . A hint about what’s going on is gained by making a different starting guess for the values of the parameters: * set k1=240 k2=.47 k3=1e-7 k4=.00086 Fit results show the same reduced χ2 , but with radically different parameters: REDUCED chi-squared= K1= 0.2419768 K4= 0.8801518E-03
3.790896 chi-squared= 363.9260 K2= 20.40159 K3= 0.9285606E-07
There is a curve of (L, R, C, V ) values that produces exactly the same resonance curve. Substituting L → αL
(C.6)
R → αR
(C.8)
C → C/α
(C.7)
V
(C.9)
→ αV
(C.10)
into Eq. C.1 results in an unchanged function (the α cancels out). We can rewrite our function with just three parameters: I=r
1 + Q2
I0
f0 f
−
f f0
2
(C.11)
where: I0 = V /R p L/C Q = R 1 √ f0 = 2π LC
(C.12) (C.13) (C.14) (C.15)
184
Appendix C: fit & plot Example
This reparametrization of the function is both required and easier, in the sense that initial guesses for the parameters can be immediately obtained from the raw data: I0 ∼ .012 A is the current at resonance, f0 ∼ 17600 Hz is the resonant frequency, and the√quality factor Q is f0 /∆f , where ∆f ∼ 19500 − 15800 = 3700 is the full width at I = I0 / 2 ∼ .0085 A. * set f(x)=k1/sqrt(1+k2^2*(k3/x-x/k3)^2) * set k1=.012 k2=4.7 k3=17.6e3 * fit Enter list of Ks to vary, e.g. K1-K3,K5 k1-k3 FIT finished with change in Ks implying 4 significant figures 2 iterations used REDUCED chi-squared= 3.751816 chi-squared= 363.9261 K1= 0.1186070E-01 K2= 4.772110 K3= 17605.02 Display covariance/curvature matrices? No, Screen, File [N,S,F] s Results of a fit to the function f(x)=K1/SQRT(1+K2^2*(K3/X-X/K3)^2) at 16:05 on 13-MAY-2009 found a REDUCED chi-square of 3.751816 using the following parameter estimates K1= 0.1186070E-01 K2= 4.772110 and 100 data points at row 2 in the file: LRC.dat
and a chi-squared of K3=
363.9261
17605.02
Below are displayed the covariance matrix = inverse curvature matrix and the curvature matrix = .5 * Hessian = .5 * matrix of 2nd partials of chi-square see Numerical Recipes 15.5-15.6 for a discussion of these matrices Matrix form:
K1 K2 K3 K1 K2 K3
COVARIANCE 0.52E-09 0.25E-06 -.68E-06
MATRIX: 0.25E-06 0.14E-03 -.77E-03
-.68E-06 -.77E-03 0.16E+02
CURVATURE MATRIX: 0.13E+11 -.24E+08 -.24E+08 0.52E+05 -.60E+03 0.15E+01
-.60E+03 0.15E+01 0.63E-01
We now have reasonable results for errors: p δk1 = .52 × 10−9 = 2.3 × 10−5 A p δk2 = .14 × 10−3 = .012 p .16 × 102 = 4 Hz δk3 =
(C.16) (C.17) (C.18)
185
Appendix C: fit & plot Example Another way of estimating errors is to bootstrap: * boots The results: MEANS 1.187386E-02 4.77916 STANDARD DEVIATIONS 2.891049E-05 1.747892E-02
17605.1 8.34518
are similar to those above. With the largish reduced χ2 one might want to enlarge y-error values until reduced χ2 ≈ 1. . . fit calls that a fudge. Note that once the y-errors have been modified the only way to return to the unmodified data is to re-read the file. * fudge * fit Enter list of Ks to vary, e.g. K1-K3,K5 k1-k3 FIT finished with change in chi-square= 9.8335266E-02 1 iterations used REDUCED chi-squared= 0.9858394 chi-squared= 95.62643 K1= 0.1186207E-01 K2= 4.772794 K3= 17605.02 Display covariance/curvature matrices? No, Screen, File [N,S,F] s Results of a fit to the function f(x)=K1/SQRT(1+K2^2*(K3/X-X/K3)^2) at 16:15 on 13-MAY-2009 found a REDUCED chi-square of 0.9858394 using the following parameter estimates K1= 0.1186207E-01 K2= 4.772794 and 100 data points at row 2 in the file: LRC.dat
and a chi-squared of K3=
95.62643
17605.02
Below are displayed the covariance matrix = inverse curvature matrix and the curvature matrix = .5 * Hessian = .5 * matrix of 2nd partials of chi-square see Numerical Recipes 15.5-15.6 for a discussion of these matrices Matrix form:
K1 K2 K3 K1 K2 K3
COVARIANCE 0.20E-08 0.94E-06 -.26E-05
MATRIX: 0.94E-06 0.51E-03 -.29E-02
-.26E-05 -.29E-02 0.60E+02
186
Appendix C: fit & plot Example
Once the error bars have been expanded (fudge) this new fit gives a new covariance matrix from which fudged errors can be determined: p .20 × 10−8 = 4.5 × 10−5 A (C.19) δk1 = p δk2 = .51 × 10−3 = .023 (C.20) p .60 × 102 = 7.7 Hz (C.21) δk3 =
What is the source of the uncomfortably large reduced χ2 ? If you look at the plot you’ll notice that there are glitches at the two frequencies where I = .003 A — that is where the ammeter has automatically switched scales. A scale switch can result in different systematic error and it certainly results in a different voltage burden3 . If we just fit to the data with I > .003 A we can avoid those effects. This data lies between rows 22 and 88 in the file LRC.dat; we can selectively have fit read this data: * set row=22 npoint=67 * read * print Note that when you read all 21 rows of unwanted stuff in the file is echoed to your screen. The command print displays the data that is inside the program, allowing eyeballs to confirm that fit’s current data has I > .003 A. The fit result is: REDUCED chi-squared= K1= 0.1174801E-01
0.6370068 chi-squared= 40.76844 K2= 4.633543 K3= 17598.74
One can achieve almost the same thing by telling fit to limit its analysis to a subset block of the data it holds. When the data was originally read in, it was stored in adjacent cells labeled 1,2,3,. . . ,npoint. We can tell fit to analyze any subset block of that data by changing which cell holds the “first” data point (ibegin) and the length of the run of wanted data. * set ibegin=21 npoint=67 In this case we should not re-read the data, as the intent is to use the data already held in fit. (Earlier in this example changed the y-errors when we commanded fudge, so here we really do need to re-read the unmodified data rather than just use a subset of the modified data.) We can use plot to display this new fitted curve with the entire dataset. All we need to do is change the function, set the newly found parameters k1 -k3 via copy&paste from fit, and then redraw the plot: * clear * set f(x)=k1/sqrt(1+k2^2*(k3/x-x/k3)^2) * set K1= 0.1174801E-01 K2= 4.633543 3
see page 83
K3=
17598.74
187
Appendix C: fit & plot Example * border * dpoint * fcurve Resonance .012
RMS Current (A)
.010 .008 .006 .004 .002 10000
15000
20000 25000 Frequency (Hz)
30000
Mathematica The program Mathematica can also find these nonlinear fits; of course it uses a different set of commands and data structures. You can follow along with this example if you grab a copy of the (slightly modified4 ) data file. At the linux % prompt type: % cp /usr/local/physics/help/LRCm.dat . and a copy of the datafile should appear in your current directory (which is called “.”). There are two version of Mathematica a command line version (math) and a browser version5 (mathematica). The commands are the same in the two versions; the main difference is an absence of user friendly features—like arrow-key line editing, on-line help, drop-down menus—in the kernel version (math). I every much prefer the no-frills math version, but it is an acquired taste. In a terminal, cd to the directory that contains the file LRCm.dat. Start Mathematica by typing to the linux % prompt: % mathematica or, if you’re hard-core: % math 4 5
the first row of text has been removed i.e., with a Graphical User Interface—GUI, pronounced “gooey”—with menus and line editing
188
Appendix C: fit & plot Example
The file LRC.com that you previously copied also contains the Mathematica commands discussed below. Note that in mathematica to execute a command you need to simultaneously hit Shift and Enter, whereas in math it’s just the usual Enter. In Mathematica data is stored in set notation: xy={{x1 , y1 }, {x2 , y2 }, . . . , {xN , yN }} with the errors stored in a different set (which must of course be in the same order): ey={δy1 , δy2 , . . . , δyN } To read the data into Mathematica in this format we need to do some rearrangement: xyey=ReadList["LRCm.dat",{{Number,Number},Number}] xy=Part[xyey,All,1] ey=Part[xyey,All,2] The following line must be typed in exactly (including case) as shown: nlm = NonlinearModelFit[xy, k1/Sqrt[1+k2^2*(k3/x-x/k3)^2],{{k1,.012},{k2,4.7}, {k3,17600}},{x},Weights->1/ey^2,VarianceEstimatorFunction->(1 &)] I hope the parts make some sense to you: notice particularly how the initial parameter guesses have been entered (e.g., {k1,.012}). The result should look at bit like: FittedModel[ √
.0118607 ] 1 + 22.7729 ≪ 1 ≫2
The symbol nlm now contains all the information about the fit; proper requests can extract that information: In[5]:= nlm["ParameterConfidenceIntervalTable"] Out[5]= Estimate Standard Error Confidence Interval 0.0118153 k1 0.0118607 0.000022857 0.011906
k2
k3
4.77209
17605.
0.0116262
4.74902 4.79517
3.98307
17597.1 17612.9
In[6]:= nlm["CovarianceMatrix"] -10 -7 -7 -7 Out[6]= {{5.2244 10 , 2.46067 10 , -7.06763 10 }, {2.46067 10 , 0.000135169,
>
-0.00078076}, {-7.06763 10
In[7]:= nlm["ANOVATable"] Out[7]=
DF
-7 , -0.00078076, 15.8649}}
SS
MS 6
Appendix C: fit & plot Example Model
3
1.88937 10
629790.
Error
97
363.926
3.75182
Uncorrected Total
100
6 1.88974 10
Corrected Total
99
377421.
189
ParameterConfidenceIntervalTable gives us the fitted parameters with usual errors. As with fit, these estimates are the square root of the diagonal elements of the CovarianceMatrix. The quality of the fit is available in the ANOVATable: the Error row, the SS column is the χ2 and the MS column is the reduced χ2 . Mathematica needs some help to do plots with error bars. Additionally the error bars and the data need to be put together in exactly the right way: Needs["ErrorBarPlots‘"] bar=Map[ErrorBar,ey] xyEB=Transpose[{xy,bar}] ErrorListPlot[xyEB] Show[ErrorListPlot[xyEB], Plot[nlm[x], {x, 10000, 30000}], Frame->True]
190
Appendix C: fit & plot Example