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