5
$\begingroup$

Good day to all!

I have some data that I would like to fit with a normal distribution:

Ftab = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,2.15693*^18,5.87542*^18,9.75208*^18,1.35966*^19,1.74086*^19,2.14985*^19,2.57867*^19,3.06021*^19,3.56339*^19,4.08717*^19,4.63014*^19,5.14949*^19,6.30563*^19,9.59017*^19,1.29848*^20,1.64222*^20,1.98359*^20,2.31568*^20,2.6649*^20,2.93886*^20,3.12365*^20,3.30426*^20,3.47828*^20,3.64334*^20,3.79713*^20,3.93726*^20,4.08463*^20,4.28421*^20,4.46295*^20,4.61719*^20,4.74343*^20,4.83849*^20,4.89959*^20,4.89155*^20,4.8162*^20,4.70771*^20,4.56704*^20,4.3959*^20,4.19674*^20,3.97271*^20,3.72758*^20,3.46562*^20,3.16687*^20,2.8565*^20,2.54801*^20,2.24653*^20,1.95684*^20,1.68313*^20,1.42884*^20,1.19657*^20,9.88005*^19,7.96794*^19,6.3389*^19,4.97365*^19,3.84731*^19,2.93283*^19,2.20236*^19,1.62852*^19,1.1853*^19,8.48836*^18,5.91329*^18,4.08807*^18,2.77304*^18,1.84785*^18,1.20917*^18,7.76709*^17,4.89578*^17,3.02707*^17,1.83532*^17,1.09078*^17,6.36243*^16,3.63769*^16,2.037*^16,1.11682*^16,5.99341*^15,3.14729*^15,1.61677*^15,8.12254*^14,3.98982*^14,1.91568*^14,8.88572*^13,4.09217*^13,1.8515*^13,8.2291*^12,3.59248*^12,1.5403*^12,6.48546*^11,2.68139*^11,1.08849*^11,4.33808*^10,1.69724*^10,6.5182*^9,2.45709*^9,9.09066*^8,3.30082*^8,1.1809*^8,4.14603*^7,1.42842*^7,4.82907*^6,1.60192*^6,521400.,166763.,52889.9,16659.9,5212.2,1619.74,499.999,153.326,46.7099,14.1375,4.25138,1.27029,0.37715,0.111272,0.0326245,0.00950618,0.00275293,0.000792364,0.000227188,0.0000647492,0.000018344,5.16638*^-6,1.44746*^-6,4.04993*^-7,1.13245*^-7,3.16491*^-8,8.84104*^-9,2.48832*^-9,6.89887*^-10,1.92817*^-10,5.3877*^-11,1.50516*^-11,4.19478*^-12,1.1663*^-12,3.23529*^-13,8.95459*^-14,2.47316*^-14,6.81648*^-15,1.87536*^-15,5.15688*^-16,1.4177*^-16,3.89675*^-17,1.07096*^-17,2.94315*^-18,8.08787*^-19,2.22263*^-19,6.10835*^-20,1.67907*^-20,4.61637*^-21,1.26962*^-21,3.49314*^-22,9.61425*^-23,2.64971*^-23,7.32514*^-24,2.03164*^-24,5.65294*^-25,1.57813*^-25,4.42091*^-26,1.24274*^-26,3.50629*^-27,9.93014*^-28,2.8206*^-28,8.03604*^-29,2.29661*^-29,6.58427*^-30,1.89381*^-30,5.46852*^-31,1.58502*^-31,4.61092*^-32,1.34633*^-32,3.94568*^-33,1.16067*^-33,3.42722*^-34,1.01591*^-34,3.02326*^-35,9.03278*^-36,2.7097*^-36,8.16238*^-37,2.46913*^-37,7.50123*^-38,2.28876*^-38,7.01369*^-39,2.15558*^-39,6.63493*^-40,2.04544*^-40,6.31601*^-41,1.95371*^-41,6.05436*^-42,1.87969*^-42,5.84708*^-43,1.82244*^-43,5.69175*^-44,1.78132*^-44,5.58675*^-45,1.75599*^-45,5.53156*^-46,1.74646*^-46,5.52685*^-47,1.75316*^-47,5.57461*^-48,1.777*^-48,5.67974*^-49,1.82043*^-49,5.85113*^-50,1.8926*^-50,6.20322*^-51,2.03881*^-51,6.71987*^-52,2.22122*^-52,7.36367*^-53,2.44869*^-53,8.1692*^-54,2.73435*^-54,9.18288*^-55,3.0944*^-55,1.04633*^-55,3.55043*^-56,1.20901*^-56,4.13181*^-57,1.41721*^-57,4.87905*^-58,1.68614*^-58,5.79854*^-59,1.99629*^-59,6.90152*^-60,2.3963*^-60,8.35665*^-61,2.92717*^-61,1.03004*^-61,3.64151*^-62,1.29343*^-62,4.60958*^-63,1.64609*^-63,5.89022*^-64,2.11209*^-64,7.58939*^-65,2.73295*^-65,9.86276*^-66,3.56716*^-66,1.29306*^-66,4.69784*^-67,1.71071*^-67,6.24407*^-68,2.28447*^-68,8.37801*^-69,3.08*^-69,1.13508*^-69,4.19358*^-70,1.55323*^-70,5.7676*^-71,2.14725*^-71,8.0159*^-72,3.00073*^-72,1.12647*^-72,4.24074*^-73,1.60107*^-73,6.06226*^-74,2.30213*^-74,8.76822*^-75,3.34958*^-75,1.28346*^-75,4.93351*^-76,1.90206*^-76,7.34281*^-77,2.83763*^-77,1.09777*^-77,4.25153*^-78,1.64839*^-78,6.39831*^-79,2.48641*^-79,9.67365*^-80,3.76821*^-80,1.46972*^-80,5.73981*^-81,2.24458*^-81,8.78933*^-82,3.44642*^-82,1.35326*^-82,5.32112*^-83,2.09529*^-83,8.26251*^-84,3.26937*^-84,1.30314*^-84,5.17437*^-85,2.04684*^-85,8.1092*^-86,3.21771*^-86,1.2788*^-86,5.09376*^-87,2.03636*^-87,8.15419*^-88,3.27042*^-88,1.31379*^-88,5.28637*^-89,2.13063*^-89,8.60167*^-90,3.47851*^-90,1.40912*^-90,5.7181*^-91,2.32444*^-91,9.46562*^-92,3.8613*^-92,1.57788*^-92,6.45927*^-93,2.64891*^-93,1.08826*^-93,4.4791*^-94,1.84692*^-94,7.62982*^-95,3.15773*^-95,1.30778*^-95,5.41761*^-96,2.24526*^-96,9.30955*^-97,3.86189*^-97,1.60282*^-97,6.65567*^-98,2.76519*^-98,1.14936*^-98,4.77753*^-99,1.98643*^-99,8.26825*^-100,3.4443*^-100,1.43547*^-100,5.98552*^-101,2.49705*^-101,1.04227*^-101,4.3527*^-102,1.81875*^-102,7.6038*^-103,3.18078*^-103,1.33134*^-103,5.57575*^-104,2.33657*^-104,9.79784*^-105,4.11124*^-105,1.72629*^-105,7.25367*^-106,3.05037*^-106,1.28386*^-106,5.40825*^-107,2.28023*^-107,9.62247*^-108,4.06431*^-108,1.71825*^-108,7.2709*^-109,3.07963*^-109,1.30564*^-109,5.54077*^-110,2.35364*^-110,1.00079*^-110,4.25977*^-111,1.81498*^-111,7.74111*^-112,3.30513*^-112,1.42425*^-112,6.16014*^-113,2.62352*^-113,1.1185*^-113,4.77361*^-114,2.04521*^-114,8.7969*^-115,3.78806*^-115,1.63307*^-115,7.04855*^-116,3.04587*^-116,1.31808*^-116,5.71268*^-117,2.47975*^-117,1.07809*^-117,4.6944*^-118,2.04736*^-118,8.94306*^-119,3.90966*^-119,1.71011*^-119,7.48426*^-120,3.27728*^-120,1.4359*^-120,6.2948*^-121,2.76118*^-121,1.21189*^-121,5.32222*^-122,2.33877*^-122,1.02837*^-122}

I have tried three different approaches:

pars = FindDistributionParameters[Ftab,NormalDistribution[a1, a2],ParameterEstimator->"MethodOfMoments"];

(*{a1->2.61159*^19,a2->9.27508*^19}*)


epstab = Table[eps, {eps, 0, 50,.1}];

Ffit = FindFit[ Transpose@Join[{epstab}, {Ftab}], a1 PDF[NormalDistribution[a2, a3], eps],{a1, a2, a3}, eps]
(*{a1->1.,a2->1.,a3->1.} *)



nlm = NonlinearModelFit[Ftab, k Exp[-(eps - m1)^2/(2 s^2)], {m1, s, k}, eps];
(*{m1->1.,s->0.5,a3->1.} *)

However, these are not even close to the actual data.

I have tried even with the truncated data

Ftabsel=Select[Ftab, # >= 10^-21 &];

Pretty much the same.

$\endgroup$
7
  • $\begingroup$ While you only give a one-dimensional list of values, your data is (apparently) data = Transpose[{Range[Length[Ftab]], Ftab}];. Therefore you have a regression that you want to perform with a normal curve shape as opposed to a random sample from a normal distribution. In that case FindDistributionParameters is not at all what you want. $\endgroup$
    – JimB
    Mar 23, 2022 at 4:21
  • $\begingroup$ More specifically, Transpose[{epstab, Ftab}] with epstab = Table[eps, {eps, 0, 50,.1}]; $\endgroup$ Mar 23, 2022 at 4:29
  • 1
    $\begingroup$ Sorry, I should read more carefully. That's exactly what you have in the question. (The mixing up of regression and fitting a probability distribution tends to set me off.) $\endgroup$
    – JimB
    Mar 23, 2022 at 4:34
  • $\begingroup$ I tried FindFit which is supposed to do nonlinear regression, right? $\endgroup$ Mar 23, 2022 at 4:43
  • $\begingroup$ Theoretically normal dist is never zero. $\endgroup$ Mar 23, 2022 at 4:54

4 Answers 4

6
$\begingroup$

Consider the following truncated data set:

epstab = Table[eps, {eps, 0, 50, .1}];
data = Transpose@Join[{epstab}, {Ftab}];
data = Select[data, 10 < #[[1]] < 20 &];

Now fit a normal/Gaussian shaped curve to that data:

nlm1 = NonlinearModelFit[data, a1 Exp[-(x - b1)^2/c1], {{a1, 5 10^20}, {b1, 15}, {c1, 2}}, x];

nlm1["EstimatedVariance"]^0.5
(* 1.29432*10^19 *)

Show[ListPlot[data, PlotRange -> All], Plot[nlm1[x], {x, 10, 20}]]

Data and fit with 1 Gaussian-shaped curve

Now plot the residuals vs the predictors to see more clearly where the lack-of-fit is located:

ListPlot[Transpose[{data[[All, 1]], nlm1["FitResiduals"]}], 
 PlotRange -> {All, {-5 10^19, 5 10^19}}, Frame -> True,
 FrameLabel -> {"x", "Residual"}]

Residuals vs predictors for 1 gaussian-shaped curve

If that fit is good enough with respect to your fit criteria, then you are done. But if not, suppose we add in another Gaussian-shaped curve:

nlm2 = NonlinearModelFit[data, a1 Exp[-(x - b1)^2/c1] + a2 Exp[-(x - b2)^2/c2],
  {{a1, 5 10^20}, {a2, 10^20}, {b1, 15}, {b2, 13.5}, {c1, 2}, {c2, 0.26}}, x];

nlm2["EstimatedVariance"]^0.5
(* 4.47305*10^18 *)

This measure of goodness-of-fit is about one-third the size of the previous model. The residuals vs the predictors looks like the following:

Show[ListPlot[data, PlotRange -> All], Plot[nlm2[x], {x, 10, 20}]]

Data and fit for 2 Gaussian-shaped curves

ListPlot[Transpose[{data[[All, 1]], nlm2["FitResiduals"]}], 
 PlotRange -> {All, {-5 10^19, 5 10^19}}, Frame -> True,
 FrameLabel -> {"x", "Residual"}]

Residuals vs predicted values for 2 Gaussian-shaped curves

Is that adequate? I don't know.

$\endgroup$
2
  • $\begingroup$ Say, I have 100 sets of data and want to repeat what you did. Is there a way to pick the endpoints 10 and 20, as well as the mean and the variance in some sense automatic? $\endgroup$ Mar 23, 2022 at 23:06
  • $\begingroup$ I used EmpiricalDistribution from @Roman answer and substituted the initial guesses for a1 by Max@Ftab, for b1 by Mean[d] and for c1 by StandardDeviation[d]. At the same time, I have updated your function as follows: a1/Sqrt[2Pi]/StandardDeviation[d] Exp[-(eps - b1)^2/c1^2]. $\endgroup$ Mar 24, 2022 at 0:41
5
$\begingroup$

Edit 2:

You can also multiply the model by max of the original data to fit the original data.

Show[ListPlot[data, PlotRange -> All], 
 Plot[max nlm[x - zeros], {x, 0, 500}, PlotStyle -> Red, 
  PlotRange -> All], Frame -> True]

enter image description here Edit:

The issue is 0.'s. We can ignore them, and shift the model by the number of 0.'s

zeros = Count[data, 0.]

117

We can truncate zeros and normalize the data.

data = data/Max@data;
dataTrancated = Select[Ftab, # >= 10^-122 &];

nlm = NonlinearModelFit[dataTrancated, 
   k Exp[-(1/2) ((x - m)/s)^2], {m, s, k}, x, 
   Method -> {"NMinimize", Method -> "DifferentialEvolution"}];
Normal[nlm]

Show[ListPlot[data, PlotRange -> All], 
 Plot[nlm[x - zeros], {x, 0, 500}, PlotStyle -> Red, 
  PlotRange -> All], Frame -> True]

enter image description here Original Answer: You can normalize the data and truncate it.

Ftab = Ftab/Max@Ftab;
Ftab = Select[Ftab, # >= 10^-20 &];

nlm = NonlinearModelFit[Ftab, k Exp[-(1/2) ((x - m)/s)^2], {m, s, k}, 
  x, Method -> {"NMinimize", Method -> "DifferentialEvolution"}];
Normal[nlm]

$1.00685 e^{-0.00432893 (x-32.3828)^2}$

Show[ListPlot[Ftab], Plot[nlm[x], {x, 0, 120}, PlotStyle -> Red], 
 Frame -> True]

enter image description here

$\endgroup$
5
  • $\begingroup$ Your solution indeed works. Thanks. Due to the truncation, the peak of the fit has been shifted from the actual data to the right. Is there a way to transform it and obtain a fit for the initial data? $\endgroup$ Mar 23, 2022 at 4:16
  • $\begingroup$ @AsaturKhurshudyan What are your criteria for an adequate fit? Is there something concrete or is it "I'll know it when I see it." $\endgroup$
    – JimB
    Mar 23, 2022 at 4:23
  • $\begingroup$ @JimB That's a good question I hadn't asked to myself. I think, minimum l2 norm should work. $\endgroup$ Mar 23, 2022 at 4:27
  • $\begingroup$ That's a reasonable criterion. Now: what value of that would result in an adequate fit? $\endgroup$
    – JimB
    Mar 23, 2022 at 4:31
  • $\begingroup$ @JimB I would expect it to be small in the vicinity of the peak, e.g., 0.01. For, the total l2 norm should be quite big. $\endgroup$ Mar 23, 2022 at 4:54
2
$\begingroup$

Build a distribution that represents your data:

d = EmpiricalDistribution[Ftab -> Range[0, 50, 1/10]];

Construct a normal distribution with the same mean and standard deviation:

n = NormalDistribution[Mean[d], StandardDeviation[d]]
(*    NormalDistribution[14.8267, 1.01774]    *)

Show that the cumulative distribution functions differ by less than 3%, which is great (see Kolmogorov–Smirnov test):

Plot[CDF[d, x] - CDF[n, x], {x, 0, 50}, PlotRange -> All]

enter image description here

$\endgroup$
2
  • $\begingroup$ If Ftab is considered proportional to the counts associated with the 501 categories (assuming this set of data is a random sample from some distribution), then the Kolmogorov-Smirnov statistic is certainly a measure of the goodness-of-fit. However, because the sample size is unknown, the associated P-value cannot be determined. But if the smallest non-zero value (1.02837 * 10^-122) is associated even with a count of 1, then the P-value will be essentially zero as the sample size would have to be 1.27231 * 10^144. $\endgroup$
    – JimB
    Mar 24, 2022 at 0:28
  • $\begingroup$ @JimB That's exactly why I did not compute a P-value. What I show is, however, inspired by Kolmogorov–Smirnov as we're comparing two CDF. $\endgroup$
    – Roman
    Mar 24, 2022 at 6:10
2
$\begingroup$
data = Ftab  // Chop // N ;

(* normalize data *)
peak = Max[data]  ;
norm = data/peak ;

(* mean initial guess *)
mean =First[Flatten[ Position[norm, Max[norm]]]]  ;

(* sigma initial guess *)
{a, b} =  Select[norm, # > 0.5*norm[[mean]] &][[{1, -1}]] ;
{min, max} =Sort[ Map[First[Flatten[Position[norm, #]]] &, {a, b}]] ;
fwhm = max - min ;
sigma = fwhm/2.36 ;

(* select data *)
factor = 3 ;
{min, max} = Floor[{ mean - factor*sigma ,  mean + factor*sigma }] ;
x = Range[min, max] ;
y = Take[norm, {min, max}] ;

(* fit *)
fit = NonlinearModelFit[Transpose[{x, y}], const*Exp[-(t - m)^2/(2*s^2)], {{m, mean}, {s, sigma}, {const, 1.0}}, t] ;
fit["BestFitParameters"]

Show[ListPlot[data, PlotRange->{{min, max}, All}], Plot[peak*fit[t], {t, min, max}]]
$\endgroup$
5
  • $\begingroup$ What is 1/2.36 in the definition of sigma? $\endgroup$ Mar 28, 2022 at 1:29
  • $\begingroup$ Your solution works, thanks. However, it is shifted to the right and the horizontal axis is not with my initial epstab. $\endgroup$ Mar 28, 2022 at 1:33
  • $\begingroup$ @AsaturKhurshudyan, 2.36 is a relation between fwhm and sigma $\endgroup$
    – I.M.
    Mar 28, 2022 at 1:57
  • $\begingroup$ I don't get it, sorry $\endgroup$ Mar 28, 2022 at 2:01
  • $\begingroup$ @AsaturKhurshudyan, not sure what you mean by 'sighted to the right'. My point here is you need a good initial guess for parameters for nonlinear fit to work well. Mean, sigma and data range initial values estimation can be automated. Initial guess for sigma can be inferred from fwhm. $\endgroup$
    – I.M.
    Mar 28, 2022 at 2:05

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.