# Fitting and plotting Weibull distribution

1. How can I find the best fitted 2-parameter Weibull's distribution (PDF) to describe my data? How can i plot my data and estimated distribution all in one graph?

2. How can I show/compare results of different methods of finding both weibull's parameters (maximum likelihood, linear regression...)?

3. completely optional; are there other interesting statistical informations concerning my problem that i should ask about here but I didn't due to my lack of experience with data analysis? any recommendations are very welcome and helpfull...

data = {{18.12, 0.}, {18.36, 0.}, {18.6, 0.00749}, {18.84, 0.00749}, {19.08,
0.}, {19.32, 0.00749}, {19.56, 0.03747}, {19.8, 0.02248}, {20.04,
0.02248}, {20.28, 0.08243}, {20.52, 0.05246}, {20.76,
0.02998}, {21., 0.08993}, {21.24, 0.09742}, {21.48,
0.06745}, {21.72, 0.16487}, {21.96, 0.13489}, {22.2,
0.27728}, {22.44, 0.13489}, {22.68, 0.16487}, {22.92,
0.31475}, {23.16, 0.14988}, {23.4, 0.32224}, {23.64,
0.2473}, {23.88, 0.22482}, {24.12, 0.26978}, {24.36,
0.20983}, {24.6, 0.23981}, {24.84, 0.14239}, {25.08,
0.13489}, {25.32, 0.07494}, {25.56, 0.06745}, {25.8,
0.04496}, {26.04, 0.07494}, {26.28, 0.08243}, {26.52,
0.02998}, {26.76, 0.02248}, {27., 0.01499}, {27.24,
0.02248}, {27.48, 0.02998}, {27.72, 0.02248}, {27.96,
0.00749}, {28.2, 0.}, {28.44, 0.}, {28.68, 0.}, {28.92,
0.00749}, {29.16, 0.00749}, {29.4, 0.}, {29.64, 0.}, {29.88, 0.}}

• Look up FindDistributionParameters[]. – J. M.'s technical difficulties Jul 14 '15 at 10:11
• tried that. i get that one or more data points are not in support of weibull distr... :-( – gd1980 Jul 14 '15 at 10:15
• Related: 63074, 61424 – dionys Jul 14 '15 at 10:16
• Should've mentioned that to begin with… it didn't work even with a good guess? – J. M.'s technical difficulties Jul 14 '15 at 10:21
• yeah, guess i should, sorry... i'm a newbie when it comes to mathematica... what do you mean by 'good guess'? or could it be sth wrong with my data? – gd1980 Jul 14 '15 at 10:36

Your problem is that you are not fitting raw data to a distribution, you are fitting the emperical PDF of that distribution (probably in terms of values, percentages pairs). That won't work as the functions you are using (I guess EstimatedDistributionor FindDistributionParameters) expect the raw measurement data, not frequencies.

To deal with your specific data you need either FindFit or NonlinearModelFit and fit that to the PDF of the Weibull distribution (and I don't think your data is OK for a two-parameter Weibull):

f[x_] =
NonlinearModelFit[
data,
{PDF[WeibullDistribution[α, β, μ], x], α > 0 && β > 0},
{α, β, μ}, x
]["BestFit"];

Show[
ListPlot[data],
Plot[f[x], {x, Min[data[[All, 1]]], Max[data[[All, 1]]]}]
]


This is a perfect place to use WeightedData.

dist = EstimatedDistribution[WeightedData @@ Transpose[data], WeibullDistribution[a, b]]

Show[ListPlot[data], Plot[PDF[dist, x], {x, 0, 30}]]


It is worth noting that observations with zero weights are ignored because they do not contribute to the likelihood.

• @Andy I believe what Karsten means here is that in my FindFit all numbers are treated equally, although there is more data in the higher frequencies. I guess that's true, but I suppose that the difference in modelling assumptions/methods between NonLinearModelFit and EstimatedDistribution plays a role too. – Sjoerd C. de Vries Jul 14 '15 at 23:16
• @SjoerdC.deVries I've unintentionally made things too complicated by trying to point out that the zeros get ignored. I was hoping to not get into any more detail than that :) Removing all reference to NonlinearModelFit – Andy Ross Jul 14 '15 at 23:21

As @Sjoerd C. de Vries states: one really wants to use the raw data. However, your data is binned with (I assume) the bin midpoints along with the associated relative frequency. And ideally you should account for the binning (although in this case it doesn't make much difference in the estimates and because as stated by others the fit is not hot anyway). If you had many fewer bins, then you would definitely need to account for the binning to get appropriate estimates.

To get the original integer frequencies divide by the smallest relative frequency - in this case it works because I'm sure you'll verify that the smallest positive frequency is 1 and corresponds to the relative frequency 0.07494. Further "proof" is that all of the other relative frequencies turn into numbers awfully close to integers.)

Here is the code to obtain the maximum likelihood estimates accounting for the binning of the data:

(* Bin mid-point and frequency *)
data = {{18.6, 1}, {18.84, 1}, {19.32, 1}, {19.56, 5}, {19.8,
3}, {20.04, 3}, {20.28, 11}, {20.52, 7}, {20.76, 4}, {21,
12}, {21.24, 13}, {21.48, 9}, {21.72, 22}, {21.96, 18}, {22.2,
37}, {22.44, 18}, {22.68, 22}, {22.92, 42}, {23.16, 20}, {23.4,
43}, {23.64, 33}, {23.88, 30}, {24.12, 36}, {24.36, 28}, {24.6,
32}, {24.84, 19}, {25.08, 18}, {25.32, 10}, {25.56, 9}, {25.8,
6}, {26.04, 10}, {26.28, 11}, {26.52, 4}, {26.76, 3}, {27,
2}, {27.24, 3}, {27.48, 4}, {27.72, 3}, {27.96, 1}, {28.92,
1}, {29.16, 1}};

(* Width of each bin *)
binWidth = 24/100;

(* Define log likelihood of binned data *)
logLikelihood[k_, lambda_] := Sum[data[[i, 2]] Log[
Exp[-((data[[i, 1]] - binWidth/2)/lambda)^k] -
Exp[-((data[[i, 1]] + binWidth/2)/lambda)^k]],
{i, Length[data[[All, 1]]]}]

(* Find maximum likelihood estimates *)
sol = FindMaximum[{logLikelihood[k, lambda],
k > 0 && lambda > 0}, {{k, 10}, {lambda, 24}}]


which results in

{-1910., {k -> 14.0869, lambda -> 24.205}}


A plot of the fit is given by

(* Plot results *)
p1 = Plot[(k/lambda) (x/lambda)^(k - 1) Exp[-(x/lambda)^k] /.
sol[[2]], {x, Min[data[[All, 1]]] - binWidth,
Max[data[[All, 1]]] + binWidth}];
n = Total[data[[All, 2]]];
p2 = ListPlot[
Table[{data[[i, 1]], data[[i, 2]]/(n binWidth) }, {i,
Length[data[[All, 1]]]}], AxesOrigin -> {18, 0}];
Show[{p2, p1}]


ADDITION: Estimates of variance and covariance can be obtained also (although good estimates of variance for a poor fit doesn't do one much good):

(* First partial derivatives *)
pLogL = D[logLikelihood[k, lambda], {{k, lambda}}];

(* Estimates of variances and the covariance of k and lambda *)
var = -Inverse[D[pLogL, {{k, lambda}}] /. sol[[2]]]
(* {{0.179681, 0.0108254}, {0.0108254, 0.00598091}} *)

(* Check to see that the first partial derivatives are close to zero *)
pLogL /. sol[[2]]
(* {2.13724*10^-7, -9.42956*10^-9} *)