fitting by fortan 3

[Tony1984]

Hey
i need to use fortan to make a program for fitting data,
in my project i have an experimental data and theoritical one , i need to fit this 2 kinds of data's by fortran, suppose we have a seperate file called "DATA" contains 2 colums one for the expiremtal and the other for the theoritical one,
Notice: Im dealing with a huge Number of data, some times up to 100 thousands of experimental and theoritical data needes to be fit.
Regards


File "DATA" contian:
Exp Theory

5.293 5.639
5.346 5.692
5.399 5.745
5.454 5.8
5.509 5.855
5.564 5.91
5.62 5.966
5.676 6.022
5.733 6.079
5.791 6.137
5.849 6.195
5.908 6.254
5.967 6.313
6.027 6.373
6.088 6.434
6.149 6.495
6.211 6.557
6.273 6.619
6.336 6.682

 

[gummibaer]

What do you mean by 'fitting' ?

There are a lot of fortran callable software packages around that offer statistical software and all kinds of statistical tools to analyse your data, compute regressions and all such. Google for 'fortran statistics library' and you will find more references to such packages than you will ever need.

For the data you added, I'd say there is an offset of roughly 0.345 between the two columns and that is it.

Norbert


The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[Tony1984]

Dear Gummibaer
the data above is not a real data i just write it down as an example,
the fitting procedure includes ; least square method, generlized least square method ,and bayes method, so if you familiar of one of these methods please let me know.

Thanks

 

[mikrom]

Fortran doesn't have build in function for this. You have to program it self, or search if there is a procedure available which you can use.

I'm not sure if I undestand your problem right:
Your data file has only 2 columns: experimental measured value and theoretical value. I guess, you want to fit the measured value and compare them with theoretical values. But you have not independent variable in your data file. Your measured value must depend of something, maybe of time ... IMO you need to create your datafile something like this:

t1, exp_t1, theory_t1
t2, exp_t2, theory_t2
...
...
tN, exp_tN, theory_tN

That is: you have in the 1-st column independent variable - for example time values
in the 2-nd column experimental velues, which depend on 1-st column - for example values measured at given time
and if you want, you can have the 3-th column which contain theoretical values - for example computed according to a given formula.

The least squares method is described for here:

http://en.wikipedia.org/wiki/Least_squares

It seems that you have only one dimensional case. When you only need linear fitting, i.e. you have to find the coefficients a, b of the linear function y = a*x + b, which aproximate best your data, then it's simple to derive the formulas self and to program it.

 

[gummibaer]

Look here:

http://www.itl.nist.gov/div898/software/reeves/homepage.htm

and if this does not appeal to you, check here for alternatives:

http://fortranwiki.org/fortran/show/Libraries

Norbert


The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[ArkM]

The only case when no need for independent variable(s): if both experimental and theoretical data are defined on the uniform grid. No matter what's this grid origin or interval.

 

[Tony1984]

Dear mikrom
both types of data are function of energy.

 

[mikrom]

both types of data are function of energy

And what you want to do? To fit a curve which describes the dependency between experimental values and energy, or to fit a curve which describes the dependency between experimental and theoretical values ?

 

[mikrom]

Notice: Im dealing with a huge Number of data, some times up to 100 thousands of experimental and theoritical data needes to be fit.

Linear fit seems to be adequate for your data. It's very simple to program without dealing with large arrays:
linear_fit.f95

[COLOR=#a020f0]program[/color] linear_fit
  [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
  [COLOR=#0000ff]! Fitting line y = a*x + b to data points (Least Squares Method)[/color]
  
  [COLOR=#2e8b57][b]integer[/b][/color] :: n 
  [COLOR=#2e8b57][b]double precision[/b][/color] :: x, y, [COLOR=#804040][b]&[/b][/color] [COLOR=#0000ff]! data points[/color]
                      Sx,   [COLOR=#804040][b]&[/b][/color] [COLOR=#0000ff]! Sx  = x1 + x2 + ..[/color]
                      Sy,   [COLOR=#804040][b]&[/b][/color] [COLOR=#0000ff]! Sy  = y1 + y2 + ..[/color]
                      Sxx,  [COLOR=#804040][b]&[/b][/color] [COLOR=#0000ff]! Sxx = x1*x1 + x2*x2 + ..[/color]
                      Sxy,  [COLOR=#804040][b]&[/b][/color] [COLOR=#0000ff]! Sxy = x1*y1 + x2*y2 + ...[/color]
                      a, b, [COLOR=#804040][b]&[/b][/color] [COLOR=#0000ff]! coefficients of line: y = a*x + b[/color]
                      Sres, [COLOR=#804040][b]&[/b][/color] [COLOR=#0000ff]! Sum of squared residuals[/color]
                      rr,   [COLOR=#804040][b]&[/b][/color] [COLOR=#0000ff]! coefficient of determination: rr = Sc/Sd[/color]
                      Sf,   [COLOR=#804040][b]&[/b][/color] [COLOR=#0000ff]! Sum of squared fiffrences betwen fit and mean [/color]
                      Sd,   [COLOR=#804040][b]&[/b][/color] [COLOR=#0000ff]! Sum of squared fiffrences betwen data and mean [/color]
                      y_mean  [COLOR=#0000ff]! y_mean = Sy/n[/color]

  [COLOR=#0000ff]! open file[/color]
  [COLOR=#804040][b]open[/b][/color] ([COLOR=#804040][b]unit[/b][/color] [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]10[/color], [COLOR=#804040][b]file[/b][/color] [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]'linear_fit_data.txt'[/color], [COLOR=#804040][b]status[/b][/color] [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]'unknown'[/color])

  n [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]
  x [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]; y [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]
  Sx [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]; Sy [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]; Sxx [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]; Sxy [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]
  [COLOR=#0000ff]! read data and compute sums[/color]
  [COLOR=#804040][b]do[/b][/color] [COLOR=#804040][b]while[/b][/color] ([COLOR=#ff00ff].true.[/color])
    [COLOR=#804040][b]read[/b][/color] ([COLOR=#ff00ff]10[/color], [COLOR=#804040][b]*[/b][/color], [COLOR=#804040][b]end[/b][/color][COLOR=#804040][b]=[/b][/color][COLOR=#804040][b]99[/b][/color]) x, y
    [COLOR=#0000ff]!write(*,*) x, y[/color]
    n [COLOR=#804040][b]=[/b][/color] n [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]
    Sx [COLOR=#804040][b]=[/b][/color] Sx [COLOR=#804040][b]+[/b][/color] x
    Sy [COLOR=#804040][b]=[/b][/color] Sy [COLOR=#804040][b]+[/b][/color] y
    Sxx [COLOR=#804040][b]=[/b][/color] Sxx [COLOR=#804040][b]+[/b][/color] x[COLOR=#804040][b]*[/b][/color]x
    Sxy [COLOR=#804040][b]=[/b][/color] Sxy [COLOR=#804040][b]+[/b][/color] x[COLOR=#804040][b]*[/b][/color]y
  [COLOR=#804040][b]end do[/b][/color]
[COLOR=#804040][b]  99 continue[/b][/color] [COLOR=#0000ff]! end of file[/color]

  [COLOR=#0000ff]! compute the line coefficients[/color]
  a [COLOR=#804040][b]=[/b][/color] (n[COLOR=#804040][b]*[/b][/color]Sxy [COLOR=#804040][b]-[/b][/color] Sx[COLOR=#804040][b]*[/b][/color]Sy)[COLOR=#804040][b]/[/b][/color](n[COLOR=#804040][b]*[/b][/color]Sxx [COLOR=#804040][b]-[/b][/color] Sx[COLOR=#804040][b]**[/b][/color][COLOR=#ff00ff]2[/color])
  b [COLOR=#804040][b]=[/b][/color] (Sy [COLOR=#804040][b]-[/b][/color] a[COLOR=#804040][b]*[/b][/color]Sx)[COLOR=#804040][b]/[/b][/color]n 

  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color], [COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Total number of data points = '[/color], n
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color])  [COLOR=#ff00ff]'Line computed: y = '[/color], a, [COLOR=#ff00ff]'*x + '[/color], b
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color], [COLOR=#804040][b]*[/b][/color])
  
  [COLOR=#0000ff]! now reading the data again to determine quality of the fit[/color]
  [COLOR=#804040][b]rewind[/b][/color]([COLOR=#ff00ff]10[/color])
  x [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]; y [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]
  Sres [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]; Sf [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]; Sd [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]
  y_mean [COLOR=#804040][b]=[/b][/color] Sy[COLOR=#804040][b]/[/b][/color]n
  [COLOR=#804040][b]do[/b][/color] [COLOR=#804040][b]while[/b][/color] ([COLOR=#ff00ff].true.[/color])
    [COLOR=#804040][b]read[/b][/color] ([COLOR=#ff00ff]10[/color], [COLOR=#804040][b]*[/b][/color], [COLOR=#804040][b]end[/b][/color][COLOR=#804040][b]=[/b][/color][COLOR=#804040][b]999[/b][/color]) x, y
    [COLOR=#0000ff]!write(*,*) x, y[/color]
    Sres [COLOR=#804040][b]=[/b][/color] Sres [COLOR=#804040][b]+[/b][/color] (a[COLOR=#804040][b]*[/b][/color]x [COLOR=#804040][b]+[/b][/color] b [COLOR=#804040][b]-[/b][/color] y)[COLOR=#804040][b]**[/b][/color][COLOR=#ff00ff]2[/color]
    Sf [COLOR=#804040][b]=[/b][/color] Sf [COLOR=#804040][b]+[/b][/color] (a[COLOR=#804040][b]*[/b][/color]x [COLOR=#804040][b]+[/b][/color] b [COLOR=#804040][b]-[/b][/color] y_mean)[COLOR=#804040][b]**[/b][/color][COLOR=#ff00ff]2[/color]
    Sd [COLOR=#804040][b]=[/b][/color] Sd [COLOR=#804040][b]+[/b][/color] (y [COLOR=#804040][b]-[/b][/color] y_mean)[COLOR=#804040][b]**[/b][/color][COLOR=#ff00ff]2[/color]
  [COLOR=#804040][b]end do[/b][/color]
[COLOR=#804040][b]  999 continue[/b][/color] [COLOR=#0000ff]! end of file[/color]
  
  [COLOR=#0000ff]! compute coefficient of determination[/color]
  rr [COLOR=#804040][b]=[/b][/color] Sf[COLOR=#804040][b]/[/b][/color]Sd
  
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color], [COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Quality of the fit:'[/color]  
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color], [COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Sum of squared residuals     = '[/color], Sres
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color], [COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Coefficient of determination = '[/color], rr

  [COLOR=#0000ff]! close file[/color]
  [COLOR=#804040][b]close[/b][/color]([COLOR=#ff00ff]10[/color])
[COLOR=#a020f0]end program[/color] linear_fit

I have read the file twice to determine quality of the fit, i.e. sum of squares and the coefficient of determination. If you don't need these values, you have to read the file only once.

Here is the output:

$ gfortran linear_fit.f95 -o linear_fit

$ linear_fit
 Total number of data points =           19
 Line computed: y =   0.99999999999996958      *x +   0.34600000000017855     

 Quality of the fit:
 Sum of squared residuals     =   1.85446285405266586E-027
 Coefficient of determination =   0.99999999999993949

It seems that the line adequately fits the given data set, because
Sum of squares -> 0 and Coefficient of determination -> 1

 

[salgerman]

Except, as mentioned before, that Tony seems to have provided two sets of Y's and no X's; so, either X's need to be provided or a constant step may be assumed, i.e., X's = 0,1,2,3,4...

Assuming X's are 0,1,2,3...18 and realizing the data is pretty straight, a first approximation to the line would be:

b = first data point = 5.64
m = (last - first) / span = (6.68 - 5.64)/18 = 0.05777

Fit: y = 0.05778*x + 5.64

Or, running the program above:

Total number of data points = 19
Line computed: y = 0.05782456140350794 *x + 5.625894736842114

Quality of the fit:
Sum of squared residuals = 0.0013445614035087767
Coefficient of determination = 0.9992950232171091

 

[mikrom]

But maybe Tony only wants to look at the dependecy between experimental value and theoretical value - for example to adjust his formula for theoretical compuation.
Then from the result above, it seems to be: theoretical_value = experimental_value + 0.346, what was already estimated in the first answer by gummibär.