Vtgyjfy

I have started quite recently to use Mathematica and I don´t have experience in coding.

What I need to do is to fit a peak which probably is a combination of two Voigt or Lorentzian. I have tried with code already present in the forum but without big success. Could you help me in writing this? Maybe with comments in the code so I can understand better what we are doing.

Data

data = {{19.4, 16672.}, {19.41, 16642.}, {19.42, 16778.}, {19.43, 

  16857.}, {19.44, 16833.}, {19.45, 17086.}, {19.46, 17129.}, {19.47, 

  17405.}, {19.48, 17483.}, {19.49, 17308.}, {19.5, 17884.}, {19.51, 

  17950.}, {19.52, 18202.}, {19.53, 18473.}, {19.54, 19021.}, {19.55, 

  19279.}, {19.56, 20040.}, {19.57, 20399.}, {19.58, 21412.}, {19.59, 

  22354.}, {19.6, 23334.}, {19.61, 24399.}, {19.62, 25724.}, {19.63, 

  27133.}, {19.64, 28825.}, {19.65, 30078.}, {19.66, 32224.}, {19.67, 

  33907.}, {19.68, 36299.}, {19.69, 39152.}, {19.7, 41980.}, {19.71, 

  45181.}, {19.72, 49547.}, {19.73, 55438.}, {19.74, 62094.}, {19.75, 

  69884.}, {19.76, 80306.}, {19.77, 92448.}, {19.78, 107115.}, {19.79,

   126574.}, {19.8, 148842.}, {19.81, 175298.}, {19.82, 

  205953.}, {19.83, 240900.}, {19.84, 278834.}, {19.85, 

  322364.}, {19.86, 365952.}, {19.87, 411105.}, {19.88, 

  457658.}, {19.89, 500221.}, {19.9, 544824.}, {19.91, 

  583862.}, {19.92, 619383.}, {19.93, 650362.}, {19.94, 

  672886.}, {19.95, 690179.}, {19.96, 695603.}, {19.97, 

  692265.}, {19.98, 677707.}, {19.99, 657226.}, {20., 

  630722.}, {20.01, 599184.}, {20.02, 558854.}, {20.03, 

  514989.}, {20.04, 469037.}, {20.05, 421656.}, {20.06, 

  370503.}, {20.07, 324609.}, {20.08, 278435.}, {20.09, 

  233750.}, {20.1, 195167.}, {20.11, 160965.}, {20.12, 

  131452.}, {20.13, 108026.}, {20.14, 88341.}, {20.15, 

  71993.}, {20.16, 59909.}, {20.17, 51054.}, {20.18, 44365.}, {20.19, 

  39526.}, {20.2, 36292.}, {20.21, 34308.}, {20.22, 32666.}, {20.23, 

  31599.}, {20.24, 30743.}, {20.25, 29621.}, {20.26, 29034.}, {20.27, 

  28213.}, {20.28, 27597.}, {20.29, 27485.}, {20.3, 26921.}, {20.31, 

  26588.}, {20.32, 26337.}, {20.33, 25705.}, {20.34, 26199.}, {20.35, 

  25321.}, {20.36, 25017.}, {20.37, 25011.}, {20.38, 24566.}, {20.39, 

  24232.}, {20.4, 24005.}}

My starting point is

data1 = Rest@Transpose[Rescale /@ (Transpose@data)];

peakfunc[A_, [Mu]_, [Sigma]_, x_] = A^2 E^(-((x - [Mu])^2/(2 [Sigma]^2)));



Clear[model, modelvalue]

model[data_, n_] := 

 Module[{dataconfig, modelfunc, objfunc, fitvar, fitres}, 

  dataconfig = {A[#], [Mu][#], [Sigma][#]} & /@ Range[n];

  modelfunc = (peakfunc[##, fitvar] & @@@ dataconfig // Total);

  objfunc = 

   Total[((Sqrt[data[[All, 2]]])/

       data[[All, 

         1]]) (data[[All, 2]] - (modelfunc /. fitvar -> # &) /@ 

         data[[All, 1]])^2];

  FindMinimum[objfunc, Join[{}, Flatten@dataconfig]]]

modelvalue[data_, n_] /; NumericQ[n] := 

 If[n >= 1, model[data, n][[1]], 0]

fitres = ReleaseHold[

   Hold[{Round[n], model[data1, Round[n]]}] /. 

    FindMinimum[modelvalue[data1, Round[n]], {n, 5}, 

      Method -> "PrincipalAxis"][[2]]] // Quiet



With[{n = 2}, 

 resfunc = 

  peakfunc[A[#], [Mu][#], [Sigma][#], x] & /@ Range[n] /. 

   model[data1, n][[2]]]



Show@{Plot[Evaluate[resfunc], {x, 0, 1}, 

   PlotStyle -> ({Directive[Dashed, Thick, 

         ColorData["Rainbow"][#]]} & /@ 

      Rescale[Range[Length[resfunc]]]), PlotRange -> All, 

   Frame -> True, Axes -> False], 

  Plot[Evaluate[Total@resfunc], {x, 0, 1}, 

   PlotStyle -> Directive[Thick, Red], PlotRange -> All, 

   Frame -> True, Axes -> False], 

  Graphics[{PointSize[.003], Black, Point@data1}]}

This is already not very clear to me, as you can see the fit is not very good on the right side (I expect two profiles very close as the yellow and green curve in the second picture). I would like also to know how to evaluate if the fit is good enough or not.

Many thanks!

Here is a fit using Fit and a basis of functions. I am not sure is this a solution OP wants -- may be only one Exp term is desired.

(The data variable data is taken from the question.)

(* Taken from the question. *)

Clear[peakfunc]

peakfunc[A_, [Mu]_, [Sigma]_, x_] = A^2 E^(-((x - [Mu])^2/(2 [Sigma]^2)));



(* Generate a basis of functions. *)

bFuncs = Flatten@

   Table[peakfunc[1, m, s, x], {m, 19.4, 20.4, 0.2}, {s, 0.1, 0.4, 0.1}];

Length[bFuncs]    

(* 24 *)



(* Load the QRMon package. *)

Import["https://raw.githubusercontent.com/antononcube/MathematicaForPrediction/master/MonadicProgramming/MonadicQuantileRegression.m"]



(* Do a fit with the basis, plot it, and plot the relative errors. *)

qmon =

  QRMonUnit[data]⟹

   QRMonEchoDataSummary⟹

   QRMonFit[bFuncs]⟹

   QRMonPlot⟹

   QRMonErrorPlots;

(* Take the regression function from the monad object and show it. *)

qFunc = (qmon⟹QRMonTakeRegressionFunctions)["mean"];

qFunc[x]



(* -2.445*10^7 E^(-50. (-20.4 + x)^2) + 

 3.41622*10^10 E^(-12.5 (-20.4 + x)^2) - 

 1.63757*10^11 E^(-5.55556 (-20.4 + x)^2) - 

 2.91489*10^11 E^(-3.125 (-20.4 + x)^2) + 

 1.83408*10^7 E^(-50. (-20.2 + x)^2) + 

 6.05021*10^10 E^(-12.5 (-20.2 + x)^2) - 

 4.04765*10^11 E^(-5.55556 (-20.2 + x)^2) + 

 8.90334*10^11 E^(-3.125 (-20.2 + x)^2) - 

 1.47152*10^7 E^(-50. (-20. + x)^2) + 

 5.14204*10^10 E^(-12.5 (-20. + x)^2) - 

 1.34254*10^11 E^(-5.55556 (-20. + x)^2) + 

 1.54544*10^11 E^(-3.125 (-20. + x)^2) + 

 1.35117*10^7 E^(-50. (-19.8 + x)^2) + 

 1.77025*10^10 E^(-12.5 (-19.8 + x)^2) - 

 2.38068*10^10 E^(-5.55556 (-19.8 + x)^2) - 

 1.63691*10^11 E^(-3.125 (-19.8 + x)^2) - 

 1.75233*10^7 E^(-50. (-19.6 + x)^2) - 

 1.74451*10^10 E^(-12.5 (-19.6 + x)^2) + 

 2.91813*10^11 E^(-5.55556 (-19.6 + x)^2) - 

 6.2293*10^11 E^(-3.125 (-19.6 + x)^2) + 

 3.084*10^7 E^(-50. (-19.4 + x)^2) - 

 2.02534*10^10 E^(-12.5 (-19.4 + x)^2) + 

 1.21184*10^11 E^(-5.55556 (-19.4 + x)^2) + 

 2.04278*10^11 E^(-3.125 (-19.4 + x)^2) *)

Given the large number of data points and the smoothness of the resulting curve, the most accurate fit will be using interpolation:

Show[ListPlot[data], Plot[Interpolation[data, x], {x, 19.4, 20.4}], ImageSize -> Large]

The next best fit is probably the answer by @AntonAntonov.

And notice that I am using the term "fit" meaning "a description of" rather than "an explanation of". Without a theoretical expectation for a particular curve form all one can do is cobble together some curve forms so that one can make reasonable predictions without needing all of the original data. (But computing is cheap these days: sometimes it's best just to use the original data as the summary. Why have an approximation when it's unnecessary?)

A slightly more parsimonious fit (at the expense of a less accurate fit) is to fit a mixture of a Gaussian and two Cauchy's (Cauchy <=> Lorentzian):

(* Scale the data so that the sum of the product of the bin width and the heights are 1 *)

d = data;

c = 0.01 Total[d[[All, 2]]];

d[[All, 2]] = d[[All, 2]]/c;



(* Fit a model that is a weighted mixture of Gaussian and two Cauchy-shaped curves *)

nlm = NonlinearModelFit[d, {a0 + wN PDF[NormalDistribution[μ1, σ1], x] +

    wC1 PDF[CauchyDistribution[μ2, σ2], x] +

    (1 - wN - wC1) PDF[CauchyDistribution[μ3, σ3], x],

   a0 > 0 && 0 < wC1 < 1 && 0 < wN < 1 && σ1 > 0 && 0 < wN + wC1 < 1 && σ2 > 0 && σ3 > 0},

  {{a0, d[[1, 2]]}, {wN, 0.8}, {wC1, 0.15}, {μ1, 19.9}, {σ1, 0.08},

   {μ2, 19.8}, {σ2, 0.08}, {μ3, 20}, {σ3, 0.08}}, x]



nlm["BestFitParameters"]

(* {a0 -> 0.051424552604261736`,w -> 0.7802442962949592`,μ1 -> 19.9626895200524`,

   σ1 -> 0.08312603130697435`,μ2 -> 19.905467082972322`,σ2 -> 0.19658464223570066`} *)



Show[ListPlot[data], Plot[c nlm[x], {x, 19.4, 20.4}, PlotStyle -> Red], ImageSize->Large]

How good is the fit? The two summaries of fit that I would consider are the root mean square error and a plot of the residuals vs. the predicted:

(* Root mean square error in units of counts *)

c nlm["EstimatedVariance"]^0.5

(* 1995.1131663245405` *)



ListPlot[c Transpose[{nlm["PredictedResponse"], nlm["FitResiduals"]}],

 Frame -> True, FrameLabel -> (Style[#, Bold, 18] &) /@ {"Predicted", "Residual"}]

One can see from the root mean square and the residuals vs. predicted plot that most of the predictions are within plus or minus 4,000 (in terms of the original counts). There certainly is still something potentially left to describe as one sees a pattern in the residuals.

Is this fit good enough? While I've suggested two fit summaries, it's up to you to decide on whether those metrics are small enough for your particular use of the predictions. (Alternatively one can simply state how good the fit is with those metrics and then let someone else decide if the fit is adequate for the desired objective.)

I also want to point out that despite using two curve forms that also just happen to be the shape of popular probability distributions, there is no probabilistic interpretation as if the fit were that of a probability distribution. In other words, scaling the resulting curve to integrate to 1 and then taking random samples or calculating means and variances is completely unwarranted.