1
$\begingroup$

I have the following set of data:

dataNorm =
{{5., 0.}, {10., 0.00615385}, {15., 0.0184615}, {20., 0.0246154}, {25., 0.0461538}, 
{30., 0.0707692}, {35., 0.138462}, {40., 0.135385}, {45., 0.147692}, {50., 0.107692}, 
{55., 0.12}, {60., 0.0553846}, {65., 0.0369231}, {70., 0.0307692}, {75., 0.0123077}, 
{80., 0.0153846}, {85., 0.0184615}, {90., 0.00615385}, {95., 0.00307692}, {100., 0.00615385}}

which values are coming from experiments, creating "bins" from 5 to 100 (with step of 5) and computing the relative frequency of each case. The total number of measurements is 325. The objective is to fit these values with a LogNormal distribution:

model1 = PDF[LogNormalDistribution[a, b], t];
fit1 = FindFit[dataNorm, model1, {a,b}, t]
func1 = Function[{t}, Evaluate[model1 /. fit1]];
(* {a -> 3.74074, b -> 0.146623} *)     

The resulting plot appears as:

Show[Plot[func1[t], {t, 0, 100}, PlotStyle -> ColorData[112, "ColorList"][[1]]], ListPlot[dataNorm], PlotRange -> All]

Plot of first fit

The result is very off, even if the shape and position of the resulting distribution seems to be all right. Modifying the argument in the LogNormal distribution in order to have some sort of scaling and shifting, the result is:

model2 = PDF[LogNormalDistribution[a, b], c t - d];
fit2 = FindFit[dataNorm, model2, {a, b, c, d}, t ] 
func2 = Function[{t}, Evaluate[model2 /. fit2]];
(* {a -> 2.94948, b -> 0.146413, c -> 0.208136, d -> -9.86194} *)

Modified distribution

Question 1: I admit that I do not know much about properties of the LogNormal distribution, but why is it necessary to modify the argument of the function to make it fit?

Then, in order to check if the solution is a proper probability distribution, the CDF is computed:

model2CDF = CDF[LogNormalDistribution[a, b], c t - d] func2CDF = Function[{t}, Evaluate[model5C /. fit5]];

Cumulative density function

The cumulative density function tends to 1 but func2CDF[100] = 0.999397 and not 1. In addition, the integral of the PDF is Integrate[func5[t], {t, 0, 100}] = 4.80164. Even if the function for the fit is truncated, such that

model2Truncated = PDF[TruncatedDistribution[{0, 100}, LogNormalDistribution[a, b]], c t - d]

the CDF do not reach 1 at 100 and the integration is the same.

Considering that the cumulative frequency of the data that was provided before is

{0., 0.00615385, 0.0246154, 0.0492308, 0.0953846, 0.166154, 0.304615,
0.44, 0.587692, 0.695385, 0.815385, 0.870769, 0.907692, 0.938462,
0.950769, 0.966154, 0.984615, 0.990769, 0.993846, 1.}

Question 2: am I missing something or there is an error in the procedure that I used to compute the distribution function?

$\endgroup$
  • $\begingroup$ You should avoid "regression" functions (such as NonlinearModelFit and FindFit) when you have (or think you have) a random sample from a probability distribution. Having the raw data would be better and then using FindDistributionParameters would be appropriate (assuming you've chosen an appropriate probability distribution). Alternatively when one only has the resulting histogram, then one needs to use the actual counts rather than the relative frequencies. $\endgroup$ – JimB Sep 3 '18 at 1:57
  • $\begingroup$ I agree completely with JimB. Anyway, since a pdf is always normalized in the sense Integral f[x] dx =1, the ordinates must be divided by the bin width (the "dx"). If you use dataNorm[[All, 2]] /= 5 before the fit, the result is better, but not the best, as JimB pointed out. $\endgroup$ – Vito Vanin Sep 3 '18 at 2:08
2
$\begingroup$
(* I assume that the first element is the midpoint of the bin and the 
   second is the relative frequency based on 325 observations.
   So transform back to the original counts *)
data = dataNorm;
data[[All, 2]] = Round[325 data[[All, 2]]];

(* Create the log of the likelihood for the histogram summary *)
cdf[x_, μ_, σ_] := 
 CDF[LogNormalDistribution[μ, σ], x]
halfWidth = 2.5;
logL[μ_, σ_, shift_] := 
 Sum[data[[i, 2]] Log[cdf[data[[i, 1]] + halfWidth - shift, μ, σ] - 
     cdf[data[[i, 1]] - halfWidth - shift, μ, σ]], {i, Length[data]}]

(* Make guesses as to initial values *)
shift0 = 2.5;
sampleMean = Total[(data[[All, 1]] - shift0) data[[All, 2]]]/Total[data[[All, 2]]];
sampleVariance = Total[(data[[All, 1]] - shift0)^2 data[[All, 2]]]/
   Total[data[[All, 2]]] - sampleMean^2;
μ0 = Log[sampleMean^2/Sqrt[sampleMean^2 + sampleVariance]];
σ0 = Sqrt[2] Sqrt[Log[Sqrt[sampleMean^2 + sampleVariance]/sampleMean]];
(* Get maximum likelihood estimates *)
mle = FindMaximum[{logL[μ, σ, shift], {σ > 0 && shift <= 2.5}}, 
  {{shift, shift0}, {μ, μ0}, {σ, σ0}}]
(* {-825.554, {shift -> -39.1882, μ -> 4.42897, σ -> 0.182086}} *)

(* Plot the results *)
Show[ListPlot[Transpose[{data[[All, 1]], data[[All, 2]]/(5*325)}], 
  PlotRange -> {{0, 120}, All}],
 Plot[PDF[LogNormalDistribution[μ, σ], x - shift] /. 
   mle[[2]], {x, shift /. mle[[2]], 120}]]

Lognormal distribution fit

(* Now any good analysis requires (in my world) estimates of precision.
   Here is an approach to estimate standard errors using the Delta Method. *)
covMatrix = -Inverse[(D[logL[μ, σ, shift], {{μ, σ, shift}, 2}]) /. mle[[2]]];
seμ = covMatrix[[1, 1]]^0.5
(* 0.229514 *)
seσ = covMatrix[[2, 2]]^0.5
(* 0.0422745 *)
seShift = covMatrix[[3, 3]]^0.5
(* 18.9102 *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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