[gnuplot] out and input of fit parameters

leonardo2887 · 10-11-2010, 08:18 AM

Hello to everybody,

I'm writing a script to fit long series of data.
This is the part where I don't exactly know how to continue:

Code:

until [ $i -gt 95 ]; do 
let j=i+2
let k=j+2
awk '{printf $1" "$2"\n"}' 'oil15_0'$i'.dat' 'oil15_0'$j'.dat' 'oil15_0'$k'.dat' >oil15_$i$j$k.dat
gnuplot <<EOF
set logscale x
f(x) = a*(p*exp(-b*x)+(1-p)*exp(-c*x))+d
f1(x)=a*exp(-b*x)+c
a=0.8
b=5
c=1
p=1
d=1
unset key
fit [0.0001:1000][0.8:2] f1(x) 'oil15_$i$j$k.dat' via a, b, c
set term svg
set out 'dls0$i.svg'
set pointsize 1
plot [0.0001:1000][0.8:2] f1(x) w l lw 2, 'oil15_$i$j$k.dat' ls 6
set term wxt
replot
EOF
let i=i+6
done

After fitting, gnuplot writes the file fit.log, here an example:

Code:

*******************************************************************************
Sun Oct 10 23:35:22 2010


FIT:    data read from 'oil15_131517.dat'
	x range restricted to [0.000100000 : 1000.00]
	y range restricted to [0.800000 : 2.00000]
        #datapoints = 479
        residuals are weighted equally (unit weight)

function used for fitting: f(x)
fitted parameters initialized with current variable values



 Iteration 0
 WSSR        : 22.9981           delta(WSSR)/WSSR   : 0
 delta(WSSR) : 0                 limit for stopping : 1e-05
 lambda	  : 0.51733

initial set of free parameter values

a               = 0.8
b               = 5
c               = 1
d               = 1
p               = 0.8

After 103 iterations the fit converged.
final sum of squares of residuals : 1.15725
rel. change during last iteration : -9.97087e-06

degrees of freedom    (FIT_NDF)                        : 474
rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 0.049411
variance of residuals (reduced chisquare) = WSSR/ndf   : 0.00244145

Final set of parameters            Asymptotic Standard Error
=======================            ==========================

a               = 0.834914         +/- 0.008833     (1.058%)
b               = 31.2275          +/- 14.74        (47.2%)
c               = 18.3346          +/- 41.12        (224.3%)
d               = 1.0012           +/- 0.002953     (0.295%)
p               = 1.25322          +/- 2.122        (169.3%)


correlation matrix of the fit parameters:

               a      b      c      d      p      
a               1.000 
b               0.408  1.000 
c              -0.296 -0.975  1.000 
d              -0.297  0.170 -0.221  1.000 
p              -0.347 -0.993  0.994 -0.189  1.000

What I'd like to do now, is to save on a separate file the values of the final set of parameter, formatted as follows:

a b c d p
0.83 31.2275 18.3346 1.0012 1.25322
values of the second cicle
third
ecc...

I thought the easiest way could be to assign to parameter values some bash parameters (l, m, ecc..), write these in the new output file, and using them again as starting parameters for the next fitting cycle. But, unfortunately, I actually don't know how to tell the script how to read the fit.log file.

Thank you in advance for any help,
Leonardo

theNbomr · 10-13-2010, 10:09 AM

The first part of your question seems easy enough to solve with something like awk or Perl.

I will rephrase the second part of your question to make sure I understand it correctly. Having acquired the coefficients of the fit function, you want to iteratively re-seed the initial values using the most recently acquired values.
From your sample script, it looks like you can use shell variables in your gnuplot script (I didn't know that; good to know). If so, I would use an approach that would assign to shell variables, one at a time, using the output of a parser that scans the fit.log. Your parser could take as a parameter, the desired coefficient name, and print its value to standard output:

Code:

a=`fitLogParser.pl a`
b=`fitLogParser.pl b`
# ... and so on...

You just need to write fitLogParser.pl (or whatever you want to call it). Do you require help to write such a parser (in any language; my weapon of choice is Perl)?
Does this make sense?

--- rod.

leonardo2887 · 10-14-2010, 03:13 AM

Yes, what you understood is right. Unfortunately I'm new to awk, and completely ignorant to perl. I would really like to learn some basics in programming for writing simple scripts... but I'm always missing time for it...

Turning to the problem again: I tried this, but it actually didn't work:

Code:

clear
for i in *.dat; do 
a = 0.8
b = 10
c = 1
#Initial set of values, before first fitting routine
gnuplot <<EOF
set logscale x
set xlabel 'Time [ms]'
set ylabel 'g(2)'
f1(x)=a*exp(-b*x)+c
a=$a
b=$b
c=$c
unset key
fit [0.0001:1000][0.95:2] f1(x) '$i' via a, b, c
set term svg
set out 'dls_$i.svg'
set pointsize 1
plot [0.0001:1000][0.95:2] f1(x) w l lw 2, '$i' ls 6
set term wxt
replot
EOF
a=`fitLogParser.pl a` #Something to copy parameter a
b=`fitLogParser.pl b`
c=`fitLogParser.pl c`
echo "$a;$b;$c" >> fitresults.txt
rm fit.log
done

It would probably be to simple to paste this few lines. Clearly, so as it is written, it doesn't work. That's probably why my outoput file is just filled with ";" and gnuplot claims that it get's no input data

Code:

gnuplot> a=
           ^
         line 0: constant expression required


gnuplot> b=
           ^
         line 0: constant expression required


gnuplot> c=
           ^
         line 0: constant expression required

I don't know how different and easy to implement perl or awk are. From my point of view, the best thing would be a simple line to write in the script something like

awk "a is the number after "a =" "

so that I can copy-paste it easily in different script, whit different number of parameters...

Thank you very much for your help!

theNbomr · 10-14-2010, 08:38 AM

You may have to simply export your shell variables:

Code:

export a=`fitLogParser.pl a`
export b=`fitLogParser.pl b`
export c=`fitLogParser.pl c`

Then, you should be able to use the exported variables in the gnuplot script:

Code:

f1(x)=$a*exp(-$b*x)+$c
# These are redundant
# a=$a
# b=$b
# c=$c
unset key
fit [0.0001:1000][0.95:2] f1(x) '$i' via $a, $b, $c

Here is a little bit of code that I used to test this (I just entered it at the commandline)

Code:

export start=2
export finish=6
gnuplot << EOF
plot [$start*pi:$finish*pi] sin(x)
EOF

Is your fitLogParser.pl doing anything? In order for the script to work as you ultimately want, it should be extracting the final fit coefficients from the gnufit output file, and printing the appropriate one to stdout. You can test the principle by having it simply print a random value if you haven't got the parser working yet.

--- rod.

theNbomr · 10-14-2010, 09:09 AM

Here's a version of the fit.log parser.

Code:

#! /usr/bin/perl -w
#
#	LQleonardo2887.pl
#	Usage:
#		LQleonardo2887.pl fit.log  a
#
use strict;
my %coefficients;

	# Find all lines that look like :
	#  a               = 0.834914         +/- 0.008833     (1.058%)
	# Use the trailing ')' to distinguish from other similar lines.
	#
	# Save all coefficients, then print the value of the 
	# one specified  on the commandline;

	open( FIT, shift );
	while(<FIT>){ 
		if( $_ =~ m/^([a-z])\s+=\s+([0-9+-.]+).+\)$/ ){
			$coefficients{$1} = $2;
			# print "$1 = $2\n";
		} 
	}
	close FIT;
	print $coefficients{$ARGV[0]},"\n";

There is no error checking, but it seems to work against your sample data.

--- rod.