forum.alglib.net http://forum.alglib.net/ |
|
Non-linear curve fit using a split gaussian function http://forum.alglib.net/viewtopic.php?f=2&t=4310 |
Page 1 of 1 |
Author: | NicD [ Fri Dec 06, 2019 8:19 am ] |
Post subject: | Non-linear curve fit using a split gaussian function |
I am trying to implement a non-linear curve fit for asymmetric peaks using a split gaussian function. It uses separate variances for the left and right side of the peak (see https://arxiv.org/pdf/1405.4995.pdf). A symmetric gaussian is straight forward, but I am stuck with implementing the conditional statement of the split function using different variances (c[2] or c[3]) for x<=c[1] and x>c[1], respectively. The code below shows the two sides of the gaussian in the callback, but I don't know how to go about using one for the part of the x array <= c[1], and the other one for x array > c[1]. Any ideas how this can be implemented? Nic Code: #include "stdafx.h"
#include <stdlib.h> #include <stdio.h> #include <math.h> #include "interpolation.h" using namespace alglib; void function_cx_1_func(const real_1d_array &c, const real_1d_array &x, double &func, void *ptr) { // c[0]: scale factor // c[1]: center position mu // c[2]: variance left side (sigma1) // c[3]: variance right side (sigma2) // for x <= c[1]: func = c[0]*sqrt(2.0/M_PI)/(c[2]+c[3])*exp(-0.5*pow(x[0]-c[1],2.0)/pow(c[2],2.0)); // for x > c[1]: func = c[0]*sqrt(2.0/M_PI)/(c[2]+c[3])*exp(-0.5*pow(x[0]-c[1],2.0)/pow(c[3],2.0)); } int main(int argc, char **argv) { real_2d_array x = "[[0.0],[0.1],[0.2],[0.3],[0.4],[0.5],[0.6],[0.7],[0.8],[0.9],[1.0],[1.1],[1.2],[1.3],[1.4],[1.5]]"; real_1d_array y = "[0.00, 0.00, 0.01, 0.05, 0.12, 0.31, 0.63, 1.10, 1.65, 2.10, 2.28, 1.38, 0.31, 0.03, 0.00, 0.00]"; real_1d_array c = "[1.0, 1.0, 0.25, 0.10]"; double epsx = 0.000001; ae_int_t maxits = 0; ae_int_t info; lsfitstate state; lsfitreport rep; double diffstep = 0.0001; lsfitcreatef(x, y, c, diffstep, state); lsfitsetcond(state, epsx, maxits); alglib::lsfitfit(state, function_cx_1_func); lsfitresults(state, info, c, rep); } |
Page 1 of 1 | All times are UTC |
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group http://www.phpbb.com/ |