forum.alglib.net :: View topic - Nonlinear least squares fitting

private void button2_Click(object sender, EventArgs e)
{
double[,] x = new double[,] { { 0.0000 }, { 0.4900 }, { 1.0000 }, { 1.4700 }, { 1.9600 }, { 2.4500 }, { 2.9400 }, { 3.4300 }, { 3.9200 }, { 4.4100 }, { 4.9000 } };
double[] y = new double[] { 0.0000, 0.1003, 0.1375, 0.1581, 0.1743, 0.1856, 0.1952, 0.2027, 0.2096, 0.2147, 0.2213, };
double[] c = new double[] { 0.2, 0.5, 0.7 };
double epsf = 0;
double epsx = 0.000001;
int maxits = 0;
int info;
alglib.lsfitstate state;
alglib.lsfitreport rep;
double diffstep = 0.0001;
// Fitting without weights
// X - array[0..N-1,0..M-1], points (one row = one point)
// Y - array[0..N-1], function values.
// C - array[0..K-1], initial approximation to the solution,
// N - number of points, N>1
// M - dimension of space
// K - number of parameters being fitted
// DiffStep- numerical differentiation step; should not be very small or large; large = loss of accuracy small = growth of round-off errors
alglib.lsfitcreatef(x, y, c, diffstep, out state);
alglib.lsfitsetcond(state, epsf, epsx, maxits);
//alglib.lsfitfit(state, function_cx_1_func, null, null); // using function values only
//alglib.lsfitfit(state, function_cx_1_func, function_cx_1_grad, null, null); // Gradient
alglib.lsfitfit(state, function_cx_1_func, function_cx_1_grad, function_cx_1_hess, null, null); // Gradient + Hessian
alglib.lsfitresults(state, out info, out c, out rep);

label2.Text = alglib.ap.format(c, 7);
label3.Text = info.ToString();
}

public static void function_cx_1_func(double[] c, double[] x, ref double func, object obj)
{
////c[0] == s
////c[1] == a
////c[2] == n
////y = (s (1-exp(-a x^n)))
//func = (1 - Math.Exp(-c[0] * Math.Pow(x[0], c[1]))); //(1-exp(-a x^n))
func = c[0] * (1 - Math.Exp(-c[1] * Math.Pow(x[0], c[2]))); // s*(1-exp(-a*pow(x,n)))
}

public static void function_cx_1_grad(double[] c, double[] x, ref double func, double[] grad, object obj)
{
////c[0] == s
////c[1] == a
////c[2] == n
////(d)/(ds) (s (1-exp(-a x^n))) = 1-e^(-a x^n)
////(d)/(da) (s (1-exp(-a x^n))) = s x^n e^(-a x^n)
////(d)/(dn) (s (1-exp(-a x^n))) = a s x^n log(x) e^(-a x^n)
grad[0] = 1 - Math.Exp(-c[1] * Math.Pow(x[0],c[2]));//1-e^(-a x^n)
grad[1] = c[0] * Math.Pow(x[0],c[2]) * Math.Exp(-c[1] * Math.Pow(x[0],c[2])); //s x^n e^(-a x^n)
grad[2] = c[1] * c[0] * Math.Pow(x[0], c[2]) * Math.Log10(x[0]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2])); //a s x^n log(x) e^(-a x^n)

}

public static void function_cx_1_hess(double[] c, double[] x, ref double func, double[] grad, double[,] hess, object obj)
{
////c[0] == s
////c[1] == a
////c[2] == n
////(d^2)/(ds^2) (s (1-exp(-a x^n))) = 0
////(d^2)/(da^2) (s (1-exp(-a x^n))) = s x^(2 n) (-e^(-a x^n))
////(d^2)/(dn^2) (s (1-exp(-a x^n))) = -a s x^n log^2(x) e^(-a x^n) (a x^n-1)
func = c[0] * (1 - Math.Exp(-c[1] * Math.Pow(x[0], c[2])));
grad[0] = 1 - Math.Exp(-c[1] * Math.Pow(x[0], c[2]));//1-e^(-a x^n)
grad[1] = c[0] * Math.Pow(x[0], c[2]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2])); //s x^n e^(-a x^n)
grad[2] = c[1] * c[0] * Math.Pow(x[0], c[2]) * Math.Log10(x[0]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2])); //a s x^n log(x) e^(-a x^n)
// [c0,c0]...[c0,c2]
//(d^2)/(ds^2) (s (1-exp(-a x^n))) = 0
//(d)/(ds)(d)/(da) (s (1-exp(-a x^n))) = x^n e^(-a x^n)
//(d)/(ds)(d)/(dn) (s (1-exp(-a x^n))) = a x^n log(x) e^(-a x^n)
hess[0,0] = 0;
hess[0, 1] = Math.Pow(x[0], c[2]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2]));
hess[0, 2] = c[1] * Math.Pow(x[0], c[2]) * Math.Log10(x[0]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2]));
//[c1,c0]...[c1,c2]
//(d)/(da)(d)/(ds) (s (1-exp(-a x^n))) = x^n e^(-a x^n)
//(d^2)/(da^2) (s (1-exp(-a x^n))) = s x^(2 n) (-e^(-a x^n))
//(d)/(da)(d)/(dn) (s (1-exp(-a x^n))) = s x^n log(x) (-e^(-a x^n)) (a x^n-1)
hess[1, 0] = Math.Pow(x[0], c[2]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2]));
hess[1, 1] = c[0] * Math.Pow(x[0],(2 * c[2])) * (-Math.Exp(-c[1] * Math.Pow(x[0], c[2])));
hess[1, 2] = c[0] * Math.Pow(x[0], c[2]) * Math.Log10(x[0]) * (-Math.Exp(-c[1] * Math.Pow(x[0], c[2]))) * (c[1] * Math.Pow(x[0],(c[2] - 1)));
//[c2,c0]...[c2,c2]
//(d)/(dn)(d)/(ds) (s (1-exp(-a x^n))) = a x^n log(x) e^(-a x^n)
//(d)/(dn)(d)/(da) (s (1-exp(-a x^n))) = s x^n log(x) (-e^(-a x^n)) (a x^n-1)
//(d^2)/(dn^2) (s (1-exp(-a x^n))) = -a s x^n log^2(x) e^(-a x^n) (a x^n-1)
hess[2, 0] = c[1] * Math.Pow(x[0], c[2]) * Math.Log10(x[0]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2]));
hess[2, 1] = c[0] * Math.Pow(x[0], c[2]) * Math.Log10(x[0]) * (-Math.Exp(-c[1] * Math.Pow(x[0], c[2]))) * (c[1] * Math.Pow(x[0], (c[2] - 1)));
hess[2, 2] = -c[1] * c[0] * Math.Pow(x[0], c[2]) * Math.Pow(Math.Log10(x[0]),2) * (-Math.Exp(-c[1] * Math.Pow(x[0], c[2]))) * (c[1] * Math.Pow(x[0], (c[2] - 1)));
}

Author:	stl [ Sat Apr 02, 2016 12:26 pm ]
Post subject:	Nonlinear least squares fitting
Hello, i am implementing curve fitting in C# and i tried curve fitting by function values, gradient and i am trying Hessian right now. But i am not really sure about Hessian values, I am not so good in math and this is the first time for me implementing some huge math problem. My function is: y = s(1-exp(-apow(x,n))) where s, a, n are unknown variables which i am looking for. Could someone check my Hessian, please? I hope it is correct, but i need to be sure. Also i have question about tweaking the calculation, by which parameters i can influence precision of the results? Thank you very much! Code: private void button2_Click(object sender, EventArgs e) { double[,] x = new double[,] { { 0.0000 }, { 0.4900 }, { 1.0000 }, { 1.4700 }, { 1.9600 }, { 2.4500 }, { 2.9400 }, { 3.4300 }, { 3.9200 }, { 4.4100 }, { 4.9000 } }; double[] y = new double[] { 0.0000, 0.1003, 0.1375, 0.1581, 0.1743, 0.1856, 0.1952, 0.2027, 0.2096, 0.2147, 0.2213, }; double[] c = new double[] { 0.2, 0.5, 0.7 }; double epsf = 0; double epsx = 0.000001; int maxits = 0; int info; alglib.lsfitstate state; alglib.lsfitreport rep; double diffstep = 0.0001; // Fitting without weights // X - array[0..N-1,0..M-1], points (one row = one point) // Y - array[0..N-1], function values. // C - array[0..K-1], initial approximation to the solution, // N - number of points, N>1 // M - dimension of space // K - number of parameters being fitted // DiffStep- numerical differentiation step; should not be very small or large; large = loss of accuracy small = growth of round-off errors alglib.lsfitcreatef(x, y, c, diffstep, out state); alglib.lsfitsetcond(state, epsf, epsx, maxits); //alglib.lsfitfit(state, function_cx_1_func, null, null); // using function values only //alglib.lsfitfit(state, function_cx_1_func, function_cx_1_grad, null, null); // Gradient alglib.lsfitfit(state, function_cx_1_func, function_cx_1_grad, function_cx_1_hess, null, null); // Gradient + Hessian alglib.lsfitresults(state, out info, out c, out rep); label2.Text = alglib.ap.format(c, 7); label3.Text = info.ToString(); } public static void function_cx_1_func(double[] c, double[] x, ref double func, object obj) { ////c[0] == s ////c[1] == a ////c[2] == n ////y = (s (1-exp(-a x^n))) //func = (1 - Math.Exp(-c[0] * Math.Pow(x[0], c[1]))); //(1-exp(-a x^n)) func = c[0] * (1 - Math.Exp(-c[1] * Math.Pow(x[0], c[2]))); // s(1-exp(-apow(x,n))) } public static void function_cx_1_grad(double[] c, double[] x, ref double func, double[] grad, object obj) { ////c[0] == s ////c[1] == a ////c[2] == n ////(d)/(ds) (s (1-exp(-a x^n))) = 1-e^(-a x^n) ////(d)/(da) (s (1-exp(-a x^n))) = s x^n e^(-a x^n) ////(d)/(dn) (s (1-exp(-a x^n))) = a s x^n log(x) e^(-a x^n) grad[0] = 1 - Math.Exp(-c[1] * Math.Pow(x[0],c[2]));//1-e^(-a x^n) grad[1] = c[0] * Math.Pow(x[0],c[2]) * Math.Exp(-c[1] * Math.Pow(x[0],c[2])); //s x^n e^(-a x^n) grad[2] = c[1] * c[0] * Math.Pow(x[0], c[2]) * Math.Log10(x[0]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2])); //a s x^n log(x) e^(-a x^n) } public static void function_cx_1_hess(double[] c, double[] x, ref double func, double[] grad, double[,] hess, object obj) { ////c[0] == s ////c[1] == a ////c[2] == n ////(d^2)/(ds^2) (s (1-exp(-a x^n))) = 0 ////(d^2)/(da^2) (s (1-exp(-a x^n))) = s x^(2 n) (-e^(-a x^n)) ////(d^2)/(dn^2) (s (1-exp(-a x^n))) = -a s x^n log^2(x) e^(-a x^n) (a x^n-1) func = c[0] * (1 - Math.Exp(-c[1] * Math.Pow(x[0], c[2]))); grad[0] = 1 - Math.Exp(-c[1] * Math.Pow(x[0], c[2]));//1-e^(-a x^n) grad[1] = c[0] * Math.Pow(x[0], c[2]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2])); //s x^n e^(-a x^n) grad[2] = c[1] * c[0] * Math.Pow(x[0], c[2]) * Math.Log10(x[0]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2])); //a s x^n log(x) e^(-a x^n) // [c0,c0]...[c0,c2] //(d^2)/(ds^2) (s (1-exp(-a x^n))) = 0 //(d)/(ds)(d)/(da) (s (1-exp(-a x^n))) = x^n e^(-a x^n) //(d)/(ds)(d)/(dn) (s (1-exp(-a x^n))) = a x^n log(x) e^(-a x^n) hess[0,0] = 0; hess[0, 1] = Math.Pow(x[0], c[2]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2])); hess[0, 2] = c[1] * Math.Pow(x[0], c[2]) * Math.Log10(x[0]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2])); //[c1,c0]...[c1,c2] //(d)/(da)(d)/(ds) (s (1-exp(-a x^n))) = x^n e^(-a x^n) //(d^2)/(da^2) (s (1-exp(-a x^n))) = s x^(2 n) (-e^(-a x^n)) //(d)/(da)(d)/(dn) (s (1-exp(-a x^n))) = s x^n log(x) (-e^(-a x^n)) (a x^n-1) hess[1, 0] = Math.Pow(x[0], c[2]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2])); hess[1, 1] = c[0] * Math.Pow(x[0],(2 * c[2])) * (-Math.Exp(-c[1] * Math.Pow(x[0], c[2]))); hess[1, 2] = c[0] * Math.Pow(x[0], c[2]) * Math.Log10(x[0]) * (-Math.Exp(-c[1] * Math.Pow(x[0], c[2]))) * (c[1] * Math.Pow(x[0],(c[2] - 1))); //[c2,c0]...[c2,c2] //(d)/(dn)(d)/(ds) (s (1-exp(-a x^n))) = a x^n log(x) e^(-a x^n) //(d)/(dn)(d)/(da) (s (1-exp(-a x^n))) = s x^n log(x) (-e^(-a x^n)) (a x^n-1) //(d^2)/(dn^2) (s (1-exp(-a x^n))) = -a s x^n log^2(x) e^(-a x^n) (a x^n-1) hess[2, 0] = c[1] * Math.Pow(x[0], c[2]) * Math.Log10(x[0]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2])); hess[2, 1] = c[0] * Math.Pow(x[0], c[2]) * Math.Log10(x[0]) * (-Math.Exp(-c[1] * Math.Pow(x[0], c[2]))) * (c[1] * Math.Pow(x[0], (c[2] - 1))); hess[2, 2] = -c[1] * c[0] * Math.Pow(x[0], c[2]) * Math.Pow(Math.Log10(x[0]),2) * (-Math.Exp(-c[1] * Math.Pow(x[0], c[2]))) * (c[1] * Math.Pow(x[0], (c[2] - 1))); } Thank you for this library, I found it as the easiest to understand and implement unlike other ones. So as a student i am considering small donation to this project (via paypal if it is possible).

Author:	Sergey.Bochkanov [ Mon Apr 04, 2016 1:09 pm ]
Post subject:	Re: Nonlinear least squares fitting
Hello! First, your code does NOT perform Hessian-based optimization :) you created optimizer object with lsfitcreatef function call, so it works in "function value + numerical differentiation" mode. Even gradient values are not used. You should call some other constructor function (lsfitcreatefg for gradient-based optimization, lsfitcreatefgh for Hessian-based optimization). Second, it is hard to check your Hessian analytically - it is really huge. However, you can easily perform numerical check. Just calculate Hessian matrix using your code at some arbitrary point (use non-trivial values for coefficients). And then calculate numerical Hessian (use one of numerical differentiation formulae, some of the basic ones), and compare results. If they are dissimilar, then you have an error. If they are similar... well, either everything is correct - or you made "even number of mistakes" - choose a point where incorrect formula gives correct results purely by chance alone. As for the donation... I really appreciate your intentions! But I think the best donation would be to recommend library to someone else in a search for numerical software :) do it next time you see someone trying to solve some numerical problem, and everyone will be happy :)

Author:	stl [ Mon Apr 04, 2016 6:02 pm ]
Post subject:	Re: Nonlinear least squares fitting
Thank you for help, i forgot about other constructors! So i modified code and tried to run gradient and hessian fitting and i am getting back same array as initial "approximation to the solution" (c). I found that there was no iterations in algorithm (rep.iterationscount) and as info is "4 gradient norm is no more than EpsG". I checked gradients and they should be correct, i have also searched in examples and manual, but i have no clue how to solve this. What am i missing? Thanks.

Author:	Sergey.Bochkanov [ Tue Apr 05, 2016 8:38 am ]
Post subject:	Re: Nonlinear least squares fitting
I just noticed that your gradient calculation code does not set func parameter. It should calculate function value + gradient, but it calculates just gradient. It confuses optimizer. Try fixing it. Does it solve your problem?

Author:	stl [ Tue Apr 05, 2016 9:36 am ]
Post subject:	Re: Nonlinear least squares fitting
Yes, i modified gradient calculation while i was adding new constructors for gradient and hessian. I have tried different values for epsf, epsx, maxits and change CheapFG to true or false. I am still getting same results and same info state (4). Actual code (hessian calculation is not included) is: Code: private double[] LeastSquaresFitting() { double[,] x = new double[,] { { 0.0000 }, { 0.4900 }, { 1.0000 }, { 1.4700 }, { 1.9600 }, { 2.4500 }, { 2.9400 }, { 3.4300 }, { 3.9200 }, { 4.4100 }, { 4.9000 } }; double[] y = new double[] { 0.0000, 0.1003, 0.1375, 0.1581, 0.1743, 0.1856, 0.1952, 0.2027, 0.2096, 0.2147, 0.2213, }; double[] c = new double[] { 0.2, 0.5, 0.7 }; double epsf = 0; double epsx = 0.000001; //0.000001; int maxits = 0; int info; alglib.lsfitstate state; alglib.lsfitreport rep; double diffstep = 0.0001; c = new double[] { 0.2, 0.5, 0.7 }; int Mode = 2; //1- function values, 2- gradient, 3-hessian switch (Mode) { case 1: alglib.lsfitcreatef(x, y, c, diffstep, out state); alglib.lsfitfit(state, function_cx_1_func, null, null); break; case 2: alglib.lsfitcreatefg(x, y, c, true, out state); alglib.lsfitfit(state, function_cx_1_func, function_cx_1_grad, null, null); break; case 3: alglib.lsfitcreatefgh(x, y, c, out state); alglib.lsfitfit(state, function_cx_1_func, function_cx_1_grad, function_cx_1_hess, null, null); break; default: alglib.lsfitcreatef(x, y, c, diffstep, out state); alglib.lsfitfit(state, function_cx_1_func, null, null); break; } alglib.lsfitsetcond(state, epsf, epsx, maxits); alglib.lsfitresults(state, out info, out c, out rep); LogEvent("Linearization ended: " + rep.iterationscount.ToString() + " iterations; status " + info.ToString() + ";"); return c; } public static void function_cx_1_func(double[] c, double[] x, ref double func, object obj) { ////c[0] == s c[1] == a c[2] == n y = (s (1-exp(-a x^n))) func = c[0] * (1 - Math.Exp(-c[1] * Math.Pow(x[0], c[2]))); // s(1-exp(-apow(x,n))) } public static void function_cx_1_grad(double[] c, double[] x, ref double func, double[] grad, object obj) { ////c[0] == s c[1] == a c[2] == n ////(d)/(ds) (s (1-exp(-a x^n))) = 1-e^(-a x^n) ////(d)/(da) (s (1-exp(-a x^n))) = s x^n e^(-a x^n) ////(d)/(dn) (s (1-exp(-a x^n))) = a s x^n log(x) e^(-a x^n) func = c[0] * (1 - Math.Exp(-c[1] * Math.Pow(x[0], c[2]))); // s(1-exp(-apow(x,n))) grad[0] = 1 - Math.Exp(-c[1] * Math.Pow(x[0],c[2])); //1-e^(-a x^n) grad[1] = c[0] * Math.Pow(x[0],c[2]) * Math.Exp(-c[1] * Math.Pow(x[0],c[2])); //s x^n e^(-a x^n) grad[2] = c[1] * c[0] * Math.Pow(x[0], c[2]) * Math.Log10(x[0]) * Math.Exp(-c[1] * Math.Pow(x[0], c[2])); //a s x^n log(x) e^(-a x^n) } Thank you very much. :)

forum.alglib.net http://forum.alglib.net/

Nonlinear least squares fitting http://forum.alglib.net/viewtopic.php?f=2&t=3745	Page 1 of 1

Author:	stl [ Tue Apr 05, 2016 12:19 pm ]
Post subject:	Re: Nonlinear least squares fitting
Thank you for pointing logarithm error. I have corrected that and set upper and lower bounds as constraints. I have tried Lower bound from 0.01 to 0.2 and upper from 1.0 to 2.0 but still getting same results and info (4) everytime.

Author:	Sergey.Bochkanov [ Wed Apr 06, 2016 9:11 am ]
Post subject:	Re: Nonlinear least squares fitting
Hello! Here are my findings: 1. if you want to experiment with setcond(), you should insert it BEFORE call to lsfitfit(). There is no sense in setting stopping conditions after optimization session is done :) 2. optimization algorithm stalls because one of x[] has exactly zero value. However, derivative for x^n implicitly assumes that x is non-zero, so for zero x it results in attempt to calculate Ln(0). Which, in turn, results in NaN value, which confuses optimizer. You should either (a) replace zero x by small non-zero value, or (b) fix your derivative formula to handle borderline cases correctly. 3. I do not know why you got no success with finite differentiation version of algorithm. In my experiments it worked correctly because it avoids triggering this Log-related bug. In any case, it is a good option for fast prototyping. I recommend you to use as much as possible.

Author:	stl [ Wed Apr 06, 2016 7:50 pm ]
Post subject:	Re: Nonlinear least squares fitting
Hi Sergey, thank you very much for your time! When i was adding corect constructors i also moved setcond() before lsfitfit(). I have replaced zero value by 0.0001 and it works. So right now i am able to calculate function parameters by all three linearization modes. Thank you very much for your support.

Page 1 of 1	All times are UTC
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group http://www.phpbb.com/