Prof Darren Wilkinson Professor of Stochastic Modelling

School of Mathematics & Statistics

Newcastle University

Basic probability distributions with Maple

This computer practical provides an introduction to the use of Maple for doing plots and calculations relating to simple probability distributions. You should work through the practical doing the exercises as you go along. Don't try and split the exercises from the rest of the practical - just do them as you get to them in the same worksheet, and then just print off the whole lot when you are finished. When you have completed the practical, you must follow the instructions given at the end in order to get your mark for this session.

Getting started

Load up Maple, and type the following at the Maple command prompt ( > ).

with(stats);

This command loads up the Maple stats library package, which adds some functions to Maple which are useful for probability and statistics. All Maple commands must end with a colon (:) or a semi-colon (;). A colon is used when we do not want to see the results of a command (this will be useful later).

Having this practical on-line means that you can copy-and-paste commands from the web browser directly into Maple. Practice this now. First highlight the next command

statplots[histogram]([random[normald](100)]);

with the mouse, then select Edit - Copy in the web browser. Now switch back to Maple, and select Edit - Paste, and then press Return. This should help save lots of typing - and typos! You can probably guess what this command does, but if you can't don't worry about it - all will be explained later. If there is any part of this practical you don't understand, you can get help on the commands used using the Maple on-line help. Simply type the command you want help for at the Maple prompt preceded by a ?. Try this now by typing

?statplots

You don't need to read this now - we will come back to it later.

The Poisson distribution

The Maple function statevalf knows all kinds of things about all of the standard probability distributions. However, it has lots of options, and so is a bit messy to use directly. What we normally do is define our own simple functions using it. This is easiest to understand in the context of an example. Let's define a function called poispmf representing the PMF of the Poisson distribution. We want this to have two arguments: x, the value of interest, and lambda, the parameter of the distribution. We can declare this as follows.

poispmf:=(x,lambda)->statevalf[pf,poisson[lambda]](x);

Here pf stands for discrete Probability mass Function. You can find out more about statevalf by entering ?statevalf. Do this now, and make sure you understand basically how it works. So, now suppose we want to know the probability of X=0 when X is Poisson with mean 5, we can enter

poispmf(0,5);

and we will get the answer. This is just a shorthand for typing statevalf[pf,poisson[5]](0), but poispmf(0,5) is much easier, so it's always worth declaring functions for complex expressions.

Unfortunately, Maple isn't very good at plotting discrete distributions, but we can do it using the following command - see if you can figure out how it works.

plot([seq([i,poispmf(i,5)],i=0..20)],style=point);

If we wanted to know the probability that X is strictly less than 3, we could do this by brute force as follows.

sum(poispmf(r,5),r=0..2);

Fortunately, statevalf also knows about Discrete Cumulative Distribution Functions. We can define a CDF function, evaluate it at 2 and plot it as follows.

poiscdf:=(x,lambda)->statevalf[dcdf,poisson[lambda]](x);
poiscdf(2,5);
plot([seq([i,poiscdf(i,5)],i=0..20)],style=point);

Now, we can go back to the PMF, and verify that it sums up to one.

evalf(sum(poispmf(r,5),r=0..infinity));

We can also calculate the expected value by getting Maple to do the relevant sum for us.

evalf(sum(r*poispmf(r,5),r=0..infinity));

Exercises:

1. If X is Poisson with mean 8, what is the probability that X=10?
2. What is the probability that 10 < X < 20?

The exponential distribution

Maple can be an extremely useful tool for all sorts of computations relating to continuous distributions. Let's start by looking at the PDF of the exponential distribution. The exponential distribution is a family of continuous probability distributions that is indexed by a single parameter, lambda, that changes the scale of the distribution.

exppdf:=(x,lambda)->statevalf[pdf,exponential[lambda]](x);
plot(exppdf(x,1),x=-1..5);

Here pdf represents a continuous Probability Density Function. This gives us a plot for the case lambda=1, but we can put multiple graphs on the same plot if we want:

plot([exppdf(x,1),exppdf(x,2)],x=-1..5);

If we are feeling really ambitious, we can load up the plots library package, and look at an animation for how the PDF changes as the exponential parameter varies.

with(plots):
animate(exppdf(x,lambda),x=0..2,lambda=1..10,frames=50);

When you first type this command, it will just bring up a plot. To get the animation to play, click on it with the Right mouse button and then select Animation->Play. Of course, we can also plot the Cumulative Distribution Function as follows.

expcdf:=(x,lambda)->statevalf[cdf,exponential[lambda]](x);
plot(expcdf(x,2),x=-1..10);

So, for X distributed exponentially with parameter 2, we can read cumulative probabilities off the above graph. In particular, the probability that X is less than 1 is given by

expcdf(1,2);

We can verify this by integrating the PDF as follows.

evalf(int(exppdf(x,2),x=0..1));

Indeed, we can verify that the PDF integrates to 1 with

evalf(int(exppdf(x,2),x=0..infinity));

We can even get Maple to calculate the mean and variance of X using

mu:=evalf(int(x*exppdf(x,2),x=0..infinity));
evalf(int(((x-mu)^2)*exppdf(x,2),x=0..infinity));

Exercises:

3. Suppose that X is exponential with parameter 0.2. Calculate P(3<X<6) both by using the CDF, and by integrating the PDF.
4. Also calculate E(X).

The normal distribution

The normal distribution is probably the most important distribution in the whole of probability theory - you will learn more about it in Semester 2. It has two parameters (mu and sigma), which represent its mean and standard deviation, but this doesn't really change things much. Note that Maple refers to the normal distribution as normald - NOT normal, which in Maple means something else entirely.

normpdf:=(x,mu,sigma)->statevalf[pdf,normald[mu,sigma]](x);
plot([normpdf(x,0,1),normpdf(x,0.5,0.3)],x=-3..3);

This defines a normal PDF, plots the density for a couple of different parameter values. To understand how the distribution changes for different parameter values, consider the following.

animate(normpdf(x,m,2),x=-10..10,m=-5..5,frames=50);
animate(normpdf(x,0,sigma),x=-10..10,sigma=1..5,frames=50);

This produces an animation for the PDF as mu varies, and another as sigma varies (the animate command requires the plots package to be loaded, but we did this earlier). Follow the instructions from earlier in order to get the animations to play.

Now let's look at the CDF.

normcdf:=(x,mu,sigma)->statevalf[cdf,normald[mu,sigma]](x);
plot(normcdf(x,0,1),x=-4..4);

This allows us to read off cumulative probabilities at a given x-value. However, we often want to do the reverse; that is, read off the x value corresponding to a given probability. We could get Maple to solve this problem for us by brute force. For example, to find the 95% point of the standard normal, we could use

solve(normcdf(x,0,1)=0.95,x);

Fortunately we don't need to resort to this for standard distributions, as Maple also knows about Inverse CDFs, which we can define and plot as follows.

normicdf:=(p,mu,sigma)->statevalf[icdf,normald[mu,sigma]](p);
plot(normicdf(p,0,1),p=0..1);

We can now use this function to find the 5% and 95% points of the standard normal.

normicdf(0.05,0,1);
normicdf(0.95,0,1);

Exercises:

5. Suppose that X is normal with mean 4 and standard deviation 5. Verify that the mean and median (50% point) are both 4, and that the variance is 25.
6. Calculate the interquartile range (IQR).

Non-standard distributions

Maple is very happy to work with any PDF you care to use. For example, the following function defines a density on (0,infinity).

pdf:=(x)->x*exp(-x);
plot(pdf(x),x=0..8);

In fact, this is a very standard distribution, which you'll find out more about in MAS271 next year.

Exercises:

7. Using the PDF above, use Maple to verify that it is a PDF by integrating it over (0,infinity).
8. Calculate the mean, median and variance of the distribution.
9. What are the 5% and 95% points of the distribution?
10. What is the probability that a value from this distribution will lie between 2 and 4?

This distribution is fairly straightforward to work with analytically, so you can check your answers when you are finished by re-doing some of the above calculations by hand!

Randomly generated data

We can also get Maple to generate (pseudo) random data from any of the distributions it knows about. We can generate 1000 observations from a Normal with mean 5 and standard deviation 2 as follows.

normdat:=[random[normald[5,2]](1000)]:

Note that by ending the command with a colon rather than a semi-colon, we do not have to see all 1000 observations on the screen. Now that we have our data, we can use statplots to produce some graphical summaries.

statplots[histogram](normdat);
statplots[boxplot](normdat);

Alternatively, we can use describe to obtain all of the usual descriptive statistics.

describe[mean](normdat);
describe[range](normdat);
describe[quartile[1]](normdat);
describe[quartile[2]](normdat);
describe[quartile[3]](normdat);
describe[standarddeviation](normdat);

Exercises:

11. Generate 1000 observations from an exponential distribution with parameter 0.3.
12. Plot a histogram of these values, and calculate their upper quartile.
13. Check that this is close to the theoretical 75% point of the distribution, which you can obtain from the inverse CDF.

So what?

Obviously, this short practical has only illustrated some of the functionality provided by the Maple stats library. It has barely scratched the surface of possibilities for using Maple to assist probabilistic and statistical computation. However, by combining the knowledge of the stats library gained in this practical together with more general knowledge of the Maple language gained from the MAS102 Maple labs, you should be able to begin to use Maple to help you in your probability and statistics courses.

You should probably "bookmark" this page (or at least remember how to get to it), as it will prove very useful in future months and years for reminding you of the basic syntax for some of the more commonly used stats library commands. Also, when you do bivariate distributions in MAS271, you may want to come back here and look at my page on the bivariate normal distribution.

When you are finished...

When you are finished you should print out your worksheet. Then in a pen which is not black (preferably a brightly coloured felt-tip), write the exercise numbers in the margin next to the commands you used to solve each exercise (eg. write a "1" in the margin next to the command you used to solve exercise 1). You should staple this to the back of your solutions to the last lot of exercises, and hand it in with that on the 13th of December (or before).

Darren J Wilkinson

darren.wilkinson@ncl.ac.uk

http://www.staff.ncl.ac.uk/d.j.wilkinson/

	Professor of Stochastic Modelling
	School of Mathematics & Statistics
	Newcastle University


		`darren.wilkinson@ncl.ac.uk`
		`http://www.staff.ncl.ac.uk/d.j.wilkinson/`