UC SCI 111, vM+ 
Notebook 6

<< Statistics`DataManipulation`
<<Statistics`ContinuousDistributions`
<< Graphics`Graphics`
<< 
      Graphics`Colors`
<< 
      Statistics`DiscreteDistributions`
<< 
      Statistics`DescriptiveStatistics`
Off[
        General::spell1];

How to let Mathematica do your statistics.

Conducting an experiment: RandomArray, Random.

Various samples can be obtained using the Random function. To produce 100 random integers between 30 and 40 one types:

Table[Random[Integer, {30, 40}], {100}]

Exercise 3. 
Produce 100 random real numbers between 3 and 4.

The distribution underlying the sample you obtained above is the uniform distribution. Any distribution can be obtained from the uniform distribution and there are two Mathematica packages that help perform this task for you. Discrete distributions are in  Statistics`DiscreteDistributions`  and  continuous distributions are in  Statistics`ContinuousDistributions`. With these packages it is easy to produce samples.

ndist = NormalDistribution[0, 5]

ndist is now nothing but a symbol that represents a normal distribution with mean 0 and standard deviation 5. It can be used with the RandomArray function.

RandomArray[ndist, 100]

?RandomArray

Exercise 4.
a) Produce a sample of size 100 with the Poisson distribution. Set λ=2.
b) Write a function

	MyRandomArray[distribution,n]
	
in Mathematica, showing that there is nothing mysterious about the command RandomArray.
Hint: read help about Random.

The definition Mathematica makes for the geometric distribution is different from the one in the lecture notes, as we can see from the following example:

gdist = GeometricDistribution[0.4];
sample = RandomArray[gdist, 200];

Mean[sample] // N

Mean[gdist]

So we get a mean for the sample here and a theoretical mean, equal to 1.5, according to Mathematica. In the lecture notes it says that the mean equals 1/p . In this case 1/p equals 2.5. So there is a difference of 1. To produce a sample with the geometric distribution as in the lecture notes, add 1 to all elements of the sample set. So:

gdist = GeometricDistribution[0.4];
sample = RandomArray[gdist, 2000] + 1;

Mean[sample] // N

You see that the last outcome corresponds better to the theoretical outcome predicted in the lecture notes.

Exercise 5.
With the definition of the lecture notes the standard deviation of a geometric distribution is

	((Sqrt[1 - p])/p)

See also Exercise 1. With Mathematica's definition of the geometric distribution, what is the standard deviation?

The empirical cumulative distribution function.

If we have a sample s, that is a subset of ℝ, we can define the empirical cumulative distribution function relative to the sample. For every x ∈ ℝ it is defined as the fraction of elements of the sample s that is smaller than or equal to x. So if the sample is

sample = {1, 7, 0, 1, 2, 0, 0, 12, 1, 2, 1, 0, 0, 4, 0, 1, 0, 1, 3, 
          4, 2, 0, 0, 2, 3, 4, 0, 1, 1, 0, 2, 3, 0, 0, 2, 0, 0, 3, 8, 1, 1, 
          2, 0, 0, 1, 2, 2, 2, 8, 2, 1, 2, 0, 1, 0, 1, 0, 5, 3, 0, 0, 1, 5, 
          3, 4, 1, 1, 1, 4, 4, 3, 0, 2, 1, 3, 0, 0, 2, 0, 0, 2, 0, 0, 0, 3, 
          2, 0, 3, 1, 0, 1, 0, 1, 3, 0, 0, 0, 3, 2, 0, 0, 1, 0, 1, 0, 0, 1, 
          0, 2, 4, 0, 0, 1, 2, 0, 1, 1, 3, 0, 0, 1, 1, 2, 0, 0, 1, 0, 0, 0, 
          0, 1, 0, 1, 3, 2, 1, 4, 0, 2, 4, 0, 0, 0, 1, 0, 0, 1, 0, 0, 5, 4, 
          1, 0, 3, 2, 0, 0, 0, 0, 4, 2, 5, 1, 0, 1, 1, 7, 0, 3, 0, 1, 0, 0, 
          8, 0, 0, 0, 2, 1, 0, 1, 1, 1, 2, 0, 0, 1, 0, 2, 3, 2, 0, 0, 0, 2, 
          0, 1, 1, 1, 8};

We can find ECDF(3.5) by looking at the number of elements that are smaller then or equal to 3.5.

count = 0;
Do[If[sample[[i]] <= 3.5, count = count + 1], {i, 1, 
        Length[sample]}]

count

Length[sample]

So 178 out of 200 elements are smaller than or equal to 3.5. Therefore ECDF(3.5) is (178/200) = 0.89.  As it would be quite a hassle to calculate the ECDF in the way we just did when you want to plot it, here is a function that calculates the ECDF of a sample directly.

ECDF[sample_][x_] := 
      (Plus@@(If[# <= x, 1, 0]&/@sample))/Length[sample];

Example:

ECDF[sample][3.5]

A plot is just as easy:

Plot[ECDF[sample][x], {x, 0, 7}];

The sample was produced with the geometric distribution with parameter 0.4.

Exercise 6.
What is the theoretical prediction of ECDF[3.5]?

As you saw in the previous exercise the ECDF corresponds to the CDF. The larger the sample, the closer to the CDF the ECDF will be. So one way to check whether a sample corresponds to a certain random variable is to compare the CDF of the random variable with the ECDF of the sample.

\(
\(ndist = NormalDistribution[0, 5];\)\)
sample = RandomArray[ndist, 100];
Plot[{ECDF[sample][x], CDF[ndist, x]}, {x, -7, 7}, 
        PlotStyle -> {Red, Green}];

Here the ECDF is plotted in Red, and the CDF is plotted in Green.

Other comparisons: mean, standard deviation and the histogram.

One can also compare numerical parameters associated to both a sample and a random variable. If these parameters are almost equal one can at least conclude that the random variable might underlie a certain sample. It is especially convenient to compare a set of parameters that completely determine the random variable. Thus, good candidates are the mean and the standard deviation. Take the last sample:

ndist = NormalDistribution[0, 5];
sample = RandomArray[ndist, 100];

Suppose we didn't know that the sample came from a normal distribution with mean 0 and standard deviation 5, but that we did suspect already that it came from a normal distribution. We could do the following

mn = Mean[sample]

sd = StandardDeviation[sample]

The normal distribution with these parameters indeed fits better than the one we used to produce our sample. That is when we plot the ECDF with the CDF corresponding to mn and sd, we find that they fit better according to our eyes.

ndist2 = NormalDistribution[mn, sd];
Plot[{ECDF[sample][x], CDF[ndist2, x]}, {x, -7, 7}, 
        PlotStyle -> {Red, Green}];

Plot[{ECDF[sample][x], CDF[ndist, x]}, {x, -7, 7}, 
        PlotStyle -> {Red, Green}];

The question whether a certain sample has some underlying random variable is the subject of statistical tests. This is a subject which we shall not pursue in these notes.

Another efficient comparison is the comparison between the (theoretical) Probability Density Function (abbreviated in the sequel as PDF) and the histogram of a sample. To produce a histogram, don't use the built-in facilities of mathematica but instead the following routine. The code of this is somewhat involved, and we do not explain the theoretical considerations leading to it, so please don't try to understand it....

\(\(histogram[sample_List, \ width_?NumberQ, opts___]\  := \ 
        Module[
	{s, n, min, \ max, \ relfreq, maxfreq, \ \ mids, \ 
            hbasis, \ dt, \ i, dps, opts2}, 
		s\  = \ 
            Sort[sample]; 
		n\  = \ Length[s]; 
		min\  = \ 
            s\ [\([1]\)]; 
		max\  = \ 
            s[\([n]\)]; 
		dps = \(DisplayFunction /. {opts}\) /. 
              Options[histogram]; 
		relfreq\  = \ 
            BinCounts[\ 
                sample, \ {min, \ max, \ width}]/\((n\ *
                  width)\)\ ; 
		maxfreq\  = \((9/8)\)*\ 
              Max[relfreq]; 
		mids = \ 
            N[Table[\ 
                min + \ width \((\  .5\  + \ i)\), \ {i, \ 0, \ 
                  Length[relfreq]\  - \ 1}], 
              2]; 
		hbasis\  = \ \((2/5)\)\ width; 
		dt\  = \ 
            Table[
			Rectangle[
				{mids[\([i]\)] - hbasis, 
                  0}, \ {mids[\([i]\)] + hbasis, \ 
                  relfreq[\([i]\)]}
				], 
			{i, \ 
                Length[mids]}
		]; 
		If[{opts} === {}, 
            opts2 = Rule[DisplayFunction, $DisplayFunction], 
            opts2 = opts]; 
		Show[
            Graphics[
			Join[{Red}, 
                dt]], 
			opts, \ 
			Axes\  -> \ 
              True, 
			AxesOrigin\  -> \ {Floor[min], 0}, \ 
            PlotRange\  -> \ {0, maxfreq}
		]
];\)
\)
\(\(histogram[sample_List, opts___]\  := \ 
        Module[
		{n, \ width, s}, 
		n\  = \ Length[sample]; 
		s = 
            Sort[sample]; 
		width = 
            N[2*\((s[\([Ceiling[3  n/4]]\)] - 
                      s[\([Ceiling[n/4]]\)])\)/\((n^\((1/
                          3)\))\)]; 
		histogram[sample, width, 
            opts]];\)	\)

.. just produces a nice histogram.

gr1 = histogram[sample, DisplayFunction -> Identity];
gr2 = 
        Plot[PDF[ndist, x] // Evaluate, {x, -10, 10}, 
          PlotStyle -> {Green}, DisplayFunction -> Identity];
Show[{gr1, gr2}, DisplayFunction -> $DisplayFunction];

If we increase the size of the sample, the graph of the PDF fits the histogram much better.

ndist = NormalDistribution[0, 5];
sample = 
        RandomArray[ndist, 
          1000];

gr1 = histogram[sample, DisplayFunction -> Identity];
gr2 = 
        Plot[PDF[ndist, x] // Evaluate, {x, -10, 10}, 
          PlotStyle -> {Green}, DisplayFunction -> Identity];
Show[{gr1, gr2}, DisplayFunction -> $DisplayFunction];

Projects.

A sample.

Exercise 7.
Suppose that a random variableY has a PDF (e^(-x)). The sample space of Y is the interval [0,∞).
a) What is the distribution function of Y ?
b) Show how to apply Exercise 2) in order to obtain random data corresponding with the distribution of Y.
c) Produce a sample of 500 drawings of this random variable.
Hint: Use Random[] to obtain a sample from the uniform distribution on [0,1]. 
d) Plot the density function of Y together with a histogram of the data you produced.

The central limit theorem.

The aim of this project is to illustrate the central limit theorem.  
Let (Subscript[X, 1])and (Subscript[X, 2]) be two random variables with distributions ndist1 and ndist2.

ndist1 = NormalDistribution[0, 5];
ndist2 = NormalDistribution[20, 5];

Exercise 8.
a) Create a sample of size 1000 based on the random variable (((Subscript[X, 1]) + (Subscript[X, 2]))/2) . 
b) Suppose (Subscript[X, 1])has (Subscript[f, 1])as its PDF and (Subscript[X, 2]) (Subscript[f, 2]). Does  (((Subscript[X, 1]) + (Subscript[X, 2]))/2) have  (((Subscript[f, 1]) + (Subscript[f, 2]))/2) as its PDF? Illustrate your answer with a picture based on the sample you obtained in a).
☺ c) Calculate the PDF of  (((Subscript[X, 1]) + (Subscript[X, 2]))/2).

Let
 
	(Subscript[S, n])  =  (1/n) (∑_(i = 1)^n)(Subscript[X, i]), 

for  i.i.d. random variables (Subscript[X, i]) with mean μ and standard deviation σ.   Let 

	\(TraditionalForm\`S\&~\_n = \(\@n\) \((\(S\_n - μ\)\/σ)\)\)

Take for the (Subscript[X, i]) a Poisson distribution with parameter 1:

pdist = PoissonDistribution[1]

d) Generate samples of size 500 for the \(TraditionalForm\`\(S\_i\)\&~\), i=1,...10.
e) Compute their means and standard deviations.
f) Plot the ECDF of these samples.
g) What do you conclude about the applicability of the central limit theorem?

☺ Advanced Graphics.

Sometimes one wants to combine several graphical objects in one picture. To learn you how to do this we provide some background information in this appendix. 

 Plot is a Mathematica function, i.e. : you put something in Plot - the arguments - and you get something back - the result. Let's see how this works.

\(f[x_] := x^2; \ 
gr1 = Plot[f[x], {x, \(-2\), 2}]; \)

The variable gr1 now contains the result that we are looking for: but in what form?

FullForm[gr1]

We see that gr1 is a big object, but it starts with Graphics. We can check this by asking for the Head of the expression.

Head[gr1]

Everything included inside the Graphics is a big list of so-called graphics primitives (Line, Circle, Point, Disk, Rectangle etc..)  and options for the Graphics, such as ImageSize, Axes, AspectRatio. Thus, it is perfectly possible to use gr1 to reproduce a modified version of the plot we made above:

Show[gr1, ImageSize -> {100, 100}, Background -> VenetianRed];

Suppose that we wanted in addition to have a Banana-colored disk in the picture. And also some nice Tomato-colored rays. We could type the following:

\(n = 20; 
p1 = {1, 3}; 

    gr2 = Graphics[{
          Join[{Tomato, Thickness[0.02]}, 
            Table[r1 = 1; 
	r2 = 2; 
	
              Line[{p1 + 
                    r1*{Cos[\((\((k*Pi*2)\)/n)\)], 
                        Sin[\((\((k*Pi*2)\)/n)\)]}, 
                  p1 + r2*{Cos[\((\((k*Pi*2)\)/n)\)], 
                        Sin[\((\((k*Pi*2)\)/n)\)]}}], 
	{k, 1, n}]], {
            Banana, Thickness[0.02], Disk[p1 + {0, 0}, r1]}}]; \)

You can look at FullForm[gr2] to see that it is a similar to FullForm[gr1]. They are both valid Graphics objects and so we can show them in one picture.

Show[{gr2, gr1}, Axes -> True, AspectRatio -> Automatic, 
      ImageSize -> {100, 100}];

In this way we can combine any number of graphics objects. 

There are also a lot of other ways to manipulate Graphics. Some of what can be done is in the demos section of the help files. Enthusiasts tend to run into Beethoven  there. The picture of Beethoven is not the result of some mathematical formula, it was created with a CAD/CAM program, and then imported into Mathematica.

We have already seen some options that can be passed to Graphics objects: AspectRatio, Axes, ... Have a look in the help files if you want to know exactly what they mean. Here we will explain one option that is particularly useful: the DisplayFunction option. This controls whether graphics are going to be displayed or not. The two values that you'll use are Identity and $DisplayFunction. If you write

gr1 = 
        Plot[Sin[x], {x, -Pi, Pi}, 
          DisplayFunction -> Identity];

you will not see very much. But, you can force Mathematica to show gr1 by typing

Show[gr1, DisplayFunction -> $DisplayFunction, 
      ImageSize -> {200, 100}];

In the notebook we repeatedly suppress the display of graphs until we have gathered all graphs and then conclude by showing all of them in one picture. Click here to go to an example. 

Some of you who have already been experimenting with Mathematica will remark that Plots and ParametricPlots can be combined in a much more simple fashion, i.e. you don't have to write

gr1 = Plot[x, {x, 0, 1}, DisplayFunction -> Identity];
gr2 = Plot[x^2, {x, 0, 1}, DisplayFunction -> Identity];
gr3 = Plot[x^3, {x, 0, 1}, DisplayFunction -> Identity];
gr4 = Plot[x^4, {x, 0, 1}, DisplayFunction -> Identity];
gr5 = Plot[x^5, {x, 0, 1}, DisplayFunction -> Identity];
Show[{gr1, gr2, gr3, gr4, gr5}, DisplayFunction -> $DisplayFunction, 
        ImageSize -> {100, 100}, Axes -> False];

, since

Plot[{x, x^2, x^3, x^4, x^5}, {x, 0, 1}, ImageSize -> {100, 100}, 
      Axes -> False];

will do just as well. The first method, though, can also combine one Plot and one ParametricPlot.