Good day all!
I use C# version of ALGLIB 3.13.0 (source code generated 2017-12-29).
My experimental data are approximated by the sum of the Gaussians (1 - 50) and the background line.
For approximation I use lsfitcreatef method, which works nice for me at least for ~ 30 Gaussians in data. Example of aproximation for 1 Gaussian + line is in the attachement.
To increase the speed and reliability of convergence, I tried to use lsfitcreatefg method but failed.
For clarity, here is my code:
Code:
static void Main(string[] args)
{
string filepath = @"D:\data.txt";
StreamReader reader = new StreamReader(filepath);
double[] spectr = new double[45];
for (int counter = 0; counter < spectr.Length; counter++)
{
spectr[counter] = Double.Parse(reader.ReadLine()); // y data
}
reader.Close();
reader.Dispose();
int start_channel = 2107;
double[,] x = new double[spectr.Length, 1];
for (int i = 0; i < x.Length; i++)
{
x[i, 0] = (double)(i + start_channel); // x data
}
double[] s = new double[] { 10000, 1000, 1000, 1 }; // scale
double[] c = new double[] { 18000, 2127, 3747, -3.71 }; // starting values
double[] bndl = new double[] { 0, 0, Double.NegativeInfinity, Double.NegativeInfinity };
double[] bndu = new double[] { +20000, +2150, Double.PositiveInfinity, Double.PositiveInfinity };
double epsx = 0.0000001;
int maxits = 0;
int info;
alglib.lsfitstate state;
alglib.lsfitreport rep;
double diffstep = 0.000001;
alglib.lsfitcreatef(x, spectr, c, diffstep, out state);
alglib.lsfitsetbc(state, bndl, bndu);
alglib.lsfitsetcond(state, epsx, maxits);
alglib.lsfitsetscale(state, s);
alglib.lsfitfit(state, function_func, null, null);
alglib.lsfitresults(state, out info, out c, out rep);
System.Console.WriteLine("{0}", info);
System.Console.WriteLine("{0}", alglib.ap.format(c, 6));
/* alglib.lsfitcreatefg(x, spectr, c, true, out state);
alglib.lsfitsetbc(state, bndl, bndu);
alglib.lsfitsetcond(state, epsx, maxits);
alglib.lsfitsetscale(state, s);
alglib.lsfitfit(state, function_func, function_grad, null, null);
alglib.lsfitresults(state, out info, out c, out rep);
System.Console.WriteLine("{0}", info);
System.Console.WriteLine("{0}", alglib.ap.format(c, 6)); */
}
public static void function_func(double[] c, double[] x, ref double func, object obj)
{
func = c[0] * 1 / (Math.Pow(2 * Math.PI, 0.5) * Math.Pow(1.72, 0.5)) * Math.Exp(-(c[1] - x[0]) * (c[1] - x[0]) / (2 * 1.72 * 1.72)) + c[2] + c[3] * x[0] ;
}
public static void function_grad(double[] c, double[] x, ref double func, double[] grad, object obj)
{
grad[0] = 1 / (Math.Pow(2 * Math.PI, 0.5) * Math.Pow(1.72, 0.5)) * Math.Exp(-(c[1] - x[0]) * (c[1] - x[0]) / (2 * 1.72 * 1.72)) + c[2] + c[3] * x[0];
grad[1] = (c[0] * 1 / (Math.Pow(2 * Math.PI, 0.5) * Math.Pow(1.72, 0.5)) * Math.Exp(-(c[1] - x[0]) * (c[1] - x[0]) / (2 * 1.72 * 1.72)) + c[2] + c[3] * x[0]) *
((x[0] - c[1]) / (1.72 * 1.72));
grad[2] = 1;
grad[3] = x[0];
}
{
string filepath = @"D:\data.txt";
StreamReader reader = new StreamReader(filepath);
double[] spectr = new double[45];
for (int counter = 0; counter < spectr.Length; counter++)
{
spectr[counter] = Double.Parse(reader.ReadLine()); // y data
}
reader.Close();
reader.Dispose();
int start_channel = 2107;
double[,] x = new double[spectr.Length, 1];
for (int i = 0; i < x.Length; i++)
{
x[i, 0] = (double)(i + start_channel); // x data
}
double[] s = new double[] { 10000, 1000, 1000, 1 }; // scale
double[] c = new double[] { 18000, 2127, 3747, -3.71 }; // starting values
double[] bndl = new double[] { 0, 0, Double.NegativeInfinity, Double.NegativeInfinity };
double[] bndu = new double[] { +20000, +2150, Double.PositiveInfinity, Double.PositiveInfinity };
double epsx = 0.0000001;
int maxits = 0;
int info;
alglib.lsfitstate state;
alglib.lsfitreport rep;
double diffstep = 0.000001;
alglib.lsfitcreatef(x, spectr, c, diffstep, out state);
alglib.lsfitsetbc(state, bndl, bndu);
alglib.lsfitsetcond(state, epsx, maxits);
alglib.lsfitsetscale(state, s);
alglib.lsfitfit(state, function_func, null, null);
alglib.lsfitresults(state, out info, out c, out rep);
System.Console.WriteLine("{0}", info);
System.Console.WriteLine("{0}", alglib.ap.format(c, 6));
/* alglib.lsfitcreatefg(x, spectr, c, true, out state);
alglib.lsfitsetbc(state, bndl, bndu);
alglib.lsfitsetcond(state, epsx, maxits);
alglib.lsfitsetscale(state, s);
alglib.lsfitfit(state, function_func, function_grad, null, null);
alglib.lsfitresults(state, out info, out c, out rep);
System.Console.WriteLine("{0}", info);
System.Console.WriteLine("{0}", alglib.ap.format(c, 6)); */
}
public static void function_func(double[] c, double[] x, ref double func, object obj)
{
func = c[0] * 1 / (Math.Pow(2 * Math.PI, 0.5) * Math.Pow(1.72, 0.5)) * Math.Exp(-(c[1] - x[0]) * (c[1] - x[0]) / (2 * 1.72 * 1.72)) + c[2] + c[3] * x[0] ;
}
public static void function_grad(double[] c, double[] x, ref double func, double[] grad, object obj)
{
grad[0] = 1 / (Math.Pow(2 * Math.PI, 0.5) * Math.Pow(1.72, 0.5)) * Math.Exp(-(c[1] - x[0]) * (c[1] - x[0]) / (2 * 1.72 * 1.72)) + c[2] + c[3] * x[0];
grad[1] = (c[0] * 1 / (Math.Pow(2 * Math.PI, 0.5) * Math.Pow(1.72, 0.5)) * Math.Exp(-(c[1] - x[0]) * (c[1] - x[0]) / (2 * 1.72 * 1.72)) + c[2] + c[3] * x[0]) *
((x[0] - c[1]) / (1.72 * 1.72));
grad[2] = 1;
grad[3] = x[0];
}
For lsfitcreatef method all is fine. Parameters values change from c=[18000, 2127, 3747, -3.71] to c = [10564, 2127, 624, -0.248]. Very good result. Picture with this result is in the attachement.
If I remove the commenting for lsfitcreatefg method in the my code and comment lsfitcreatef method instead the result is c=[17999.99, 2126.99, 3747, -3.71].
The resulting values of the parameters remained practically unchanged.
I checked the accuracy of the computation of gradient functions - all right.
Could you explain what do I do wrong?
Regards
| Attachments: |
|
Exp_vs_Approx.PNG [ 38.07 KiB | Viewed 5791 times ] |