Dear all,

Is it possible to send me an example for fitting of a non-linear equation to a matrix of data? Thanks

Dear all,

Is it possible to send me an example for fitting of a non-linear equation to a matrix of data? Thanks

Hello pl try to impose the following code for fgauss (user-supplied function) with the routine mrqmin(inturn using mrqcof,covsrt and gaussj) which after sandwiching ,it fits the rquired model NR(Ch15)
Also this fg class below requires a little more debugging before use especially so for the accessing of vector/objects .Thanks


#include "nr3.h"
#include <math.h>
#define myran
#define Ran

using namespace std;
void fgauss(const Doub x, VecDoub_I &a, Doub &y, VecDoub_O &dyda);
void Nonlinear(void);
class fg{

private:
const Doub x,y;
VecDoub_I a;
VecDoub_O dyda;


public :
void fgauss(void);
void fgauss(const Doub x, VecDoub_I &a, Doub &y, VecDoub_O &dyda) {
Int i,na=a.size();
Doub fac,ex,arg;
y=0.;
for (i=0;i<na-1;i+=3) {
arg=(x-a[i+1])/a[i+2];
ex=exp(-SQR(arg));
fac=a[i]*ex*2.*arg;
y += a[i]*ex;
dyda[i]=ex;
dyda[i+1]=fac/a[i+2];
dyda[i+2]=fac*arg/a[i+2];
}
};
};
int main()
{
Ran myran(12345678); // For debugging, use the same seed every time
//Ran myran(time(0)); // For exploration, use different seeds

Int n =10;

VecDoub x(n), y(n);
VecDoub_I a(n);
VecDoub_O dyda(n);

for (int i = 0; i < n; i++) {
x[i] = myran.doub();
y[i] = myran.doub();
a[i] = myran.doub();
dyda[i] = myran .doub();
// fgauss( Doub , VecDoub_I , Doub , VecDoub_O );
fgauss( x[i],y[i],a[i],dyda[i]);


}
return 0;
}

The following is the Algorith ref from Text,

The Example below can be tried for more clarity on non Linear type of Matrix problems




Quassi_Newton algorithm


Step:0)

Chose x^0 and B.0(x^0 should be in the doman of atttraction of x*^8).
For k=0,1,2 ....., do step 1 - 4 untill || f^k|| < E,

where E is sometolerence chosen on the basis of rounding errors

and || || denotes The Euclidean norm.

Step: 1) Solve B.(k)s^k =-f^k for s^k

Step: 2) Choose t(k) from a reasonably small set of t values on the basis of some apprximate line such that
||f(x^k+t(k)s^k||= mtin||f(x^k + ts^k)||
If x^k is close to x^8 ,then t(k) can be taken as unity.

Step:3)

Compute x^(k+1),y^k,p^k from
x^(k+1= x^k +t(k)s^k,

y^k= f^(k+1) - f^k,
p^k =(y^k/t(k))-B(k)s^k.

Step:4)
Obtain the next approximation to the Jacobian as
B(k+1) = B(k) + triangleB(k),
where
triangleB(k)=p^ks^k^t/s^k^ts^k
is a solution of triangle B(k)s^k.


Use this algorithm to solve the following two type of problems

1)

fi=(3-klxi)xi+1-xi-1-2xi+1,
i=l(l)n.


2)
fi=(kl+k2xi^2)xi+l-k(3){xj+xj^2) {while (j=i-r1,(j!= i)) and tends to (i+r2)}.
i=l(l)n

For (1)
n k(l)
--- ----

5 0.1

5 0.5

10 0.5

10 0.5

20 0.5

600 0.5

600 2.0

In this problem,

xl , i< l or i > n

is regarded as zero ; also the parameters

k1,k2,k3 allow the amount of nonlinearity to be varied.
Take

a) the initial estimate (for the solution) xi = -l,i=l(l)n,

b) the final tolerance on || f|| = 10^(-6)

c) he step length t(k) =l and

d) the initial approximation B(0)(The Jocobian) ,computed
using
the forward differences(with 0.001 as the increment of each variable)

Thanks Kutta,
I think that your code can be implemented for non-linear optimization of a matrix of data and I am trying to use it in my program. However, the code wrote is basically refers to Ed2 b/c of fgaus, gaussj, covsrt, mrqcof,mrqmin routines. Is it a way to use the new routines in 3rd Ed so that everything would be more consistent with nr3.h? Thanks again for your email.

Ed2...fgaus, gaussj, covsrt, mrqcof,mrqmin routines --->3rd Ed...



//
// Demonstration of use of Numerical Recipes Fitmrq.
//
// The original function is the sum of two Gaussian functions,
// and the data points are corrupted by noise.
//
// The Numerical Recipes distribution conveniently supplies a
// function that can be fed to Fitmrq to try to model data that
// consists of samples from sum of any number of Gaussian-type functions.
//
// I corrupt the data values with additive Gaussian noise.
//
// davekw7x
//


// Always include nr3.h first
#include "../code/nr3.h"

// For Normaldev
#include "../code/ran.h"
#include "../code/gamma.h"
#include "../code/deviates.h"

// For Fitmrq
#include "../code/gaussj.h"
#include "../code/fitmrq.h"

// The fgauss() function is in the file fit_examples.h
#include "../code/fit_examples.h"

//
// Simplify appearance of expressions with x*x by using an in-line function
//
inline Doub sqr(const Doub & x)
{
return x * x;
}

int main()
{
int i;

Doub noise = 0.1; // Standard deviation of generated noise
Normaldev ran(0, noise, 8654321); // Gaussian noise generator
//Normaldev ran(0, noise, time(0)); // Gaussian noise generator

Int num_points = 101;

// In order to use fgauss on K Gaussians, we need
// some vectors with size 3*K. So, for two Gaussians
// we have size = 6
//
// The two groups of three parameters contain amplitude, mean,
// and standard deviation of a Gaussian function that is
// used to generate the data points for testing.
//
// This is the actual array of parameters for the two Gaussian
// functions:
// [first group], [second group]
const Doub actual_d[] = { 5.0, 2.0, 3.0, 2.0, 5.0, 3.0 };

//
// This is the "guess" that we supply to fitmrq:
//
// Try it with a "pretty good" guess:
const Doub guess_d[] = { 4.8, 2.2, 2.8, 2.2, 4.8, 3.2 };

//
// Just for kicks, you might want to try it with a "not so pretty
// good" guess. Maybe you will get lucky; maybe not.
// The point is that Levenberg-Marquardt requires an initial
// guess of some kind.
//
//const Doub guess_d[] = { 2.0, 2.0, 2.0, 1.0, 1.0, 1.0 };

//
// The number of elements in the array is needed to initialize
// a VecDoub object from that array.
//
const Int ma = sizeof(actual_d) / sizeof(actual_d[0]);

VecDoub actual(ma, actual_d);
VecDoub guess(ma, guess_d);
VecDoub x(num_points);
VecDoub y(num_points);
VecDoub sd(num_points);

//
// Fill array y with values for sum of two gaussians, corrupted
// by additive gaussian noise.
//
// The sd array contains estimates of the standard deviations
// of the elements of y. We use a "measurement error" factor
// equal to the standard deviation of the noise distribution.
//
Doub xstart = 0.0;
Doub xmax = 10.0;
Doub deltax = xmax / (num_points - 1);
Doub xx;
for (i = 0, xx = xstart; i < num_points; i++, xx += deltax) {
x[i] = xx;
y[i] = ran.dev() + (actual[0]*exp(-sqr((xx - actual[1]) / actual[2])) +
actual[3]*exp(-sqr((xx - actual[4]) / actual[5])));
sd[i] = noise * y[i];
}

//
// Use the constructor to bind the data
//
Fitmrq mymrq(x, y, sd, guess, fgauss); // Use default tolerance 1.0e-3

//
// Now do the fit
//
mymrq.fit();

cout << scientific << setprecision(5);
cout << setw(10) << "Actual"
<< setw(13) << "Fitted"
<< " "
<< setw(13) << "Err. Est."
<< endl;

for (i = 0; i < actual.size(); i++) {
cout << setw(13) << actual[i]
<< setw(13) << mymrq.a[i]
<< " +- " << sqrt(mymrq.covar[i][i])
<< endl;
}
cout << endl;

cout << fixed;
cout << "Number of points = " << num_points
<< ", Chi-squared = " << mymrq.chisq << endl;

return 0;
}



Output on my Linux workstation (g++ version 4.1.2)



Actual Fitted Err. Est.
5.00000e+00 5.52170e+00 +- 6.84019e-01
2.00000e+00 2.14463e+00 +- 4.97134e-01
3.00000e+00 2.93038e+00 +- 4.86329e-01
2.00000e+00 1.78320e+00 +- 1.43392e+00
5.00000e+00 5.46802e+00 +- 9.33546e-01
3.00000e+00 2.56607e+00 +- 2.90910e-01

Number of points = 101, Chi-squared = 309.83746


Regards,

Dave

Footnotes:
1. I think that chances of getting helpful responses with fewest iterations might be improved if requestors would give some kind of problem statement in the original post. (But there are no guarantees.)

2. I make no claims as to the suitability or usability of the procedure or of the appropriateness of any assumptions about noise or anything else; my intent is to show how to invoke the tools. Interpretation of results is up to the user. (In other words: Your Mileage May Vary.)

Thanks Dave,
Sure, I will add my program and the problem statement in my future questions.

Thanks Kutta,
I think that your code can be implemented for non-linear optimization of a matrix of data and I am trying to use it in my program. However, the code wrote is basically refers to Ed2 b/c of fgaus, gaussj, covsrt, mrqcof,mrqmin routines. Is it a way to use the new routines in 3rd Ed so that everything would be more consistent with nr3.h? Thanks again for your email.
Hello Try this for more clarity of application
Thanking you
As
Kutta(C.R.Muthukumar):)

Dear all,

Is it possible to send me an example for fitting of a non-linear equation to a matrix of data? Thanks

Hello On the simplest type ,the homogenius equations (similar to non linear eqn) finds a place using cramers algorithm to obtain determinants
shown for ur perusal:-

Thanking You
As
Kutta(C.R.Muthukumar)
:)

Hi Dave and Kutta,

The example that you have provided had the fgauss function as an example which is basically a single variable function. Now, my situation is different. My function has two variables, i.e., y=f(x,eta)=> x and eta matrices should be used as an input in the function that I should write in fit_examples.h. I am confused about implementing of fitmrq.h into my program. I have attached my codes along with the req'd txt files. Note that the experimental (simulated) values are A0Sim and the calculated values are A0Calc.

Cordially,
Amir

I need help on multi-variable function as discussed above. Thanks

I need help on multi-variable function as discussed above. Thanks


Invoke the equation Z(x)=f(x) + {K(x,t)F(t,Z(t))dt}
from a to b .Here the unknown function F(t,Z(t)) may be any function of Z(t) such as[Z(t)]^3 and tan[Z(t)].
In such a case nonlinear equations are produced first instead of linear equaitions.These equations can then be solved by a suitable iterative technique such as the Newton Method with good initiall approximation.This is the required method for multivariables
However
from the contextual point Pl find an good method (SOR)included that needs a bit of fine tuning to make it workable/executable
Thanking You
As
Kutta(C.R.Muthukumar):)




#include "nr3.h"
#include <iomanip.h>
#include <iostream>
#include <fstream>
#include <math.h>

using namespace std;

//void printMatDoub(MatDoub_I & a);
//void printMatDoub(MatDoub_I & b);
//void printMatDoub(MatDoub_I & c);
//void printMatDoub(MatDoub_I & d);
//void printMatDoub(MatDoub_I & e);
//void printMatDoub(MatDoub_I & f);
//void printMatDoub(MatDoub_IO & u);

Doub printMatDoub(MatDoub_I & a , MatDoub_I & b, MatDoub_I & c, MatDoub_I & d, MatDoub_I & e, MatDoub_I & u ,MatDoub_I &f,const Doub rjac);

//typedef printMatDoub(myMatDoub);
Doub printMatDoub();
void sor(MatDoub_I &a, MatDoub_I &b, MatDoub_I &c, MatDoub_I &d, MatDoub_I &e ,MatDoub_I &f, MatDoub_IO &u, const Doub rjac)
{


const Int MAXITS=1000;
const Doub EPS=1.0e-13;
Doub anormf=0.0,omega=1.0;
Int jmax=a.nrows();
for (Int j=1;j<jmax-1;j++)
for (Int l=1;l<jmax-1;l++)
anormf += abs(f[j][l]);
for (Int n=0;n<MAXITS;n++) {
Doub anorm=0.0;
Int jsw=1;
for (Int ipass=0;ipass<2;ipass++) {
Int lsw=jsw;
for (Int j=1;j<jmax-1;j++) {
for (Int l=lsw;l<jmax-1;l+=2) {
Doub resid=a[j][l]*u[j+1][l]+b[j][l]*u[j-1][l]
+c[j][l]*u[j][l+1]+d[j][l]*u[j][l-1]
+e[j][l]*u[j][l]-f[j][l];
anorm += abs(resid);
u[j][l] -= omega*resid/e[j][l];
}
lsw=3-lsw;
}
jsw=3-jsw;
omega=(n == 0 && ipass == 0 ? 1.0/(1.0-0.5*rjac*rjac) :
1.0/(1.0-0.25*rjac*rjac*omega));
}
if (anorm < EPS*anormf) return;
}
throw("MAXITS exceeded");
// printMatDoub(myMatDoub);
}



int main()
{
MatDoub myMatDoub;
MatDoub_I a7,b,c,d,e,f,u;

// Define in row-major order
{
Doub a[3*2] = {1, 2, 3, 4, 5, 6};
MatDoub myMatDoub(3, 2, a);
}
{Doub b[3*2] = {1, 2, 3, 4, 5, 6};
MatDoub myMatDoub(3, 2, b);}
{
Doub c[3*2] = {1, 2, 3, 4, 5, 6};
MatDoub myMatDoub(3, 2, c);
}
{ Doub d[3*2] = {1, 2, 3, 4, 5, 6};
MatDoub myMatDoub(3, 2, d);
}
{ Doub e[3*2] = {1, 2, 3, 4, 5, 6};
MatDoub myMatDoub(3, 2, e);
}
{ Doub f[3*2] = {1, 2, 3, 4, 5, 6};
MatDoub myMatDoub(3, 2, f);
}
{ Doub u[3*2] = {1, 2, 3, 4, 5, 6};
MatDoub myMatDoub(3, 2, u);
}
//printMatDoub();

cout << " \nPrintMatDoub: ";
cout << "Number of rows = " << a.nrows()
<< ", Number of cols = " << a.ncols()<< endl;
for (int i = 0; i < a.nrows(); i++) {
cout << " ";
for (int j = 0; j < a.ncols(); j++) {
cout << "a[" << i << "][" << j << "] = " << a[i][j] << " ";
}
cout << endl;


return 0;
}
}

In the program that I have provided, the desired variable to be optimized are ksij0 and ksij1. The objective function (to be minimized) as I mentioned before are (A0Sim-A0Calc).

In the program that I have provided, the desired variable to be optimized are ksij0 and ksij1. The objective function (to be minimized) as I mentioned before are (A0Sim-A0Calc).

Hello Please try the following routine that I am sure would help you with a result though it requires some bit of initialisation,synchronsing with your three goals viz
ksij0,ksij1,the ojectivefunction which are as similar to this
Please do not get confused and a little bit of initiation will enable you to win for a concrete solution since this is from NR
However if this does not lead you so ,please try the ones under NonLinear Models whose routine name is mrqmin(with many variables).method
Thanking you
As
Kutta(C.R.Muthukumar):)

#include "nr3.h"
#define chisq
#define NR
using namespace std;

"void" NR::lfit(Vec_I_DP &x,Vec_I_Dp &y,Vec_I_DP &sig,,Vec_I_BOOL, & ia,Mat_O_DP &cover,DP &chisq,void funcs(const DP,Vec_O_DP &))
{
int i,j,k,l,m,mfit=0;
DP ym,wt,sum,sig2i;
int ndat =x.size();
Vec_Dp afunc(ma);
Mat_Dp beta(ma,l);
for(j=0;j<ma;j++)
if(ia[j]) mfit++;
if(mfit== 0) nerror("lift: no parameters to be fitted");
for(j=0;j<mfit;j++){ Initialize the (symmetric) matrix,
for(k=0;k<mfit;k++)
covar[j][k]=0.0;
beta[j][0]=0.0;
}
for (i=0;i<ndat;i++){
funcs(x[i],afunc);
ym=y[i];
if(mfit <ma){
for(j=0;j<ma;j++)
if(!ia[j]);
if(!ia[j]) ym -+ a[j]*afunc[j];
}
sig2i=1.0/SQR(sig[i]);
for(j=0;l<ma;l++)
{
if((ia[l})
{
wt=afunc[i]*sig2i;
for(k=0;m=0;m<=1;m++)
if(ia[m]) covar[j][k++] += wt*afunc{m];
beta[j++][0] += ym*wt;
}
}
}
for(j=i;j<mfit;j++)
for(k=0;k<j;k++)
covar[k][j]=covar[j][k];
Mat_DP temp(mfit;j++)
for(j=0;j<mfit;k++)
temp[j][k]=covar[j][k];
gaussj(temp,beta);
for(j=o;j<mfit;j++)
for (k=0;k<mfit;k++)
covar[j][k] =temp[j][k];
for(j=0,l=0;l<ma;l++)

if(ia[l]) a[l]=beta[j++][0];
chis=0.0;
for((i=0;i<nda;i++)
{
funcs(x[i],afunc);
sum=0.0;
for(j=j=0;j<ma;j++) sum += a[j]*afunc[j];
chisq += SQR(y[i]-sum/sig[i]);
}
covsrt(covar,ia,mfit);
}
"void" NR:: covsrt(Mat_IO_DP &covar,Vec_I_Bool &ia,const int mfit)
{
int i,j,k;
int ma=ia.size();
for(i=mfit;i<ma;i++)
for(j=0;j<i+1;j++)
covar[i][j]=covar[j][i]=0.0;
k=mfit-1;
for(j=ma-1;j>=0;j--)
{
if(ia[j]){
fo(i=0;i<ma;i++)SWAP(covar[i][k],covar[i][j]);
for(i=0;i<ma;i++)SWAP(covar[k][i],covar[j][i]);
k--;
}
}
}

int main()
{

return 0;
}

Thanks Kutta for your message. I will work on your solution but it is worth mentioning that I found MINPACK package somewhat easier relative to Numerical Recipes when it comes to optimization of multi-variable functions like the one I have enclosed. However, I always prefer to use NR for my programs and that's why I am asking Dave if he has easier solution. Thanks and warm wishes.

Thanks Kutta for your message. I will work on your solution but it is worth mentioning that I found MINPACK package somewhat easier relative to Numerical Recipes when it comes to optimization of multi-variable functions like the one I have enclosed. However, I always prefer to use NR for my programs and that's why I am asking Dave if he has easier solution. Thanks and warm wishes.
Context to the above and its pre-susequents, an example is appended below
Pl try the following for nonlinear ones of more variables.
Thanking You
As
Kutta(C.R.Muthukumar)

Context to the above and its pre-susequents, an example is appended below
Pl try the following for nonlinear ones of more variables.
Thanking You
As
Kutta(C.R.Muthukumar)

After a long ,a realistic but new C++ is appended below that will
be coordinating what has been eariler said for solving nonlinear ones via linear ones. Indeed this is programmed using Simpsons method as mentioned(PDE) .Hence please ref my old reply and start making the linear ones in sequencial, and then use the formula that has been given for NonLinear Equations which is ur present requirements of the suject at the first instant. For Example only
Thanking You
As
Kutta(C.R.Muthukumar):)

#include <fstream.h>
#include <iomanip.h>
#include <math.h>
#include <vector>

template <class Type>
Type simpson (ofstream& out, Type f(Type), Type a, Type b, Type epsilon, int level, int level_max)
{
Type c, d, e, h;
Type one_simpson,
two_simpson,
simpson_result,
left_simpson,
right_simpson;

level++;
h = b - a;
c = 0.5 * (a + b);

one_simpson = h * (f(a) + 4 * f(c) + f(b)) / 6.0;
d = 0.5 * (a + c);
e = 0.5 * (c + b);
two_simpson = h * (f(a) + 4 * f(d) + 2 * f(c) + 4 * f(e) + f(b)) / 12.0;
if (level >= level_max)
{
simpson_result = two_simpson;
out << "maximum level reached" << endl;
}
else
{
if ((fabs(two_simpson - one_simpson)) < 15.0 * epsilon)
simpson_result = two_simpson + (two_simpson - one_simpson) / 15.0;
else
{
left_simpson = simpson(out, f, a, c, epsilon/2, level, level_max);
right_simpson = simpson(out, f, c, b, epsilon/2, level, level_max);
simpson_result = left_simpson + right_simpson;
}
}
return simpson_result;
}

*************************
Example C++ prg for Example
#include "Simpson.cpp"

typedef double runtype;
typedef std::vector<runtype> Array;
typedef std::vector<Array> Matrix;

inline runtype f(runtype x) {return (1.0/ exp(x * x)); }

//inline runtype f(runtype x) {return (sqrt(4.0-pow(x,2.0)));}

//inline runtype f(runtype x) {return (sin(x));}

//inline runtype f(runtype x) {return (x*x*x+2*(x*x)+3*x+1);}


ofstream outfile("simp4.txt");

int main()
{

int n = 1000;
const runtype a = 0.0;
const runtype b = 1.0;

const runtype epsilon = 0.0000005;
int level = 0;
int maxlevel = 20;
runtype result;

result = simpson(outfile, f, a, b, epsilon, level, maxlevel);
outfile.precision(10);
outfile << "Using Adaptive Simpson's (with " << maxlevel << " runs)" << endl;
outfile << "Approximate integral = " << setw(15) << result << endl;
}

************
Output:
Using Adaptive Simpson's (with 20 runs)
Approximate integral = 0.746824134

Dave,
You provided the following sample program before:
//
// Demonstration of use of Numerical Recipes Fitmrq.
//
// The original function is the sum of two Gaussian functions,
// and the data points are corrupted by noise.
//
// The Numerical Recipes distribution conveniently supplies a
// function that can be fed to Fitmrq to try to model data that
// consists of samples from sum of any number of Gaussian-type functions.
//
// I corrupt the data values with additive Gaussian noise.
//
// davekw7x
//


// Always include nr3.h first
#include "C:\Program Files\Numerical Recipes\NR_C302\code/nr3.h"

// For Normaldev
#include "C:\Program Files\Numerical Recipes\NR_C302\code/ran.h"
#include "C:\Program Files\Numerical Recipes\NR_C302\code/gamma.h"
#include "C:\Program Files\Numerical Recipes\NR_C302\code/deviates.h"

// For Fitmrq
#include "C:\Program Files\Numerical Recipes\NR_C302\code/gaussj.h"
#include "C:\Program Files\Numerical Recipes\NR_C302\code/fitmrq.h"

// The fgauss() function is in the file fit_examples.h
#include "C:\Program Files\Numerical Recipes\NR_C302\code/fit_examples.h"

//
// Simplify appearance of expressions with x*x by using an in-line function
//
inline Doub sqr(const Doub & x)
{
return x * x;
}

int main()
{
int i;

Doub noise = 0.1; // Standard deviation of generated noise
Normaldev ran(0, noise, 8654321); // Gaussian noise generator
//Normaldev ran(0, noise, time(0)); // Gaussian noise generator

Int num_points = 101;

// In order to use fgauss on K Gaussians, we need
// some vectors with size 3*K. So, for two Gaussians
// we have size = 6
//
// The two groups of three parameters contain amplitude, mean,
// and standard deviation of a Gaussian function that is
// used to generate the data points for testing.
//
// This is the actual array of parameters for the two Gaussian
// functions:
// [first group], [second group]
const Doub actual_d[] = { 5.0, 2.0, 3.0, 2.0, 5.0, 3.0 };

//
// This is the "guess" that we supply to fitmrq:
//
// Try it with a "pretty good" guess:
const Doub guess_d[] = { 4.8, 2.2, 2.8, 2.2, 4.8, 3.2 };

//
// Just for kicks, you might want to try it with a "not so pretty
// good" guess. Maybe you will get lucky; maybe not.
// The point is that Levenberg-Marquardt requires an initial
// guess of some kind.
//
//const Doub guess_d[] = { 2.0, 2.0, 2.0, 1.0, 1.0, 1.0 };

//
// The number of elements in the array is needed to initialize
// a VecDoub object from that array.
//
const Int ma = sizeof(actual_d) / sizeof(actual_d[0]);

VecDoub actual(ma, actual_d);
VecDoub guess(ma, guess_d);
VecDoub x(num_points);
VecDoub y(num_points);
VecDoub sd(num_points);

//
// Fill array y with values for sum of two gaussians, corrupted
// by additive gaussian noise.
//
// The sd array contains estimates of the standard deviations
// of the elements of y. We use a "measurement error" factor
// equal to the standard deviation of the noise distribution.
//
Doub xstart = 0.0;
Doub xmax = 10.0;
Doub deltax = xmax / (num_points - 1);
Doub xx;
for (i = 0, xx = xstart; i < num_points; i++, xx += deltax) {
x[i] = xx;
y[i] = ran.dev() + (actual[0]*exp(-sqr((xx - actual[1]) / actual[2])) +
actual[3]*exp(-sqr((xx - actual[4]) / actual[5])));
sd[i] = noise * y[i];
}

//
// Use the constructor to bind the data
//
Fitmrq mymrq(x, y, sd, guess, fgauss); // Use default tolerance 1.0e-3

//
// Now do the fit
//
mymrq.fit();

cout << scientific << setprecision(5);
cout << setw(10) << "Actual"
<< setw(13) << "Fitted"
<< " "
<< setw(13) << "Err. Est."
<< endl;

for (i = 0; i < actual.size(); i++) {
cout << setw(13) << actual[i]
<< setw(13) << mymrq.a[i]
<< " +- " << sqrt(mymrq.covar[i][i])
<< endl;
}
cout << endl;

cout << fixed;
cout << "Number of points = " << num_points
<< ", Chi-squared = " << mymrq.chisq << endl;

return 0;
}

Here is my output with a pretty good guess vector:

Actual Fitted Err. Est.
5.00000e+000 5.52170e+000 +- 6.84019e-001
2.00000e+000 2.14463e+000 +- 4.97134e-001
3.00000e+000 2.93038e+000 +- 4.86329e-001
2.00000e+000 1.78320e+000 +- 1.43392e+000
5.00000e+000 5.46802e+000 +- 9.33546e-001
3.00000e+000 2.56607e+000 +- 2.90910e-001

Number of points = 101, Chi-squared = 309.83746

And here is the output with a bad guess.
const Doub guess_d[] = { 2.0, 2.0, 2.0, 1.0, 1.0, 1.0 }

Actual Fitted Err. Est.
5.00000e+000 3.67964e+000 +- 2.96614e+001
2.00000e+000 1.88054e+000 +- 1.36264e+001
3.00000e+000 3.67881e+000 +- 2.85860e+000
2.00000e+000 2.58572e+000 +- 3.45190e+001
5.00000e+000 3.82958e+000 +- 6.47240e+000
3.00000e+000 3.23738e+000 +- 3.21662e+000

Number of points = 101, Chi-squared = 330.11287

These results are fundamentally different whereas my experience showed that using the MINPACK codes, you'll get consistent close results. I wonder if you could elaborate on this issue of NR3. Thanks

...I wonder if you could elaborate on this issue of NR3...
My post was a response to a request for an example showing how to use some NR3 functions. I did not make, and I do not make any claims as to the suitability of my example for any particular type of problem, including the particular data in that particular example.

The whole point of the Recipes is to expose people to certain kinds of algorithms and to give people code that implements and illustrates the algorithms. Whether any particular piece of code is suitable or appropriate to any particular task is up to the user to decide.


In other words, I showed how to get numbers from a particular kind if data using a particular NR3 function. I don't have any "issue" with NR3.

Maybe some others can share their insight concerning the relative merits of whatever it is that you used or anything else that is suitable and appropriate for this particular type of problem.



Regards,

Dave

Thanks Dave for your message and your helpful guides. You made this forum the best tool for using NR3 book and codes.

Hi again Dave,

Okay. Now, let's take a look at Numerical Recipes Example Book (C++), p.252. It provides the exact same function. If we run the code given in the aforementioned book we obtain the following results with the same initial guess, i.e. const Doub guess_d[] = { 4.8, 2.2, 2.8, 2.2, 4.8, 3.2 };

Iteration # 1 chi-squared: 4342.78 alamda: 0.0001
a[0] a[1] a[2] a[3] a[4] a[5]
4.711811 1.941938 2.973030 2.226970 4.784622 3.082198

Iteration # 2 chi-squared: 163.411174 alamda: 0.000010
a[0] a[1] a[2] a[3] a[4] a[5]
4.831021 1.928210 2.957024 2.200102 4.856893 3.040741

Iteration # 3 chi-squared: 163.411174 alamda: 0.000100
a[0] a[1] a[2] a[3] a[4] a[5]
4.831021 1.928210 2.957024 2.200102 4.856893 3.040741

Iteration # 4 chi-squared: 134.672196 alamda: 0.000010
a[0] a[1] a[2] a[3] a[4] a[5]
4.887688 1.952304 2.972193 2.133448 4.902425 3.028063

Iteration # 5 chi-squared: 134.672196 alamda: 0.000100
a[0] a[1] a[2] a[3] a[4] a[5]
4.887688 1.952304 2.972193 2.133448 4.902425 3.028063

Iteration # 6 chi-squared: 121.991275 alamda: 0.000010
a[0] a[1] a[2] a[3] a[4] a[5]
4.928035 1.968719 2.982424 2.085865 4.936768 3.018292

Iteration # 7 chi-squared: 121.991275 alamda: 0.000100
a[0] a[1] a[2] a[3] a[4] a[5]
4.928035 1.968719 2.982424 2.085865 4.936768 3.018292

Iteration # 8 chi-squared: 115.210931 alamda: 0.000010
a[0] a[1] a[2] a[3] a[4] a[5]
4.957356 1.980871 2.990105 2.050750 4.962292 3.011047

Iteration # 9 chi-squared: 112.005981 alamda: 0.000001
a[0] a[1] a[2] a[3] a[4] a[5]
5.022658 2.008112 3.007475 1.972025 5.019509 2.994861

Iteration # 10 chi-squared: 107.196525 alamda: 0.000000
a[0] a[1] a[2] a[3] a[4] a[5]
5.038977 2.015662 3.012471 1.950552 5.035790 2.990326

Iteration # 11 chi-squared: 107.169242 alamda: 0.000000
a[0] a[1] a[2] a[3] a[4] a[5]
5.039349 2.015868 3.012616 1.949985 5.036243 2.990202

Iteration # 12 chi-squared: 107.169242 alamda: 0.000000
a[0] a[1] a[2] a[3] a[4] a[5]
5.039351 2.015869 3.012616 1.949984 5.036244 2.990201

Iteration # 13 chi-squared: 107.169242 alamda: 0.000000
a[0] a[1] a[2] a[3] a[4] a[5]
5.039351 2.015869 3.012616 1.949984 5.036244 2.990201

Uncertainties:
0.028491 0.012441 0.008296 0.035708 0.026387 0.007459

Expected results:
5.000000 2.000000 3.000000 2.000000 5.000000 3.000000

press return to continue with constraint

holding a[1] and a[4] constant

Iteration # 1 chi-squared: 1799.070951 alamda: 0.000100
a[0] a[1] a[2] a[3] a[4] a[5]
4.995666 2.000000 3.008596 1.993147 5.000000 3.008631

Iteration # 2 chi-squared: 110.712640 alamda: 0.000010
a[0] a[1] a[2] a[3] a[4] a[5]
4.999895 2.000000 3.002986 1.996460 5.000000 3.001130

Iteration # 3 chi-squared: 110.708858 alamda: 0.000001
a[0] a[1] a[2] a[3] a[4] a[5]
4.999913 2.000000 3.002925 1.996542 5.000000 3.001103

Iteration # 4 chi-squared: 110.708858 alamda: 0.000000
a[0] a[1] a[2] a[3] a[4] a[5]
4.999913 2.000000 3.002925 1.996543 5.000000 3.001103

Iteration # 5 chi-squared: 110.708858 alamda: 0.000001
a[0] a[1] a[2] a[3] a[4] a[5]
4.999913 2.000000 3.002925 1.996543 5.000000 3.001103

Uncertainties:
0.001013 0.000000 0.001793 0.001916 0.000000 0.000495

Expected results:
5.000000 2.000000 3.000000 2.000000 5.000000 3.000000

Whereas in the current version we obtain these results:



You see that the fitted coefficient in the old NR codes are so close to the actual values but in current version the results are so off. I wonder if you could help me with this issue.

I think my experience is not so enough in numerical computation but with your helps, I gradually will have more thoughtful procedures. I appreciate it!

...It provides the exact same function...
No it doesn't. Well, it's more-or-less the same functionality, except:

In my example the noise going into the signals is much greater.


I didn't add the noise in the y-values consistent with the the way that I generated the sigma-values, so the results didn't indicate a particularly good fit.



My hope was that people might actually look at the results and look at the code and see where a better treatment could be done.

I mean, I didn't deliberately set it up to give less-than-useful results, but maybe it just turned out that way. The whole point is to get some numbers and then scratch your head and figure out if they are meaningful. You did the "scratch your head" part, but maybe the programs were so similar that you overlooked the differences.


Anyhow, if you want to compare the reasonableness of the results with the other program:

Change

Doub noise = 0.1; // Standard deviation of generated noise
Normaldev ran(0, noise, 8654321); // Gaussian noise generator



to

Doub noise = 0.001; // Standard deviation of generated noise
Normaldev ran(0, 1, 8654321); // Gaussian noise generator




And change

y[i] = ran.dev() + (actual[0]*exp(-sqr((xx - actual[1]) / actual[2])) +
actual[3]*exp(-sqr((xx - actual[4]) / actual[5])));


to


y[i] = (actual[0]*exp(-sqr((xx - actual[1]) / actual[2])) +
actual[3]*exp(-sqr((xx - actual[4]) / actual[5])));
y[i] *= (1.0 + noise*ran.dev());


And try again.

I think the output will be more like the following:

5.00000e+00 4.95142e+00 +- 2.92716e-02
2.00000e+00 1.98081e+00 +- 1.21173e-02
3.00000e+00 2.98750e+00 +- 7.69857e-03
2.00000e+00 2.05812e+00 +- 3.48699e-02
5.00000e+00 4.95587e+00 +- 2.53770e-02
3.00000e+00 3.01329e+00 +- 7.27198e-03

Number of points = 101, Chi-squared = 88.11142



I think my experience is not so enough...

My experience may be more limited than would lead to a really good exposition of the principles, so I wasn't trying to teach how to do a meaningful analysis. I was just trying to show how to get some numbers out of the program.

Maybe I should have specifically emphasized that the chi-square value didn't indicate a particularly good fit and maybe I should have made an explicit challenge to try to improve things in the test program.

Oh, well...

I mean, no one was born knowing this stuff, and I still have a lot to learn too. (A lot!)


I hope this helps.



Regards,

Dave

Thanks Dave. I should be more critical and practice critical thinking.

Dear Dave,

Now the question arises:

What if we want to measure a single variable y as a function of more than one variable, i.e. y=f(x1,x2,...,xn). How can I use firmrq.h in this condition?

For instance consider that y=f(x0,x1) &
we have
VecDoub x0(3), VecDoub x1(5) &
MatDoub y(3,5).

where,
y[0][0]=f(x0[0],x1[0]), y[0][1]=f(x0[0],x1[1]), ..., y[2][1]=f(x0[2],x1[1]), etc.

The goal is to fit some coefficient for a non-linear equation, e.g., a[0], a[1], ..., a[5].

This situation has been vaguely discussed in sec. 15.4.4, p.798 of the NR3 book for Fitsvd but never solved for an existence example.


Would you please tell me how fitmrqmin.h works in this condition?

Hi amirvahid,

do you solve your problem of multidimensional fitting with levenberg-marquadt method with the implementation of the NR 3d edition ? I have the same problem and I don't have any solution.


Dear Dave,

Now the question arises:

What if we want to measure a single variable y as a function of more than one variable, i.e. y=f(x1,x2,...,xn). How can I use firmrq.h in this condition?

For instance consider that y=f(x0,x1) &
we have
VecDoub x0(3), VecDoub x1(5) &
MatDoub y(3,5).

where,
y[0][0]=f(x0[0],x1[0]), y[0][1]=f(x0[0],x1[1]), ..., y[2][1]=f(x0[2],x1[1]), etc.

The goal is to fit some coefficient for a non-linear equation, e.g., a[0], a[1], ..., a[5].

This situation has been vaguely discussed in sec. 15.4.4, p.798 of the NR3 book for Fitsvd but never solved for an existence example.


Would you please tell me how fitmrqmin.h works in this condition?

No. I decided to use MINPACK for this purpose. I found numerical recipes quite awkward for this purpose. It would be really helpful if the authors provide some working examples in this forum. Otherwise, we just have to stick with MINPACK approach (C++ and/or FORTRAN).