Linear curve fitting to data in a 3D matrix

Share icon

Requesting urgent help: 

 

I have two 3D matrices. Dimensions are 23 x 19 x 5. One has mass concentrations. Second has corresponding peak areas. Essentially I have 23 x 19 locations for which I need to plot five-point calibration curves and calculate linear response factors i.e. slopes for mass versus peak area plots. 

I've written the code but for some reason, it is not storing standard deviation values from the slope. Instead it just assigns it 0 for almost all 23 x 19. In some places it assigns values but they are completely random. The code is able to calculate slope for each of the 23 x 19 points but the values are incorrect. I'm not sure what mistake have I made. The slope and its standard deviation are supposed to be pretty small i.e. on the order of 1e-5 so I'm wondering if I'm not assigning a variable properly which is why the code is calculating incorrectly. 

Here goes my code: 

            wave Cnumwave,DBE,Standardmass3D,Standardpeakarea3D,calibrationmasspoints //calibrationmasspoints=5, Cnumwave=23 and DBE[0]=19
            killwaves /Z responsefactoravg,responsefactorstdv
  make/n=(numpnts(Cnumwave),DBE[0]+1)responsefactoravg
  make/n=(numpnts(Cnumwave),DBE[0]+1)responsefactorstdv
  variable i,j
 
  for (i=0;i<numpnts(Cnumwave);i+=1)
  for (j=0;j<=DBE[0];j+=1)
  make/n=(calibrationmasspoints[0])tempx
  make/n=(calibrationmasspoints[0])tempy
 
  tempx = nan
  tempy = nan
  tempy=Standardmass3D[i][j][r]
  tempx=Standardpeakarea3D[i][j][r]
 
  wave W_coef,W_sigma
  string graphName="calcurve"
  DoWindow/C $graphName
  killwaves/Z fit_tempy
  make/n=5 fit_tempy
  Display/N=graphName tempy vs tempx
  K0=0
  CurveFit/H="10"/TBOX=768 line tempy /X=tempx /D
  fit_tempy=W_coef[0]+W_coef[1]*x
 

  responsefactoravg[i][j]=W_coef[1]
  responsefactorstdv[i][j]=W_sigma[1]
   Killgraphs()
endfor
endfor

Would really appreciate some advice! 

Thanks! 

A few thoughts in no specific order:

  • The declaration statement for wave W_coef, W_sigma should be put *immediately after* the operation for CurveFit. The curve fit makes the waves, and the declaration then explicitly states that they exist.

  • Why are you making and then immediately killing graphs inside the for-loop? A graph is not needed to generate W_coef and W_sigma from a curve fit.

  • Same question for fit_tempy. It seems to be extra effort for a wave (or graph) that goes through a make, fill, and kill process each step in a loop inside a loop.

  • Since tempx and tempy are 1-D waves, the proper form to generate them should likely be

    tempy = Standardmass3D[i][j][p]
    tempx = Standardpeakarea3D[i][j][p]

    (which may be the core problem)

  • Some of the processing could probably be better handled by using /FREE waves.

I might recommend that you start this with a function call to do a linear curve fit on one (x, y) pair. Pull off a vector set from the mass and area matrixes and fit them in the function call. Convince yourself that you have coded the function to give the correct answers. Then, code to have that function as a Function/WAVE return the two fit values (w_coef, w_sigma) in a return wave. Finally, use that function inside a double loop.

Function Main()
   ...
   make/N=2/FREE rtnvalues
   for (i...)
      for (j...)
         ...
         masswave = Standardmass...[i][j][p]
         areawave = Standardarea...[i][j][p]
         rtnvalues = fit_line(masswave,areawave)
         svalue[i][j] = rtnvalues[0]
         serror[i][j] = rtnvalues[1]
      endfor
   endfor
   ...
end

Function/WAVE fit_line(mw, aw)
    wave mw, aw

    make/N=2/FREE rtnwave
    ... // do curve fit, collect results, return in rtnwave
    return rtnwave
end

 

September 3, 2019 at 06:06 am - Permalink