fitting a gamma variate curve to a set of data points in c++

Question

I have an array of values (concentration values), with each value taken at a different time point. I need to fit a gamma-variate curve (formula is in the picture below) to these values (i.e. find alpha and beta such that the curve best fits those points - all other variables are known.)

an example of the values i might get (crosses), and the curve I'd want to fit:

I have no idea how to do this. I tried to fit a simplified version of the formula, one that can be solved by using linear regression, by using matrices but I couldn't get it to work. That version of the formula (in which you only solve for one variable, alpha) looks like this:

simplified version, which would also be fine:

my attempt to solve fit the linear regression curve using matrices, using the vnl library (https://vxl.github.io/doc/release/core/vnl/html/index.html) looked like this. I was following this guy's tutorial (https://machinelearningmastery.com/solve-linear-regression-using-linear-algebra/)

  //this is the "data", m_Timing is a vector containing the time each of the data were taken at. 
  vectorVoxel = inputVectorVolumeIter.Get();

  // linear regression 
  
  //declaring the independent (y) and dependent values (x) for our linear regression 
  vnl_matrix<double> ind_var(timeSize, 1);  
  vnl_vector<double> dep_var(timeSize); 
 
  //this vector will hold the "weights" of the fitted line - although in this case there should only be 1 weight 
  vnl_vector<double> weights(1); 

  //loading the values into our independent and dependent variable holders 
  for (index = 0; index < timeSize; index++) {
    ind_var(index, 0) =  (1 + log(m_Timing[index]/tempTTP) - (m_Timing[index]/tempTTP)); 
    dep_var.put(index, log(vectorVoxel[index]));
 }
  
  //fitting the curve! 
  weights = (ind_var.transpose() * ind_var) * ind_var.transpose() * dep_var; 
 

It doesn't work - the weights vector, which should contain the coefficient (alpha) of the fitted line, just contains "null".

The code I'm working on uses the itk library (a library for medical image processsing), and I'm also using vnl for matrices. Is there any way to do this?

Thank you for reading this! I really appreciate any help/even just pointing me in the right direction because I have no idea how to proceed.

Answer 1

This is a problem which is not best suitable to solving by ITK. While you could use ITK's Optimizer infrastructure, there are better/simpler choices.

Maybe try NLOpt? Here is an example of how to use it. Also, you could look at this code which fits a polynomial to points in 3D space.