Trying to fit gamma distribution to each separate peak of time series data (chromatography).

As a peak i take local minimum-maximum-minimum part of the data each time. Since the peaks values do not start right from x-axis i add some delta (upshift) to the gamma values (limiting it to peak's smallest adjacent local minimum).

But optim() either fails to converge (with an error "L-BFGS-B needs finite values of 'fn'") or sometimes, if i give it a big start shape value(around 100), converges immediately leaving all start values almost unchanged.

My step in time series is 0.02 second, and peak's widths usually a few hundreds of points, so gamma input values are usually 2-3 seconds, probably it is too short to fit and i have to do some scaling/normalization? Or a problem in something else?

(Voltage is a y axis coordinate, Retention time is an x-axis one)

# Residuals Sum of Squares
RSS.gamma <- function(params, data){

 alpha <- params[[1]] # shape
 theta <- params[[2]] # scale
 delta <- params[[3]] # upshift

 StartNewTime <- seq(from = 0.02, length.out = length(data$RetentionTime), by = 0.02)

  with (data, sum( (Voltage - (delta + dgamma(StartNewTime, shape = alpha, scale = theta, log = FALSE)))^2 ) )

# Rightmost local minimum on the left peak's side
LeftLocalMin = max(SmoothedGraphMinima[SmoothedGraphMinima < SmoothedGraphMaxima[i]])
# Leftmost local minimum on the right peak's side
RightLocalMin = min(SmoothedGraphMinima[SmoothedGraphMinima > SmoothedGraphMaxima[i]])

DeltaMax = min( SmoothedGraph[LeftLocalMin], SmoothedGraph[RightLocalMin] ) 

InitAlpha = 5    # changed it manually in attempt to converge
InitTheta = 1 

PeakSection = Patient[LeftLocalMin:RightLocalMin, c('RetentionTime','Voltage')]  

result <- optim(par = c(InitAlpha, InitTheta, DeltaMax/2 ), RSS.gamma, data = PeakSection, method="L-BFGS-B",
              lower = c(1,0,DeltaMax/2), upper = c(Inf, Inf, DeltaMax),
              control = list(maxit = 30000, trace = 6, REPORT = 10))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.