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

n = 22;

xcum = {0.1875,0.3125,0.4375,0.5625,0.8125,1.0625,1.3125,1.5625,1.8125};
ycum = {0,2,6,13,16,19,20,21,22}/22;
datapoints = Table[{xcum[[i]],ycum[[i]]},{i,1,9}];

ratpupilradii = {0.25,0.25,0.375,0.375,0.375,0.375,0.5,0.5,0.5, 0.5,0.5,
				 0.5, 0.5, 0.75, 0.75, 0.75, 1.0,  1.0,1.0,1.25,1.5,1.75};
logratpupilradii = Log[ratpupilradii];
WW[c_] = n Log[c] + Apply[Plus,Log[ratpupilradii^(c-1)]] -
				   n Log[Apply[Plus,ratpupilradii^c]/n] - n;
uW[c_]:= (Apply[Plus,ratpupilradii^c]/n)^(1/c);
cmleW = c/.FindMinimum[-WW[c],{c,1.9}][[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,ratpupilradii]/n)/c;
WG[c_] = (c-1)Apply[Plus,logratpupilradii] - n c Log[uG[c]] - n Log[Gamma[c]] - n c;

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

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


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

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

deltaW = Apply[Plus,Table[(F[xcum[[k]]]-ycum[[k]])^2,{k,1,9}]]
deltaG = Apply[Plus,Table[(G[xcum[[k]]]-ycum[[k]])^2,{k,1,9}]]
deltaL = Apply[Plus,Table[(H[xcum[[k]]]-ycum[[k]])^2,{k,1,9}]]

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

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


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

Show[GraphicsArray[{{upper},{lower}}],PlotLabel->
   "Figure 31.6  Comparison of Gamma, lognormal and Weibull models for rat pupil data"]