2 # $Id: fit.dem,v 1.5 2004/09/25 03:39:20 sfeam Exp $
5 print "Some examples how data fitting using nonlinear least squares fit"
8 pause -1 "first plotting the pure data set (-> return)"
10 set title 'data for first fit demo'
12 set xlabel "Temperature T [deg Cels.]"
13 set ylabel "Density [g/cm3]"
16 print "now fitting a straight line to the data :-)"
17 print "only as a demo without physical meaning"
22 pause -1 "first a plot with all parameters set to zero (-> return)"
23 set title 'all fit params set to 0'
24 plot 'lcdemo.dat', l(x)
25 pause -1 "now start fitting... (-> return)"
26 fit l(x) 'lcdemo.dat' via y0, m
27 pause -1 "now look at the result (-> return)"
28 set title 'unweighted fit'
29 plot 'lcdemo.dat', l(x)
31 pause -1 "see the influence of weights for single data points (-> return)"
32 fit l(x) 'lcdemo.dat' using 1:2:3 via y0, m
33 pause -1 "now look at the result (-> return)"
34 set title 'fit weighted towards low temperatures'
35 plot 'lcdemo.dat', l(x)
37 pause -1 "now prefer the high temperature data (-> return)"
38 fit l(x) 'lcdemo.dat' using 1:2:4 via y0, m
39 pause -1 "now look at the result (-> return)"
40 set title 'bias to high-temperates'
41 plot 'lcdemo.dat', l(x)
43 print "now use real single-measurement errors to reach such a result (-> return)"
44 print "(look at the file lcdemo.dat and compare the columns to see the difference)"
45 pause -1 "(-> return)"
46 set title 'data with experimental errors'
47 plot 'lcdemo.dat' using 1:2:5 with errorbars
48 fit l(x) 'lcdemo.dat' using 1:2:5 via y0, m
49 pause -1 "now look at the result (-> return)"
50 set title 'fit weighted by experimental errors'
51 plot 'lcdemo.dat' using 1:2:5 with errorbars, l(x)
53 print "It's time now to try a more realistic model function"
56 print "density(x) is a function which shall fit the whole temperature"
57 print "range using a ?: expression. It contains 6 model parameters which"
58 print "will all be varied. Now take the start parameters out of the"
59 pause -1 "file 'start.par' and plot the function (-> return)"
61 set title 'initial parameters for realistic model function'
62 plot 'lcdemo.dat', density(x)
63 fit density(x) 'lcdemo.dat' via 'start.par'
64 pause -1 "now look at the result (-> return)"
65 set title 'fitted to realistic model function'
66 plot 'lcdemo.dat', density(x)
68 print "looks already rather nice? We will do now the following: set"
69 print "the epsilon limit higher so that we need more iteration steps"
70 print "to convergence. During fitting please hit ctrl-C. You will be asked"
71 print "Stop, Continue, Execute: Try everything. You may define a script"
72 print "using the FIT_SCRIPT environment variable. An example would be"
73 print "'FIT_SCRIPT=plot nonsense.dat'. Normally you don't need to set"
74 print "FIT_SCRIPT since it defaults to 'replot'. Please note that FIT_SCRIPT"
75 print "cannot be set from inside gnuplot."
77 pause -1 "(-> return)"
79 fit density(x) 'lcdemo.dat' via 'start.par'
80 pause -1 "now look at the result (-> return)"
81 set title 'fit with more iterations'
82 plot 'lcdemo.dat', density(x)
85 print "\nNow a brief demonstration of 3d fitting."
86 print "hemisphr.dat contains random points on a hemisphere of"
87 print "radius 1, but we let fit figure this out for us."
88 print "It takes many iterations, so we limit FIT_MAXITER to 50."
89 #HBB: made this a lot harder: also fit the center of the sphere
90 #h(x,y) = sqrt(r*r - (x-x0)**2 - (y-y0)**2) + z0
91 #HBB 970522: distort the function, so it won't fit exactly:
92 h(x,y) = sqrt(r*r - (abs(x-x0))**2.2 - (abs(y-y0))**1.8) + z0
98 set title 'the scattered points, and the initial parameter'
99 splot 'hemisphr.dat' using 1:2:3, h(x,y)
100 pause -1 "(-> return)"
102 # we *must* provide 4 columns for a 3d fit. We fake errors=1
103 fit h(x,y) 'hemisphr.dat' using 1:2:3:(1) via r, x0, y0, z0
104 set title 'the scattered points, fitted curve'
105 splot 'hemisphr.dat' using 1:2:3, h(x,y)
106 print "\n\nNotice, however, that this would converge much faster when"
107 print "fitted in a more appropriate co-ordinate system:"
108 print "fit r 'hemisphr.dat' using 0:($1*$1+$2*$2+$3*$3) via r"
109 print "where we are fitting f(x)=r to the radii calculated as the data"
110 print "is read from the file. No x value is required in this case."
111 pause -1 "(This is left as an excercise for the user). (-> return)"
112 FIT_MAXITER=0 # no limit : we cannot delete the variable once set
114 print "\n\nNow an example how to fit multi-branch functions\n"
115 print "The model consists of two branches, the first describing longitudinal"
116 print "sound velocity as function of propagation direction (upper data),"
117 print "the second describing transverse sound velocity (lower data).\n"
118 print "The model uses these data in order to fit elastic stiffnesses"
119 print "which occur differently in both branches.\n"
120 pause -1 "(-> return)"
123 set title 'sound data, and model with initial parameters'
124 plot 'soundvel.dat', vlong(x), vtrans(x)
125 # Must provide an error estimate for a 3d fit. Use constant 1
126 fit f(x,y) 'soundvel.dat' using 1:-2:2:(1) via 'sound.par'
127 #create soundfit.par, reading from sound.par and updating values
128 update 'sound.par' 'soundfit.par'
130 pause -1 "(-> return)"
131 set title 'pseudo-3d multi-branch fit to velocity data'
132 plot 'soundvel.dat', vlong(x), vtrans(x)
133 print "Look at the file 'hexa.fnc' to see how the branches are realized"
134 print "using the data index as a pseudo-3d fit"
136 print "Next we only use every fifth data point for fitting by using the"
137 print "'every' keyword. Look at the fitting-speed increase and at"
138 print "fitting result."
140 pause -1 "(-> return)"
142 fit f(x,y) 'soundvel.dat' every 5 using 1:-2:2:(1) via 'sound.par'
143 set title 'fitted only every 5th data point'
144 plot 'soundvel.dat', vlong(x), vtrans(x)
145 print "When you compare the results (see 'fit.log') you remark that"
146 print "the uncertainties in the fitted constants have become larger,"
147 print "the quality of the plot is only slightly affected."
149 print "By marking some parameters as '# FIXED' in the parameter file"
150 print "you fit only the others (c44 and c13 fixed here)."
152 pause -1 "(-> return)"
154 set title 'initial parameters'
155 plot 'soundvel.dat', vlong(x), vtrans(x)
156 fit f(x,y) 'soundvel.dat' using 1:-2:2:(1) via 'sound2.par'
157 set title 'fit with c44 and c13 fixed'
158 plot 'soundvel.dat', vlong(x), vtrans(x)
159 print "This has the same effect as specifying only the real free"
160 print "parameters by the 'via' syntax."
162 print "fit f(x) 'soundvel.dat' via c33, c11, phi0"
164 pause -1 "(-> return)"
166 set title 'initial parameters'
167 plot 'soundvel.dat', vlong(x), vtrans(x)
168 fit f(x,y) 'soundvel.dat' using 1:-2:2:(1) via c33, c11, phi0
169 set title 'fit via c33,c11,phi0'
170 plot 'soundvel.dat', vlong(x), vtrans(x)
172 print "Here comes an example of a very complex function..."
174 pause -1 "first plotting the pure data set (-> return)"
176 set xlabel "Delta [degrees]"
177 set ylabel "Reflectivity"
179 #HBB 970522: here and below, use the error column present in moli3.dat:
182 print "now fitting the model function to the data"
185 #HBB 970522: Changed initial values to something sensible, i.e.
186 # something an experienced user of fit would actually use.
187 # FIT_LIMIT is also raised, to ensure a better fit.
194 pause -1 "first a plot with all parameters set to initial values (-> return)"
195 set title 'initial parameters'
196 plot 'moli3.dat' w e, R(x)
197 pause -1 "now start fitting... (-> return)"
198 fit R(x) 'moli3.dat' u 1:2:3 via eta, tc
199 pause -1 "Hit return to continue"
200 set title 'fitted parameters'
203 #HBB 970522: added comment on result of last fit.
204 print "Looking at the plot of the resulting fit curve, you can see"
205 print "that this function doesn't really fit this set of data points."
206 print "This would normally be a reason to check for measurement problems"
207 print "not yet accounted for, and maybe even re-think the theoretic"
208 print "prediction in use."
211 print "You can have a look at all previous fit results by looking into"
212 print "the file 'fit.log' or whatever you defined the env-variable 'FIT_LOGFILE'."
213 print "Remember that this file will always be appended, so remove it"
214 print "from time to time!"
216 pause -1 "Done with fitting demo (-> return)"