gnuplot demo script: fit.dem

#
# $Id: fit.dem,v 1.5 2004/09/25 03:39:20 sfeam Exp $
#

print "Some examples how data fitting using nonlinear least squares fit"
print "can be done."
print ""

Click here for minimal script to generate this plot

set title 'data for first fit demo'
plot 'lcdemo.dat'
set xlabel "Temperature T  [deg Cels.]"
set ylabel "Density [g/cm3]"
set key below

print "now fitting a straight line to the data :-)"
print "only as a demo without physical meaning"
load 'line.fnc'
y0 = 0.0
m = 0.0
show variables

Click here for minimal script to generate this plot

set title 'all fit params set to 0'
plot 'lcdemo.dat', l(x)

Click here for minimal script to generate this plot

fit l(x) 'lcdemo.dat' via y0, m

Click here for minimal script to generate this plot

set title 'unweighted fit'
plot 'lcdemo.dat', l(x)


Click here for minimal script to generate this plot

fit l(x) 'lcdemo.dat' using 1:2:3 via y0, m

Click here for minimal script to generate this plot

set title 'fit weighted towards low temperatures'
plot 'lcdemo.dat', l(x)


Click here for minimal script to generate this plot

fit l(x) 'lcdemo.dat' using 1:2:4 via y0, m

Click here for minimal script to generate this plot

set title 'bias to high-temperates'
plot 'lcdemo.dat', l(x)

print "now use real single-measurement errors to reach such a result (-> return)"
print "(look at the file lcdemo.dat and compare the columns to see the difference)"

Click here for minimal script to generate this plot

set title 'data with experimental errors'
plot 'lcdemo.dat' using 1:2:5 with errorbars
fit l(x) 'lcdemo.dat' using 1:2:5 via y0, m

Click here for minimal script to generate this plot

set title 'fit weighted by experimental errors'
plot 'lcdemo.dat' using 1:2:5 with errorbars, l(x)

print "It's time now to try a more realistic model function"
load 'density.fnc'
show functions
print "density(x) is a function which shall fit the whole temperature"
print "range using a ?: expression. It contains 6 model parameters which"
print "will all be varied. Now take the start parameters out of the"

Click here for minimal script to generate this plot

load 'start.par'
set title 'initial parameters for realistic model function'
plot 'lcdemo.dat', density(x)
fit density(x) 'lcdemo.dat' via 'start.par'

Click here for minimal script to generate this plot

set title 'fitted to realistic model function'
plot 'lcdemo.dat', density(x)

print  "looks already rather nice? We will do now the following: set"
print  "the epsilon limit higher so that we need more iteration steps"
print  "to convergence. During fitting please hit ctrl-C. You will be asked"
print  "Stop, Continue, Execute: Try everything. You may define a script"
print  "using the FIT_SCRIPT environment variable. An example would be"
print  "'FIT_SCRIPT=plot nonsense.dat'. Normally you don't need to set"
print  "FIT_SCRIPT since it defaults to 'replot'. Please note that FIT_SCRIPT"
print  "cannot be set from inside gnuplot."
print  ""

Click here for minimal script to generate this plot

FIT_LIMIT = 1e-10
fit density(x) 'lcdemo.dat' via 'start.par'

Click here for minimal script to generate this plot

set title 'fit with more iterations'
plot 'lcdemo.dat', density(x)

FIT_LIMIT = 1e-5
print "\nNow a brief demonstration of 3d fitting."
print "hemisphr.dat contains random points on a hemisphere of"
print "radius 1, but we let fit figure this out for us."
print "It takes many iterations, so we limit FIT_MAXITER to 50."
#HBB: made this a lot harder: also fit the center of the sphere
#h(x,y) = sqrt(r*r - (x-x0)**2 - (y-y0)**2) + z0
#HBB 970522: distort the function, so it won't fit exactly:
h(x,y) = sqrt(r*r - (abs(x-x0))**2.2 - (abs(y-y0))**1.8) + z0
x0 = 0.1
y0 = 0.2
z0 = 0.3
r=0.5
FIT_MAXITER=50
set title 'the scattered points, and the initial parameter'
splot 'hemisphr.dat' using 1:2:3, h(x,y)

Click here for minimal script to generate this plot

# we *must* provide 4 columns for a 3d fit. We fake errors=1
fit h(x,y) 'hemisphr.dat' using 1:2:3:(1) via r, x0, y0, z0
set title 'the scattered points, fitted curve'
splot 'hemisphr.dat' using 1:2:3, h(x,y)
print "\n\nNotice, however, that this would converge much faster when"
print "fitted in a more appropriate co-ordinate system:"
print "fit r 'hemisphr.dat' using 0:($1*$1+$2*$2+$3*$3) via r"
print "where we are fitting f(x)=r to the radii calculated as the data"
print "is read from the file. No x value is required in this case."

Click here for minimal script to generate this plot

FIT_MAXITER=0   # no limit : we cannot delete the variable once set

print "\n\nNow an example how to fit multi-branch functions\n"
print  "The model consists of two branches, the first describing longitudinal"
print  "sound velocity as function of propagation direction (upper data),"
print  "the second describing transverse sound velocity (lower data).\n"
print  "The model uses these data in order to fit elastic stiffnesses"
print  "which occur differently in both branches.\n"

Click here for minimal script to generate this plot

load 'hexa.fnc'
load 'sound.par'
set title 'sound data, and model with initial parameters'
plot 'soundvel.dat', vlong(x), vtrans(x)
# Must provide an error estimate for a 3d fit. Use constant 1
fit f(x,y) 'soundvel.dat' using 1:-2:2:(1) via 'sound.par'
#create soundfit.par, reading from sound.par and updating values
update 'sound.par' 'soundfit.par'
print  ""

Click here for minimal script to generate this plot

set title 'pseudo-3d multi-branch fit to velocity data'
plot 'soundvel.dat', vlong(x), vtrans(x)
print  "Look at the file 'hexa.fnc' to see how the branches are realized"
print  "using the data index as a pseudo-3d fit"
print  ""
print  "Next we only use every fifth data point for fitting by using the"
print  "'every' keyword. Look at the fitting-speed increase and at"
print  "fitting result."
print  ""

Click here for minimal script to generate this plot

load 'sound.par'
fit f(x,y) 'soundvel.dat' every 5 using 1:-2:2:(1) via 'sound.par'
set title 'fitted only every 5th data point'
plot 'soundvel.dat', vlong(x), vtrans(x)
print  "When you compare the results (see 'fit.log') you remark that"
print  "the uncertainties in the fitted constants have become larger,"
print  "the quality of the plot is only slightly affected."
print  ""
print  "By marking some parameters as '# FIXED' in the parameter file"
print  "you fit only the others (c44 and c13 fixed here)."
print  ""

Click here for minimal script to generate this plot

load 'sound2.par'
set title 'initial parameters'
plot 'soundvel.dat', vlong(x), vtrans(x)
fit f(x,y) 'soundvel.dat' using 1:-2:2:(1) via 'sound2.par'
set title 'fit with c44 and c13 fixed'
plot 'soundvel.dat', vlong(x), vtrans(x)
print  "This has the same effect as specifying only the real free"
print  "parameters by the 'via' syntax."
print  ""
print  "fit f(x) 'soundvel.dat' via c33, c11, phi0"
print  ""

Click here for minimal script to generate this plot

load 'sound.par'
set title 'initial parameters'
plot 'soundvel.dat', vlong(x), vtrans(x)
fit f(x,y) 'soundvel.dat' using 1:-2:2:(1) via c33, c11, phi0
set title 'fit via c33,c11,phi0'
plot 'soundvel.dat', vlong(x), vtrans(x)

print  "Here comes an example of a very complex function..."
print  ""

Click here for minimal script to generate this plot

set xlabel "Delta [degrees]"
set ylabel "Reflectivity"
set title 'raw data'
#HBB 970522: here and below, use the error column present in moli3.dat:
plot 'moli3.dat' w e

print "now fitting the model function to the data"
load 'reflect.fnc'

#HBB 970522: Changed initial values to something sensible, i.e. 
#  something an experienced user of fit would actually use.
#  FIT_LIMIT is also raised, to ensure a better fit.
eta = 1.2e-4
tc = 1.8e-3
FIT_LIMIT=1e-10

show variables
show functions

Click here for minimal script to generate this plot

set title 'initial parameters'
plot 'moli3.dat' w e, R(x)

Click here for minimal script to generate this plot

fit R(x) 'moli3.dat' u 1:2:3 via eta, tc

Click here for minimal script to generate this plot

set title 'fitted parameters'
replot

#HBB 970522: added comment on result of last fit.
print "Looking at the plot of the resulting fit curve, you can see"
print "that this function doesn't really fit this set of data points."
print "This would normally be a reason to check for measurement problems"
print "not yet accounted for, and maybe even re-think the theoretic"
print "prediction in use."
print ""

print  "You can have a look at all previous fit results by looking into"
print  "the file 'fit.log' or whatever you defined the env-variable 'FIT_LOGFILE'."
print  "Remember that this file will always be appended, so remove it"
print  "from time to time!"
print  ""

Click here for minimal script to generate this plot

reset