|
|
Professor of Stochastic Modelling |
| School of Mathematics & Statistics | |
| Newcastle University |
f:=(x,y,rho)->(1/(2*Pi*sqrt(1-rho^2)))*exp((-1/(2*(1-rho^2)))*(x^2+y^2-2*rho*x*y));This density has both means equal to zero and both variances equal to 1. However, it does have a parameter, rho, which determines the correlation between the X and Y components.
Let's start by plotting the density for a particular value of the correlation parameter - say 0.8. This corresponds to fairly high positive correlation between the components.
plot3d(f(x,y,0.8),x=-4..4,y=-4..4,grid=[40,40]);This gives a basic plot, but there is a lot we can do with it to make it a bit more interesting. Clicking on the plot with the right mouse button gives a range of options for customising the way the plot looks. Try changing the Style to Patch and the Axes to Boxed, then select Redraw. To view the plot from another angle, click a corner of the plot with the left mouse button and drag the plot around. Then select Redraw to view the plot from the new angle.
The grid lines on the plot show the conditional distributions, which appear to be normal (because they are!). We can plot a conditional for X given Y=1 as follows:
plot(f(x,1,0.8),x=-4..4);Of course, this conditional is not normalised (it won't integrate to one). We can find out how to normalise it with:
int(f(x,1,0.8),x=-infinity..infinity);Now go back to the 3d plot and change the style to Patch and Contour. This should now show horizontal slices through the surface. Drag the plot around so that you are viewing it from above. This should give a contour plot of the surface. You can generate a contour plot directly by doing:
with(plots); contourplot(f(x,y,0.8),x=-4..4,y=-4..4,grid=[40,40]);If you want, you can try going back and redoing your plots for different values of the correlation parameter, rho.
Alternatively, we can produce an animation to show how the density varies as the value of the correlation parameter varies. This requires the plots package to be loaded, but you should have done this earlier.
animate3d(f(x,y,rho),x=-4..4,y=-4..4,rho=-0.8..0.80,frames=20,grid=[30,30]);Change the style to Patch, add the axes and drag the plot round to a convenient viewing angle before selecting Play from the Animation menu. Maple will pause periodically while it generates and stores the frames of the animation to be played. You will probably want to slow down the animation, and view it on continuous play.
To get the marginal distribution for X we can get Maple to do the integration for us.
xmargin:=int(f(x,y,0.8),y=-infinity..infinity);It can be plotted with
plot(xmargin,x=-4..4);You can verify that it integrates to one (and hence that the joint density double-integrates to one) with
volume:=int(xmargin,x=-infinity..infinity);
g:=(x,y,sx,sy,rho)->f(x/sx,y/sy,rho)/(sx*sy);As an exercise, you should plot this for some example parameter values, calculate and plot the marginals, and make sure that the new joint density does intergrate to one. Once we have this density, it is trivial to use it in order to define the density with arbitrary mean with
h:=(x,y,mx,my,sx,sy,rho)->g(x-mx,y-my,sx,sy,rho);Again, you should plot this with some example parameter values, plot some conditionals, look at the marginals, and make sure that they all integrate to one.
|
|
|
| darren.wilkinson@ncl.ac.uk | ||
| http://www.staff.ncl.ac.uk/d.j.wilkinson/ |