1
$\begingroup$

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

error

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]
$\endgroup$
4
  • 2
    $\begingroup$ "It doesn't work" is not very informative. What goes wrong? Do you get errors? Are you unhappy with the quality of the fit? Have you tried providing good initial guesses for the values of your fit parameters, rather than relying on the defaults? $\endgroup$
    – MarcoB
    Commented Nov 20, 2021 at 4:09
  • $\begingroup$ @MarcoB I put guesses: 1.173, 1.332 but still I don't get doble gaussian fit. $\endgroup$ Commented Nov 20, 2021 at 18:37
  • $\begingroup$ We need to know what is piki $\endgroup$ Commented Nov 20, 2021 at 22:32
  • $\begingroup$ @JackLaVigne I've corrected the name. $\endgroup$ Commented Nov 20, 2021 at 22:35

1 Answer 1

3
$\begingroup$

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]

Mathematica graphics

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]
 ]

Mathematica graphics

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}
 ]

Mathematica graphics

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}
 ]

Mathematica graphics

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]
  ]
 ]

Mathematica graphics

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.