Nonlinear model fit gaussian precision

Question

I wanted to fit double gaussian but it doesn't work, it doesn't fit gaussian. The error I get:

Data:data Here is code

peaks = {{398, 0.511}, {495, 0.662}, {819, 1.173}, {921, 1.332}}

lm = LinearModelFit[peaks, x, x]

co60pl1 = 
  Import[
   "/Users/co60_pl1_1200s.lst", "Table"];

nlmcopl1 = 
 NonlinearModelFit[
  Select[{lm[First[#]], Last[#]} & /@ co60pl1, 1 <= #[[1]] <= 1.5 &], 
  a3*E^(-((x - 1.173)^2/(2*b3^2))) + b3*x + c3 + 
   a4*E^(-((x - 1.332)^2/(2*b4^2))), {a3, b3, c3, a4, b4}, x, 
  MaxIterations -> 1000]

Answer 1

I see a couple of problems. The means that you are using, 1.173 and 1.332 are reasonable guesses, but I suggest letting NonlinearModelFit find those parameters. A big mistake is using the same parameter, b3, for the variance of the first Gaussian and the slope of the line. Another problem is that NonlinearModelFit is going to need some reasonable starting parameters for the model.

Here is a plot of your data.

co60pl1 = Import["co60_pl1_1200s.lst", "Table"];

ListPlot[co60pl1, PlotStyle -> Black]

Now you use peaks to shift the x axis.

peaks = {{398, 0.511}, {495, 0.662}, {819, 1.173}, {921, 1.332}}

lm = LinearModelFit[peaks, x, x]

This enables us to create a new set of data with the x-axis shifted. I'll call it shiftedData.

shiftedData = {lm[First[#]], Last[#]} & /@ co60pl1;

Next step is to limit it to where x is between 1 and 1.5. For reason's I will explain later, I also created data limited between 1.07 and 1.5.

shiftedData1 = Select[shiftedData, 1 <= #[[1]] <= 1.5 &];
shiftedData2 = Select[shiftedData, 1.07 <= #[[1]] <= 1.5 &];

Below is a plot of these three data sets. Note that shiftedData is identical to co60pl12 except with the x-axis shifted.

Show[
 ListPlot[shiftedData, PlotStyle -> Black],
 ListPlot[shiftedData1, PlotStyle -> Red],
 ListPlot[shiftedData2, PlotStyle -> Blue]
 ]

Initially I am going to use shiftedData2 and try to fit two Gaussians to it ignoring the linear portion.

I use Manipulate to get parameter values that are in the ball park.

Manipulate[
 Show[
  ListPlot[shiftedData2, PlotStyle -> Blue],
  Plot[a1*Exp[-(x - μ1)^2/(2*σ1^2)] + 
       a2*Exp[-(x - μ2)^2/(2*σ2^2)],
   {x, 1, 1.5}, PlotStyle -> Black]
  ],
 {{a1, 220}, 100, 300},
 {{μ1, 1.14}, 1.1, 1.3},
 {{σ1, 0.05}, 0.02, 0.2},
 {{a2, 160}, 100, 300},
 {{μ2, 1.3}, 1.2, 1.4},
 {{σ2, 0.04}, 0.02, 0.2}
 ]

Now refine the parameters by using NonlinearModelFit and supply the Manipulate parameters as starting values.

nlm2 = NonlinearModelFit[
  shiftedData2,
  a1*Exp[-(x - μ1)^2/(2*σ1^2)] + 
  a2*Exp[-(x - μ2)^2/(2*σ2^2)],
  {
   {a1, 220},
   {μ1, 1.14},
   {σ1, 0.05},
   {a2, 160},
   {μ2, 1.3},
   {σ2, 0.04}
   },
  x
  ]

nlm2["BestFitParameters"]
(* {a1 -> 214.277, μ1 -> 1.13994, σ1 -> 0.0532518, 
    a2 -> 155.572, μ2 -> 1.29863, σ2 -> 0.037658} *)

Next use the shiftedData1 and the parameters from nlm2 and guess the slope and intercept for the line using Manipulate.

Manipulate[
 Show[
  ListPlot[shiftedData1, PlotStyle -> Red],
  Plot[m x + b + a1*Exp[-(x - μ1)^2/(2*σ1^2)] + 
    a2*Exp[-(x - μ2)^2/(2*σ2^2)],
   {x, 1, 1.5}, PlotStyle -> Black]
  ],
 {{m, -260}, -2000, 0},
 {{b, 330}, 0, 2000},
 {{a1, 214}, 100, 300},
 {{μ1, 1.14}, 1.1, 1.3},
 {{σ1, 0.053}, 0.02, 0.2},
 {{a2, 155}, 100, 300},
 {{μ2, 1.3}, 1.2, 1.4},
 {{σ2, 0.038}, 0.02, 0.2}
 ]

Final refinement is to fit the linear and Gaussian using NonlinearModelFit and supplying it with reasonable starting values.

nlm1 = NonlinearModelFit[
  shiftedData1,
  m x + b + a1*Exp[-(x - μ1)^2/(2*σ1^2)] + 
   a2*Exp[-(x - μ2)^2/(2*σ2^2)],
  {
   {m, -260},
   {b, 330},
   {a1, 214},
   {μ1, 1.14},
   {σ1, 0.053},
   {a2, 160},
   {μ2, 1.3},
   {σ2, 0.038}
   },
  x
  ]

nlm1["BestFitParameters"]
(* {m -> -260.749, b -> 379.83, 
    a1 -> 139.05,  μ1 -> 1.14737, σ1 -> 0.0357556, 
    a2 -> 118.309, μ2 -> 1.29887, σ2 -> 0.0319747} *)

A plot of the shiftedData1 and the model shows a reasonable fit (the line from NonlinearModelFit is superimposed in Gray).

With[
 {
  m = m /. nlm1["BestFitParameters"],
  b = b /. nlm1["BestFitParameters"]
  },
 Show[
  Plot[nlm1[x], {x, 1, 1.5}, PlotStyle -> Black],
  Plot[m x + b, {x, 1, 1.5}, PlotStyle -> Gray],
  ListPlot[shiftedData1, PlotStyle -> Red]
  ]
 ]