Chap15: Modeling of Data

[amirvahid]

Dear all,

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

[kutta]

Originally Posted by amirvahid

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)

[amirvahid]

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.

[davekw7x]

Originally Posted by amirvahid

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.)

[amirvahid]

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

[kutta]

Originally Posted by amirvahid

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)

[kutta]

Originally Posted by amirvahid

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)

[amirvahid]

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

[amirvahid]

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

[kutta]

Originally Posted by amirvahid

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;
}
}

[amirvahid]

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).

[kutta]

Originally Posted by amirvahid

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;
}

[amirvahid]

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.

[kutta]

Originally Posted by amirvahid

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)

[kutta]

Originally Posted by kutta

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