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