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

n = 733;
nplotpoints = 15;

cumprob = {0,1,23,75,142,256,513,690,731,733}/733;
sizes = Table[(6k+1)/2,{k,0,9}];
datapoints = Map[Thread,Transpose[{sizes,cumprob}]];

length1 = Table[2,{1}];
length2 = Table[5,{22}];
length3 = Table[8,{52}];
length4 = Table[11,{67}];
length5 = Table[14,{114}];
length6 = Table[17,{257}];
length7 = Table[20,{177}];
length8 = Table[23,{41}];
length9 = Table[26,{2}];
length = Join[length1,length2,length3,
			  length4,length5,length6,
			  length7,length8,length9];


loglength = Log[length];
WW[c_] = n Log[c] + Apply[Plus,Log[length^(c-1)]] -
				   n Log[Apply[Plus,length^c]/n] - n;
uW[c_]:= (Apply[Plus,length^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,length]/n)/c;
WG[c_] = (c-1)Apply[Plus,loglength] - n c Log[uG[c]] - n Log[Gamma[c]] - n c;

cmleG = c/.FindMinimum[-WG[c],{c,{24,25}}][[2]];
smleG = uG[cmleG];
WG[cmleG];

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


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

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

deltaW = Apply[Plus,Table[(F[sizes[[k]]]-cumprob[[k]])^2,{k,1,10}]]
deltaG = Apply[Plus,Table[(G[sizes[[k]]]-cumprob[[k]])^2,{k,1,10}]]
deltaL = Apply[Plus,Table[(H[sizes[[k]]]-cumprob[[k]])^2,{k,1,10}]]

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

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


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

Show[GraphicsArray[{{upper},{lower}}],PlotLabel->
   "Figure 31.8  Comparison of Gamma, lognormal and Weibull models for minnow size data"]