ClearAll; Clear[f,sf];
$DefaultFont = {"Times-Roman",10};
Needs["Graphics`Graphics`"];
njags = 100;
mean[c_]:= Gamma[1+1/c];
mu2 = mean[2];
mu3 = mean[3];
mu10 = mean[10];
poly[a_,b_,f_]:= Join[Table[{a+i(b-a)/njags,f[a+i(b-a)/njags]},
				  {i,0,njags}],{{b,0},{a,0}}];
shade[a_,b_,c_,f_]:= Graphics[{GrayLevel[c],Polygon[poly[a,b,f]]}];
f[m_,z_]:= m z^(m-1) E^(-(z)^m);
f2[z_]:= f[2,z]; f3[z_]:= f[3,z]; f10[z_]:= f[10,z];
D2[z_]:= (z - mu2)^2 f2[z]; 
D3[z_]:= (z - mu3)^2 f3[z]; 
D10[z_]:= (z - mu10)^2 f10[z];

shading10 = Show[shade[0,1.5,0.65,D10],
				Axes->True,
				Ticks->{Automatic,Table[0.01i,{i,0,5}]},
				PlotRange->{{0,1.5},{0,0.045}},
				PlotLabel->"c = 10",
				AxesLabel->{"x/s","s D"}];

plotu10:= Plot[f10[z],{z,0.001,3},PlotLabel->"c = 10",
		Ticks->{Automatic,Table[i,{i,0,5}]},
		AxesLabel->{"x/s","s f"},PlotRange->{{0,1.5},{0,3.7}}];
plotl10:= Plot[D10[z],{z,0.001,1.5}];

shading3 = Show[shade[0,3,0.65,D3],
				Axes->True,
				Ticks->{Automatic,Table[0.05i,{i,0,5}]},
				PlotRange->{{0,3},{0,0.12}},
				PlotLabel->"c = 3",
				AxesLabel->{"x/s","s D"}];

plotu3:= Plot[f3[z],{z,0.001,3},PlotLabel->"c = 3",
		Ticks->{Automatic,Table[0.2i,{i,0,5}]},
		AxesLabel->{"x/s","s f"}];
plotl3:= Plot[D3[z],{z,0.001,3}];

shading2 = Show[shade[0,3,0.65,D2],
				Axes->True,
				Ticks->{Automatic,Table[0.1i,{i,0,5}]},
				PlotRange->{{0,3},{0,0.21}},
				PlotLabel->"c = 2",
				AxesLabel->{"x/s","s D"}];

plotu2:= Plot[f2[z],{z,0.001,3},PlotLabel->"c = 2",
		Ticks->{Automatic,Table[0.2i,{i,0,5}]},
		AxesLabel->{"x/s","s f"}];
plotl2:= Plot[D2[z],{z,0.001,3}];

upper10 = Show[plotu10];
lower10 = Show[shading10,plotl10];
upper3 = Show[plotu3];
lower3 = Show[shading3,plotl3];	
upper2 = Show[plotu2];
lower2 = Show[shading2,plotl2];
Show[GraphicsArray[Transpose[{{upper10,upper3,upper2},{lower10,lower3,lower2}}]],
		PlotLabel->"                                         Figure 27.1\n
		Weibull variance increases as shape parameter decreases"]