Professor of Stochastic Modelling | |
School of Mathematics & Statistics | |
Newcastle University |
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
?statplotsYou don't need to read this now - we will come back to it later.
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?
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).
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).
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!
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.
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.
darren.wilkinson@ncl.ac.uk | ||
http://www.staff.ncl.ac.uk/d.j.wilkinson/ |