ClearAll;
$DefaultFont = {"Times-Roman",10};
Needs["Statistics`DataManipulation`"];
Needs["Graphics`Graphics`"];
njags = 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 = datapoints+Join[{{0,0}},Table[{1/120,0},{k,1,15}]];
density = {0,0,0,
			  9/489,5/489,28/489,45/489,
			  165/489,90/489,118/489,17/489,
			  9/489,0,1/489,2/489};

n = 7;
F[a_,x_]:= 1 - E^(-(a x)^n);
f[a_,x_]:= n a (a x)^(n-1) E^(-(a x)^n);

error[a_] = Apply[Plus,Table[(F[a,(2k+1)/120]-cumulative[[k]])^2,{k,1,15}]];
leasterror = FindMinimum[error[a],{a,7,0,10}];
aopt = a/.leasterror[[2]];
g[x_]:= f[aopt,x];

poly[a_,b_]:= Join[Table[{a+i(b-a)/njags,g[a+i(b-a)/njags]},
				  {i,0,njags}],{{b,0},{a,0}}];
shade[a_,b_,c_]:= Graphics[{GrayLevel[c],Polygon[poly[a,b]]}];
	 


pdf:= Plot[g[x],{x,0,0.25}];
plot = Show[shade[0.09,0.12,0.65],shade[0.12,0.15,0.35],shade[0.15,0.18,0.65],pdf,
			PlotLabel->"Distribution of leaf thickness in Dicerandra linearifolia",
			GridLines->Automatic,PlotRange->{{0,0.25},{0,17.5}},
			Axes->True,AxesLabel->{" x","y"}];
			
Show[GraphicsArray[{{plot}}],PlotLabel->"Figure 8.4"]