Elliot Tanis
June 9, 2006
Expansion of Example 5.3-2  (Example_5_3-2.mws)

A. We first investigate the behaviour of observations of a Cauchy random variable.

>	randomize(): m := 500: X := [seq(evalf(tan(Pi*(rng() - 1/2))), k=1..m)]: minimum := Min(X); median := Median(X); maximum := Max(X);

B. Look at the smallest and largest observations.

>	S := sort(X): S[1..5]; S[m-4..m];

C. Construct a running average plot.

>	display({PlotRunningAverage(X)}, title=`Running Average of Cauchy Observations`, titlefont=[TIMES,BOLD,14]);

D. Look at a box-and-whisker diagram and a histogram that has 8 equal length intervals.

>	display({BoxPlot(X)}, title=`A Not Very Meaningful Box Plot\nof Cauchy Observations`, titlefont=[TIMES,ROMAN,14]);

>	display({HistogramFill(X)}, title=`A Not Very Meaningful Histogram\nof Cauchy Observations`, titlefont=[TIMES,ROMAN,14]);;

E. In order to make some sense out of these data graphically, replace all observations that are less than -5 with -5 and all observations that are greater than 5 with a 5.

>	k := 'k': xtics := [k$k=-5..5]:

>	Y := [seq(max(min(X[k], 5), -5), k = 1 .. m)]: P1 := BoxPlot(Y): P2 := plot([[0,0],[0,0]], x = -5.5 .. 5.5, y = 0 .. 1.2, xtickmarks=xtics, ytickmarks=[], labels=[``,``]): display({P1,P2}, title=`Box Plot of Modified Cauchy Observations`, titlefont=[TIMES,ROMAN,14]);

>	ytics := [0.05*k$k=1..6]:

>	P1 := HistogramFill(Y,-5..5,20): f := 1/Pi/(1 + x^2): P2 := plot(f, x=-5..5, y = 0 .. 0.35, color=blue, thickness=2, xtickmarks=xtics, ytickmarks=ytics, labels=[``,``]): display({P1,P2},title = `Modified Cauchy Observations`, titlefont=[TIMES,ROMAN,14]);

>	ytics := [0.1*k$k=1..10]:

>

P3 := EmpCDF(Y,-5..5): F := (arctan(x) + Pi/2)/Pi: P4 := plot(F,x=-5.4 ..5.4, xtickmarks=xtics, ytickmarks=ytics, color=blue, thickness=2, labels = [``,``]): display({P3, P4},title = `Empirical Distribution Function of Modified\nCauchy Observations, Cauchy Distribution Function`, titlefont=[TIMES,ROMAN,14]);

F. Now investigate the behaviour of x-bar.

>	randomize(): n := 5: m := 500:  #Try other values of n for j from 1 to m do   XX := [seq(evalf(tan(Pi*(rng() - 1/2))), k=1..n)]:   xbar[j] := Mean(XX): od: X := [seq(xbar[j], j=1..m)]: minimum := Min(X); median := Median(X); maximum := Max(X); S := sort(X): S[1..5]; S[m-4..m];

>	display({PlotRunningAverage(X)}, title=`Running Average of Cauchy Means`, titlefont=[TIMES,ROMAN,14]);

G. Make a box-and-whisker diagram and a histogram of the modified data.

>	Y := [seq(max(min(xbar[k], 5), -5), k = 1 .. m)]: P1 := BoxPlot(Y): P2 := plot([[0,0],[0,0]], x = -5.5 .. 5.5, y = 0 .. 1.2, xtickmarks=xtics, ytickmarks=[]): display({P1,P2}, title=`Box Plot of Modified Cauchy Means`, titlefont=[TIMES,ROMAN,14]);

>

P1 := HistogramFill(Y,-5..5,20): f := 1/Pi/(1 + x^2): P2 := plot(f, x=-5..5, y = 0 .. 0.35, color=blue, thickness=2, xtickmarks=xtics, ytickmarks=ytics, labels=[``,``]): tn := textplot([2.5, 0.25, `n = `||n], font=[TIMES,ROMAN,14]): display({P1,P2,tn},title = `Histogram of Modified Cauchy Means`, titlefont=[TIMES,BOLD,14]);

H. Now find the p.d.f. of x-bar theoretically. (This solution was written by John Krueger when he was a Hope College junior.)

>	p := 'p': q := 'q': t := 't': assume(p > 0): additionally(q > 0):

Let   and    have independent Cauchy distributions and let  =   be their   weighted average where  p  > 0, q  > 0. Then  

  will have the p.d.f.   *CauchyPDF( x *( p + q )/ p )),  

   will have the p.d.f.   * CauchyPDF ( x *( p  + q)/ q )),

and the p.d.f. of   is just the convolution of these two p.d.f.'s which is:

(REMARK: The integral should go from      to     and this worked in Release 3. It does not always work in Releases 4,  5, 6...)

>	#simplify(int(((p + q)/p*CauchyPDF(t*(p + q)/p)) * ((p + q)/q*CauchyPDF((x - t)*(p + q)/q)), t = -infinity .. infinity));

>	simplify(int(((p + q)/p*CauchyPDF(t*(p + q)/p)) * ((p + q)/q*CauchyPDF((x - t)*(p + q)/q)), t = -infinity .. 0) + int(((p + q)/p*CauchyPDF(t*(p + q)/p)) * ((p + q)/q*CauchyPDF((x - t)*(p + q)/q)), t = 0 .. infinity));

But this is just the Cauchy p.d.f.

>	CauchyPDF(x);

Now, if p = q  = 1 , the weighted average just comes out to be the mean of the two random variables,    and  , which means that the mean of two independently distributed Cauchy random variables has a Cauchy distribution. By fixing q  at 1 and using induction on p, it follows that the mean of any number of independently distributed Cauchy random variables has a Cauchy distribution.

We shall now animate the p.d.f. of X -bar when sampling from a Cauchy distribution and compare this p.d.f. with that of the standard normal p.d.f., N (0, 1).

>	k := 'k': xtics := [k$k=-5..5]: ytics := [0.1*k$k=1..4]:

>

normplot := plot(NormalPDF(0, 1, x), x = -5.1 .. 5.1, y = 0 .. 0.42, xtickmarks=xtics, ytickmarks=ytics, color=blue, thickness=2, labels=[``,``]): colors := [aquamarine, blue, navy, coral, cyan, brown, gold, green, gray,   khaki, magenta, maroon, orange, pink, plum, red, sienna, tan, turquoise, violet, wheat, yellow]:  for n from 1 to 22 do   P[n] := plot(CauchyPDF(x), x = -5.1 .. 5.1, y = 0 .. 0.42, xtickmarks=xtics, ytickmarks=ytics, color=colors[n], thickness=2, labels=[``,``]):   tn := textplot([2, 0.3, `n = `||n], font=[TIMES,ROMAN,14]):   Pn[n] := display({P[n],tn,normplot}): od: display([seq(Pn[j], j = 1 .. 22)], title = `N(0, 1) p.d.f. and the p.d.f. of X-bar\nWhen Sampling from a Cauchy Distribution`, insequence=true, titlefont=[TIMES,ROMAN,14]);

>

Click on the above figure and use the tool bar to animate it.