forum.alglib.net http://forum.alglib.net/ |
|
Exponential curve b*exp(-a*x) fit problem http://forum.alglib.net/viewtopic.php?f=2&t=14 |
Page 1 of 1 |
Author: | french [ Sat May 29, 2010 10:11 am ] | ||
Post subject: | Exponential curve b*exp(-a*x) fit problem | ||
Hi, I'm trying to fit exponential curve b*exp(-a*x) using LMfit. I adapted source code from example. First, I tried to fit it to set of noised data, with no success. After, I generated ideal exp function, but the result was the same. I was using lsfitnonlinearfgh and also lsfitnonlinearfg. I've changed fit to b*exp(a*x). I tried different combinations, but fitting doesn't work. Fit result is identical as starting values. When I use matlab curve fitting tool with LM algorithm (alglib use LM too), I've got pretty fine results (a = 0.003231, b = 610.1 - fit to data from attached signal.txt in range [350,1000]). What am i doing wrong? I would be grateful for any help. This is my source code: Code: using System; using System.IO; using alglib; public class ExpFit { public static int Main(string[] args) { lsfit.lsfitreport rep = new lsfit.lsfitreport(); lsfit.lsfitstate state = new lsfit.lsfitstate(); int begin = 350; int end = 1000; int n = end - begin; int m = 1; int k = 2; double a = 0.003231; double b = 610.1; // // Prepare task matrix // double[] y = new double[n]; double[,] x = new double[n, m]; double[] c = new double[k]; double[] data = LoadFile("signal.txt"); for(int i = 0; i < n; i++) { x[i,0] = i + begin; y[i] = data[i + begin]; //y[i] = b * Math.Exp(-a * x[i,0]); } double epsf = 0.0; double epsx = 0.0; int maxits = 0; int info = 0; int j = 0; c[0] = 500; c[1] = 0.001; // // Solve // lsfit.lsfitnonlinearfgh(ref x, ref y, ref c, n, m, k, ref state); //lsfit.lsfitnonlinearfg(ref x, ref y, ref c, n, m, k, false, ref state); lsfit.lsfitnonlinearsetcond(ref state, epsf, epsx, maxits); while( lsfit.lsfitnonlineariteration(ref state) ) { // // F(x) = b*exp(-a*x) // state.f = state.c[0] * Math.Exp(-(state.c[1] * state.x[0])); // // F(x) = b*exp(-a*x) // dF/da = -b*x*exp(-a*x) // dF/db = exp(-a*x) // if (state.needfg | state.needfgh) { state.g[0] = -state.c[0] * state.x[0] * Math.Exp(-(state.c[1] * state.x[0])); state.g[1] = Math.Exp(-(state.c[1] * state.x[0])); } // // F(x) = b*exp(-a*x) // d2F/da2 = b*x^2*exp(-a*x) // d2F/dadb = -x*exp(a*x) // d2F/db2 = 0 // if (state.needfgh) { state.h[0, 0] = state.c[0] * AP.Math.Sqr(state.x[0]) * Math.Exp(-(state.c[1] * state.x[0])); state.h[0, 1] = -state.x[0] * Math.Exp(-(state.c[1] * state.x[0])); state.h[1, 0] = -state.x[0] * Math.Exp(-(state.c[1] * state.x[0])); state.h[1, 1] = 0.0; } j++; } lsfit.lsfitnonlinearresults(ref state, ref info, ref c, ref rep); System.Console.Write("a: "); System.Console.Write("{0,0:F8}",c[1]); System.Console.WriteLine(); System.Console.Write("b: "); System.Console.Write("{0,0:F8}",c[0]); System.Console.WriteLine(); System.Console.Write("rms.err: "); System.Console.Write("{0,0:F3}",rep.rmserror); System.Console.WriteLine(); System.Console.Write("max.err: "); System.Console.Write("{0,0:F3}",rep.maxerror); System.Console.WriteLine(); System.Console.Write("iterations: "); System.Console.Write("{0}", j); System.Console.WriteLine(); System.Console.Write("Termination type: "); System.Console.Write("{0,0:d}",info); System.Console.WriteLine(); System.Console.WriteLine(); System.Console.WriteLine(); return 0; } public static double[] LoadFile(string path) { StreamReader sr = new StreamReader(path); uint n = uint.Parse(sr.ReadLine()); uint dummy = uint.Parse(sr.ReadLine()); double[] data = new double[n]; for (uint i = 0; i < n; i++) { data[i] = double.Parse(sr.ReadLine()); } sr.Close(); return data; } } And console output: Code: a: 0.00100000
b: 500.00000000 rms.err: 176.992 max.err: 229.619 iterations: 34450 Termination type: 2
|
Author: | Sergey.Bochkanov [ Sun May 30, 2010 7:20 am ] |
Post subject: | Re: Exponential curve b*exp(-a*x) fit problem |
Code: // // F(x) = b*exp(-a*x) // state.f = state.c[0] * Math.Exp(-(state.c[1] * state.x[0])); // // F(x) = b*exp(-a*x) // dF/da = -b*x*exp(-a*x) // dF/db = exp(-a*x) // if (state.needfg | state.needfgh) { state.g[0] = -state.c[0] * state.x[0] * Math.Exp(-(state.c[1] * state.x[0])); state.g[1] = Math.Exp(-(state.c[1] * state.x[0])); } // // F(x) = b*exp(-a*x) // d2F/da2 = b*x^2*exp(-a*x) // d2F/dadb = -x*exp(a*x) // d2F/db2 = 0 // if (state.needfgh) { state.h[0, 0] = state.c[0] * AP.Math.Sqr(state.x[0]) * Math.Exp(-(state.c[1] * state.x[0])); state.h[0, 1] = -state.x[0] * Math.Exp(-(state.c[1] * state.x[0])); state.h[1, 0] = -state.x[0] * Math.Exp(-(state.c[1] * state.x[0])); state.h[1, 1] = 0.0; } There is mistake in your derivative calculation code. When you calculate F(x,a,b), you assume that c[0]=b, c[1]=a. But when you calculate grad(F) and hess(F), you treat c[0] as a, c[1] as b. So algorithm tries to search in a wrong direction and stops (because it is NOT descent direction). |
Author: | french [ Sun May 30, 2010 11:28 am ] |
Post subject: | Re: Exponential curve b*exp(-a*x) fit problem |
Thank you very much for your reply and pointing mistake. |
Author: | Sergey.Bochkanov [ Wed Nov 28, 2012 11:47 am ] |
Post subject: | Re: Exponential curve b*exp(-a*x) fit problem |
You have an error in your function/gradient calculation code - sometimes you use State.C, and sometimes just C (local array, not fields of State structure). This local array is not updated with recent changes in parameters, so algorithm is confused and tries to make wrong steps. |
Author: | Sergey.Bochkanov [ Wed Nov 28, 2012 11:51 am ] |
Post subject: | Re: Exponential curve b*exp(-a*x) fit problem |
BTW, I recommend you to increase StpMax to 10-20.... Your algorithm has to perform too many small steps in order to reach optimum parameters. |
Page 1 of 1 | All times are UTC |
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group http://www.phpbb.com/ |