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.)

gamma-variate curve formula

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.

I need your data (numerical, not a graph) in order to check a particular method of regression and see if this method is convenient in your case. — JJacquelin, May 09 '21 at 15:15
@JJacquelin thank you for replying! the data will change every time the program is used (the program generates does mathematical processing on brain CT scans repeatedly taken over about a minute, as the patient is injected with contrast - each group of data is actually (the density measurement of) a single "pixel" from a person's brain taken at different time points over that minute. Here is example data (transformed so the regression can be done - I am not sure this is correctly done on my part): — user123, May 10 '21 at 12:59
@JJacquelin dependent variable (log of the density at those times): [3.55, 2.77, 3.13, 3.25, 3.465, 2.89, 3.58, 2.39, 3.46, 3.13, 3.61, 2.9, 2.19, 3.29] independent variable (derived from time values): 0.0, -2.1, -1.5, -1.16, 0.0, -0.0035, -0.007, -0.021, -0.038, -0.059, -0.08, -0.1, -0.14 — user123, May 10 '21 at 13:01
Cannot use it because they are 14 numbers in the first file but only 13 in the second file. — JJacquelin, May 10 '21 at 15:06
You wrote "It doesn't work" at end of your question. I would like to use the same data as yourself in the case where it doesn't work in order to reproduce the trouble that you faced and see where exactly is the difficulty. — JJacquelin, May 10 '21 at 15:13

score 0 · Accepted Answer · answered May 09 '21 at 17:14

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.

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

1 Answers1