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}];
datapoints2 = Table[{Pi xcum[[i]]^2,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];

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}]]


left = Plot[WW[c],{c,1.705,2.11},	
		    PlotRange->{{1.705,2.11},{-8.225,-7.96}},
		    PlotLabel->"(a)  Weibull distribution\n         z  =  Q(c,u(c))",
		    AspectRatio->2/GoldenRatio,
		    AxesOrigin->{1.685,-8.235},
			AxesLabel->{"c","z"}];
right = Plot[WG[c],{c,3.21,4.13},	
		    PlotRange->{{3.21,4.13},{-6.89,-6.785}},
		    PlotLabel->"(b)  Gamma distribution\n          z  =  q(c,w(c))",
		    AspectRatio->2/GoldenRatio,
		    AxesOrigin->{3.16,-6.895},
			AxesLabel->{"c","z"}];

Show[GraphicsArray[{{left,right}}],PlotLabel->
   "Figure 31.5  Functions whose maximizers yield shape parameter MLEs"]