ClearAll;
$DefaultFont = {"Times-Roman",10};
Needs["Statistics`ContinuousDistributions`"];
Needs["Graphics`Graphics`"];

n = 489;
nplotpoints = 100;

cumulative = {0,0,0,
			  9/489,14/489,42/489,87/489,
			  252/489,342/489,460/489,477/489,
			  486/489,486/489,487/489,1};
datapoints = {{0,0},{1/60,0},{2/60,0},{3/60,0},
			  {4/60,9/489},{5/60,14/489},{6/60,42/489},{7/60,87/489},
			  {8/60,252/489},{9/60,342/489},{10/60,460/489},{11/60,477/489},
			  {12/60,486/489},{13/60,486/489},{14/60,487/489},{15/60,1}};
correcteddatapoints = Drop[datapoints+Join[{{0,0}},Table[{1/120,0},{k,1,15}]],3];

thick1 = Table[4,{9}];
thick2 = Table[5,{5}];
thick3 = Table[6,{28}];
thick4 = Table[7,{45}];
thick5 = Table[8,{165}];
thick6 = Table[9,{90}];
thick7 = Table[10,{118}];
thick8 = Table[11,{17}];
thick9 = Table[12,{9}];
thick10 = Table[14,{1}];
thick11 = Table[15,{2}];
thick = N[Join[thick1,thick2,thick3,
			 thick4,thick5,thick6,
			 thick7,thick8,thick9,
			 thick10,thick11]/60];
			 logthick = Log[thick];
WW[c_] = n Log[c] + Apply[Plus,Log[thick^(c-1)]] -
				   n Log[Apply[Plus,thick^c]/n] - n;
uW[c_]:= (Apply[Plus,thick^c]/n)^(1/c);
cmleW = c/.FindMinimum[-WW[c],{c,5,6}][[2]]
smleW = uW[cmleW]
WW[cmleW]

F[x_]:=  1 - E^(-(x/smleW)^cmleW);
f[x_]:= cmleW/smleW (x/smleW)^(cmleW-1) E^(-(x/smleW)^cmleW);


uG[c_]:= (Apply[Plus,thick]/n)/c;
WG[c_] = (c-1)Apply[Plus,logthick] - n c Log[uG[c]] - n Log[Gamma[c]] - n c;

cmleG = c/.FindMinimum[-WG[c],{c,{27,28}}][[2]]
smleG = uG[cmleG]
WG[cmleG]

G[x_]:=  CDF[GammaDistribution[cmleG,smleG],x];
g[x_]:= PDF[GammaDistribution[cmleG,smleG],x];


ubar = Apply[Plus,logthick]/n
sdmle = Sqrt[Apply[Plus,(logthick-ubar)^2]/n]
dist = LogNormalDistribution[ubar,sdmle];

H[x_]:=  CDF[dist,x];
h[x_]:= PDF[dist,x];

deltaW = Apply[Plus,Table[(F[(2k+1)/120]-cumulative[[k]])^2,{k,3,15}]]
deltaG = Apply[Plus,Table[(G[(2k+1)/120]-cumulative[[k]])^2,{k,3,15}]]
deltaL = Apply[Plus,Table[(H[(2k+1)/120]-cumulative[[k]])^2,{k,3,15}]]

f[x]
g[x]
h[x]

pdfW = Plot[f[x],{x,0,1.82},PlotPoints -> nplotpoints,
			PlotStyle->Dashing[{0.01,0.01}]];
pdfG = Plot[g[x],{x,0,1.82},PlotPoints -> nplotpoints,
			PlotStyle->Dashing[{0.001,0.01}]];
pdfL = Plot[h[x],{x,0.0001,1.82},PlotPoints -> nplotpoints];

cdfW = Plot[F[x],{x,0,1.82},PlotPoints -> nplotpoints,PlotStyle->Dashing[{0.01,0.01}]];
cdfG = Plot[G[x],{x,0,1.82},PlotPoints -> nplotpoints,PlotStyle->Dashing[{0.001,0.01}]];
cdfL = Plot[H[x],{x,0.0001,1.82},PlotPoints -> nplotpoints];


upper =  Show[pdfW,pdfG,pdfL,
			PlotRange->{{0,0.26},{-0.01,15.5}},
			AxesOrigin->{-0.01,-0.035},
			AxesLabel->{"x","y"},
			PlotLabel->"(a)  y = f(x)"];
lower =  Show[cdfW,cdfG,cdfL,ListPlot[correcteddatapoints],
			PlotRange->{{0,0.26},{-0.01,1.01}},
			AxesOrigin->{-0.01,-0.035},
			AxesLabel->{"x","y"},
			PlotLabel->"(b)  y = Int(f,[0,x])"];

Show[GraphicsArray[{{upper},{lower}}],PlotLabel->
   "Figure 31.7  Comparison of Gamma, lognormal and Weibull models for leaf thickness data"]