Commit 539970d8 authored by Katter, Janike Yvonne's avatar Katter, Janike Yvonne
Browse files

Merge branch 'covarParerr' into 'develop'

Covar and parerr

See merge request !2
parents 3fb945e3 11125e73
Pipeline #38849 passed with stage
in 10 seconds
set(demos set(demos
curve1 curve1
curve2
surface1 surface1
nonlin1 nonlin1
) )
......
...@@ -13,7 +13,7 @@ ...@@ -13,7 +13,7 @@
* Homepage: https://jugit.fz-juelich.de/mlz/lmfit * Homepage: https://jugit.fz-juelich.de/mlz/lmfit
*/ */
#include "lmcurve_tyd.h" #include "lmcurve2.h"
#include <stdio.h> #include <stdio.h>
/* model function: a parabola */ /* model function: a parabola */
...@@ -27,6 +27,8 @@ int main() ...@@ -27,6 +27,8 @@ int main()
{ {
int n = 3; /* number of parameters in model function f */ int n = 3; /* number of parameters in model function f */
double par[3] = { 100, 0, -10 }; /* really bad starting value */ double par[3] = { 100, 0, -10 }; /* really bad starting value */
double parerr[3];
double covar[3*3];
/* data points: a slightly distorted standard parabola */ /* data points: a slightly distorted standard parabola */
int m = 9; int m = 9;
...@@ -41,7 +43,7 @@ int main() ...@@ -41,7 +43,7 @@ int main()
printf( "Fitting ...\n" ); printf( "Fitting ...\n" );
/* now the call to lmfit */ /* now the call to lmfit */
lmcurve_tyd( n, par, m, t, y, dy, f, &control, &status ); lmcurve2( n, par, parerr, covar, m, t, y, dy, f, &control, &status );
printf( "Results:\n" ); printf( "Results:\n" );
printf( "status after %d function evaluations:\n %s\n", printf( "status after %d function evaluations:\n %s\n",
...@@ -49,7 +51,7 @@ int main() ...@@ -49,7 +51,7 @@ int main()
printf("obtained parameters:\n"); printf("obtained parameters:\n");
for ( i = 0; i < n; ++i) for ( i = 0; i < n; ++i)
printf(" par[%i] = %12g\n", i, par[i]); printf(" par[%i] = %12g uncertainty = %12g\n", i, par[i], parerr[i]);
printf("obtained norm:\n %12g\n", status.fnorm ); printf("obtained norm:\n %12g\n", status.fnorm );
printf("fitting data as follows:\n"); printf("fitting data as follows:\n");
......
set(lib lmfit) set(lib lmfit)
set(${lib}_LIBRARY ${lib} PARENT_SCOPE) set(${lib}_LIBRARY ${lib} PARENT_SCOPE)
set(src_files lmcurve.c lmmin.c lminvert.c) set(src_files lmcurve.c lmcurve2.c lmmin.c lminvert.c)
set(inc_files lmcurve.h lmmin.h lmstruct.h lmdecls.h lmfit.hpp) set(inc_files lmcurve.h lmcurve2.h lmmin.h lmstruct.h lmdecls.h lmfit.hpp)
add_library(${lib} ${src_files}) add_library(${lib} ${src_files})
......
/* /*
* Library: lmfit (Levenberg-Marquardt least squares fitting) * Library: lmfit (Levenberg-Marquardt least squares fitting)
* *
* File: lmcurve.c * File: lmcurve2.c
* *
* Contents: Implements lmcurve_tyd(), a variant of lmcurve() that weighs * Contents: Implements lmcurve2(), a variant of lmcurve() that weighs
* data points y(t) with the inverse of the standard deviations dy. * data points y(t) with the inverse of the standard deviations dy.
* *
* Copyright: Joachim Wuttke, Forschungszentrum Juelich GmbH (2004-2013) * Copyright: Joachim Wuttke, Forschungszentrum Juelich GmbH (2004-2013)
...@@ -13,6 +13,8 @@ ...@@ -13,6 +13,8 @@
* Homepage: https://jugit.fz-juelich.de/mlz/lmfit * Homepage: https://jugit.fz-juelich.de/mlz/lmfit
*/ */
#include "lmcurve2.h"
#include <math.h>
#include "lmmin.h" #include "lmmin.h"
typedef struct { typedef struct {
...@@ -20,26 +22,30 @@ typedef struct { ...@@ -20,26 +22,30 @@ typedef struct {
const double* y; const double* y;
const double* dy; const double* dy;
double (*f)(const double t, const double* par); double (*f)(const double t, const double* par);
} lmcurve_tyd_data_struct; } lmcurve2_data_struct;
void lmcurve_tyd_evaluate( void lmcurve2_evaluate(
const double* par, const int m_dat, const void* data, double* fvec, const double* par, const int m_dat, const void* data, double* fvec,
int* info) int* info)
{ {
lmcurve_tyd_data_struct* D = (lmcurve_tyd_data_struct*)data; lmcurve2_data_struct* D = (lmcurve2_data_struct*)data;
int i; int i;
for (i = 0; i < m_dat; i++) for (i = 0; i < m_dat; i++)
fvec[i] = ( D->y[i] - D->f(D->t[i], par) ) / D->dy[i]; fvec[i] = ( D->y[i] - D->f(D->t[i], par) ) / D->dy[i];
} }
void lmcurve_tyd( void lmcurve2(
const int n_par, double* par, const int m_dat, const int n_par, double* par, double* parerr, double* covar,
const double* t, const double* y, const double* dy, const int m_dat, const double* t, const double* y, const double* dy,
double (*f)(const double t, const double* par), double (*f)(const double t, const double* par),
const lm_control_struct* control, lm_status_struct* status) const lm_control_struct* control, lm_status_struct* status)
{ {
lmcurve_tyd_data_struct data = { t, y, dy, f }; lmcurve2_data_struct data = { t, y, dy, f };
lmmin(n_par, par, m_dat, (const void*)&data, lmcurve_tyd_evaluate, lmmin2(n_par, par, NULL, covar, m_dat, NULL, (const void*)&data, lmcurve2_evaluate,
control, status); control, status);
int j;
if (parerr)
for (j = 0; j < n_par; j++)
parerr[j] = sqrt(covar[j*n_par+j]);
} }
/* /*
* Library: lmfit (Levenberg-Marquardt least squares fitting) * Library: lmfit (Levenberg-Marquardt least squares fitting)
* *
* File: lmcurve_tyd.h * File: lmcurve2.h
* *
* Contents: Declares lmcurve_tyd(), a variant of lmcurve() that weighs * Contents: Declares lmcurve2(), a variant of lmcurve() that weighs
* data points y(t) with the inverse of the standard deviations dy. * data points y(t) with the inverse of the standard deviations dy.
* *
* Copyright: Joachim Wuttke, Forschungszentrum Juelich GmbH (2004-2013) * Copyright: Joachim Wuttke, Forschungszentrum Juelich GmbH (2004-2013)
...@@ -15,23 +15,14 @@ ...@@ -15,23 +15,14 @@
#ifndef LMCURVETYD_H #ifndef LMCURVETYD_H
#define LMCURVETYD_H #define LMCURVETYD_H
#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
#define __BEGIN_DECLS extern "C" {
#define __END_DECLS }
#else
#define __BEGIN_DECLS /* empty */
#define __END_DECLS /* empty */
#endif
#include <lmstruct.h> #include "lmstruct.h"
__BEGIN_DECLS __BEGIN_DECLS
void lmcurve_tyd( LM_DLL void lmcurve2(
const int n_par, double* par, const int m_dat, const int n_par, double* par, double* parerr, double* covar,
const double* t, const double* y, const double* dy, const int m_dat, const double* t, const double* y, const double* dy,
double (*f)(double t, const double* par), double (*f)(double t, const double* par),
const lm_control_struct* control, lm_status_struct* status); const lm_control_struct* control, lm_status_struct* status);
......
...@@ -603,7 +603,7 @@ terminate: ...@@ -603,7 +603,7 @@ terminate:
if ( failure ) if ( failure )
goto no_error_estimate; goto no_error_estimate;
for (i = 0; i < m; i++) for (i = 0; i < m; i++)
fjac[j*m+i] = (wf[i] - fvec[i]) / step / S->fnorm; fjac[j*m+i] = (wf[i] - fvec[i]) / step;
x[j] = temp; /* restore */ x[j] = temp; /* restore */
} }
for (j = 0; j < n; j++) { for (j = 0; j < n; j++) {
...@@ -619,7 +619,7 @@ terminate: ...@@ -619,7 +619,7 @@ terminate:
goto no_error_estimate; goto no_error_estimate;
if (dx) if (dx)
for (j = 0; j < n; j++) for (j = 0; j < n; j++)
dx[j] = sqrt(wh2[j*n+j]); dx[j] = sqrt(wh2[j*n+j] * S->fnorm * S->fnorm / (m-n));
if (covar) if (covar)
memcpy(covar, wh2, n*n*sizeof(double)); memcpy(covar, wh2, n*n*sizeof(double));
goto end_error_estimate; goto end_error_estimate;
......
...@@ -25,15 +25,16 @@ endfunction() ...@@ -25,15 +25,16 @@ endfunction()
add_custom_target( add_custom_target(
man ALL man ALL
DEPENDS lmcurve.3 lmmin.3 lmmin2.3 lmfit.7 DEPENDS lmcurve.3 lmcurve2.3 lmmin.3 lmmin2.3 lmfit.7
) )
add_custom_target( add_custom_target(
html ALL html ALL
DEPENDS lmcurve.html lmmin.html lmmin2.html lmfit.html DEPENDS lmcurve.html lmcurve2.html lmmin.html lmmin2.html lmfit.html
) )
one_page(lmcurve 3) one_page(lmcurve 3)
one_page(lmcurve2 3)
one_page(lmmin 3) one_page(lmmin 3)
one_page(lmmin2 3) one_page(lmmin2 3)
one_page(lmfit 7) one_page(lmfit 7)
=pod
=begin html
<link rel="stylesheet" href="podstyle.css" type="text/css" />
=end html
=head1 NAME
lmcurve2 - Levenberg-Marquardt least-squares fit of a curve (t,y,dy)
=head1 SYNOPSIS
B<#include <lmcurve2.h>>
B<void lmcurve2(
const int> I<n_par>B<, double *>I<par>B<,
double *>I<parerr>B<, double *>I<covar>B<,
const int> I<m_dat>B<, constS< >double *>I<t>B<,
constS< >double *>I<y>B<, constS< >double *>I<dy>B<,
double (*>I<f>B<)( const double >I<ti>B<, const double *>I<par>B< ),
constS< >lm_control_struct *>I<control>B<,
lm_status_struct *>I<status>B<);>
B<extern const lm_control_struct lm_control_double;>
B<extern const lm_control_struct lm_control_float;>
B<extern const char *lm_infmsg[];>
B<extern const char *lm_shortmsg[];>
=head1 DESCRIPTION
B<lmcurve2()> wraps the more generic minimization function B<lmmin2()>, for use in curve fitting.
B<lmcurve2()> determines a vector I<par> that minimizes the sum of squared elements of a residue vector I<r>[i] := (I<y>[i] - I<f>(I<t>[i];I<par>)) / I<dy>[i]. Typically, B<lmcurve2()> is used to approximate a data set I<t>,I<y>,I<dy>, where I<dy> represents the standard deviation of empirical data I<y>, by a parametric function I<f>(I<ti>;I<par>). On success, I<par> represents a local minimum, not necessarily a global one; it may depend on its starting value. Users must ensure that all I<dy>[i] are positive.
Function arguments:
=over
=item I<n_par>
Number of free variables.
Length of parameter vector I<par>.
=item I<par>
Parameter vector.
On input, it must contain a reasonable guess.
On output, it contains the solution found to minimize ||I<r>||.
=item I<parerr>
Parameter uncertainties vector.
Array of length I<n_par> or B<NULL>.
On output, unless it or I<covar> is B<NULL>, it contains the weighted parameter uncertainties for the found parameters.
=item I<covar>
Covariance matrix.
Array of length I<n_par> * I<n_par> or B<NULL>.
On output, unless it is B<NULL>, it contains the covariance matrix.
=item I<m_dat>
Number of data points.
Length of vectors I<t>, I<y>, I<dy>.
Must statisfy I<n_par> <= I<m_dat>.
=item I<t>
Array of length I<m_dat>.
Contains the abcissae (time, or "x") for which function I<f> will be evaluated.
=item I<y>
Array of length I<m_dat>.
Contains the ordinate values that shall be fitted.
=item I<dy>
Array of length I<m_dat>.
Contains the standard deviations of the values I<y>.
=item I<f>
A user-supplied parametric function I<f>(ti;I<par>).
=item I<control>
Parameter collection for tuning the fit procedure.
In most cases, the default &I<lm_control_double> is adequate.
If I<f> is only computed with single-precision accuracy,
I<&lm_control_float> should be used.
Parameters are explained in B<lmmin2(3)>.
=item I<status>
A record used to return information about the minimization process:
For details, see B<lmmin2(3)>.
=back
=head1 EXAMPLE
Fit a data set y(x) with standard deviations dy(x) by a curve f(x;p):
#include "lmcurve2.h"
#include <stdio.h>
/* model function: a parabola */
double f( double t, const double *p )
{
return p[0] + p[1]*t + p[2]*t*t;
}
int main()
{
int n = 3; /* number of parameters in model function f */
double par[3] = { 100, 0, -10 }; /* really bad starting value */
double parerr[3];
double covar[3*3];
/* data points: a slightly distorted standard parabola */
int m = 9;
int i;
double t[9] = { -4., -3., -2., -1., 0., 1., 2., 3., 4. };
double y[9] = { 16.6, 9.9, 4.4, 1.1, 0., 1.1, 4.2, 9.3, 16.4 };
double dy[9] = { 4, 3, 2, 1, 2, 3, 4, 5, 6 };
lm_control_struct control = lm_control_double;
lm_status_struct status;
control.verbosity = 1;
printf( "Fitting ...\n" );
/* now the call to lmfit */
lmcurve2( n, par, parerr, covar, m, t, y, dy, f, &control, &status );
printf( "Results:\n" );
printf( "status after %d function evaluations:\n %s\n",
status.nfev, lm_infmsg[status.outcome] );
printf("obtained parameters:\n");
for ( i = 0; i < n; ++i)
printf(" par[%i] = %12g uncertainty = %12g\n", i, par[i], parerr[i]);
printf("obtained norm:\n %12g\n", status.fnorm );
printf("fitting data as follows:\n");
for ( i = 0; i < m; ++i)
printf(
" t[%1d]=%2g y=%5.1f+-%4.1f fit=%8.5f residue=%8.4f weighed=%8.4f\n",
i, t[i], y[i], dy[i], f(t[i],par), y[i] - f(t[i],par),
(y[i] - f(t[i],par))/dy[i] );
return 0;
}
=head1 COPYING
Copyright (C) 2009-2015 Joachim Wuttke, Forschungszentrum Juelich GmbH
Software: FreeBSD License
Documentation: Creative Commons Attribution Share Alike
=head1 SEE ALSO
=begin html
<a href="http://apps.jcns.fz-juelich.de/man/lmmin2.html"><b>lmmin2</b>(3)</a>
<p>
Homepage: <a href="https://jugit.fz-juelich.de/mlz/lmfit">https://jugit.fz-juelich.de/mlz/lmfit</a>
=end html
=begin man
\fBlmmin2\fR(3)
.PP
Homepage: https://jugit.fz-juelich.de/mlz/lmfit
=end man
=head1 BUGS
Please send bug reports and suggestions to the author <j.wuttke@fz-juelich.de>.
...@@ -17,12 +17,16 @@ B<lmfit> is a C library for Levenberg-Marquardt least-squares minimization and c ...@@ -17,12 +17,16 @@ B<lmfit> is a C library for Levenberg-Marquardt least-squares minimization and c
For fitting a data set {(x_i,y_i)|i=0,1,..} by a parametric curve f(x,t), see B<lmcurve>(3). For fitting a data set {(x_i,y_i)|i=0,1,..} by a parametric curve f(x,t), see B<lmcurve>(3).
For fitting a data set {(x_i,y_i+-dy_i)|i=0,1,..} by a parametric curve f(x,t), see B<lmcurve2>(3).
For generic minimization of the Eucledian norm of parametric vector, see B<lmmin2>(3). For generic minimization of the Eucledian norm of parametric vector, see B<lmmin2>(3).
For the simpler legacy API without error estimates, see B<lmmin>(3). For the simpler legacy API without error estimates, see B<lmmin>(3).
For an example how to use B<lmmin>, see the source files I<lmcurve.h> and I<lmcurve.c>. Do not patch these files; copy and modify them to create your own, differently named version of I<lmcurve_data_struct>, I<lmcurve_evaluate>, and I<lmcurve>. For an example how to use B<lmmin>, see the source files I<lmcurve.h> and I<lmcurve.c>. Do not patch these files; copy and modify them to create your own, differently named version of I<lmcurve_data_struct>, I<lmcurve_evaluate>, and I<lmcurve>.
For an example how to use B<lmmin2> for weighted data, see the source files I<lmcurve2.h> and I<lmcurve2.c>. Do not patch these files; copy and modify them to create your own, differently named version of I<lmcurve2_data_struct>, I<lmcurve2_evaluate>, and I<lmcurve2>.
=head1 COPYING =head1 COPYING
Copyright (C): Copyright (C):
...@@ -48,7 +52,7 @@ Homepage: <a href="https://jugit.fz-juelich.de/mlz/lmfit">https://jugit.fz-jueli ...@@ -48,7 +52,7 @@ Homepage: <a href="https://jugit.fz-juelich.de/mlz/lmfit">https://jugit.fz-jueli
=begin man =begin man
\fBlmcurve\fR(3), \fBlmmin\fR(3), \fBlmmin2\fR(3) \fBlmcurve\fR(3), \fBlmcurve2\fR(3), \fBlmmin\fR(3), \fBlmmin2\fR(3)
.PP .PP
Homepage: https://jugit.fz-juelich.de/mlz/lmfit Homepage: https://jugit.fz-juelich.de/mlz/lmfit
......
...@@ -58,14 +58,14 @@ On output, it contains the solution found to minimize ||I<fvec>||. ...@@ -58,14 +58,14 @@ On output, it contains the solution found to minimize ||I<fvec>||.
=item I<parerr> =item I<parerr>
Parameter error vector, either of length I<n_par>, or B<NULL>. Parameter error vector, either of length I<n_par>, or B<NULL>.
On output, unless it is B<NULL>, it contains the parameter uncertainties, On output, unless it is B<NULL>, it contains the parameter uncertainties for the unweighted model,
estimated from the diagonal elements of the covariance matrix. estimated from the diagonal elements of the covariance matrix and the sum of squared elements of I<fvec>-I<y> divided by I<m_dat>-I<n_par>.
If your data is weighed, use square root of the diagonal elements of the covariance matrix for the parameter uncertainties.
=item I<covar> =item I<covar>
Covariance matrix, stored as vector of length I<n_par>*I<n_par>, or B<NULL>. Covariance matrix, stored as vector of length I<n_par>*I<n_par>, or B<NULL>.
On output, unless it is B<NULL>, it contains the covariance matrix. On output, unless it is B<NULL>, it contains the covariance matrix.
NOT YET IMPLEMENTED!
=item I<m_dat> =item I<m_dat>
...@@ -316,6 +316,7 @@ Homepage: <a href="https://jugit.fz-juelich.de/mlz/lmfit">https://jugit.fz-jueli ...@@ -316,6 +316,7 @@ Homepage: <a href="https://jugit.fz-juelich.de/mlz/lmfit">https://jugit.fz-jueli
\fBlmmin\fR(3), \fBlmmin\fR(3),
\fBlmcurve\fR(3) \fBlmcurve\fR(3)
\fBlmcurve2\fR(3)
.PP .PP
Homepage: https://jugit.fz-juelich.de/mlz/lmfit Homepage: https://jugit.fz-juelich.de/mlz/lmfit
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment