Means and Moments and Exploring New Plotting Techniques

Introduction

OBJECTIVE: Extend the concept of the moments of a density function of a solid to applications in probability and engineering.

What do means and multiple integrals have to do with probabilities? How can the method of moments be used to help determine whether or not an object will float in an upright position? These questions will be answered as you explore this module. You will also get to practice new and interesting ways to plot functions.

Technology Guidelines

NOTE:  If you have just finished a module, restart Mathematica or Quit Kernel under the Evaluation pull-down menu before executing a new module.
TO OPEN CELLS, click on the arrow to the left of the cell or put your cursor on the right cell bracket and double click. To open all cells, highlight the bracket on the far right and select Grouping under the Cell pull-down menu and select Open All Subgroups.
TO STOP AN EXECUTION
Select the Evaluation pull-down menu and click on Abort Evaluation.
ORDER OF EXECUTION
Execute cells in the order given. Do not skip any Input cells within a given notebook.
SAVING NOTEBOOKS
You can save anytime to any directory you choose, and it is wise to save often. However, before you do your final save to reduce the size of your file, it is a good idea to delete all your output by selecting the Delete All Output selection under the Cell pull-down menu.
EXPERIENCING MAJOR PROBLEMS
Save if appropriate, then shut down Mathematica and start it up again.

Part I: Masses and Moments Applied to Probability Functions

In multivariable probability, probability density functions are defined over specified regions (domains) in the same way that density functions are defined over regions of a solid. These probability functions are designed so that the total "mass" is always 1. Probabilities are then defined to be the integrals of the density functions over a particular subregion of the domain. Consider the following example.

Suppose that a national fast-food outlet is interested in the joint behavior of the random variables x, defined as the total time between a customer's arrival at the store and departing from the service window, and y, the time that a customer waits in line before reaching the service window. Since x contains the time a customer waits in line, x must be greater than y. The relative frequency distribution of observed values of x and y can be modeled by the probability density function:  
f(x,y) = 0.04(e^(-0.1x-0.3y ))  for  0 ≤ y ≤ x < ∞  and 0 elsewhere, where x and y are measured in minutes.

Clear[f,x,y]
f[x_,y_]:=.04(E^(-.1x-.3y ))

We examine the domain of this function.

Plot[x,{x,0,50},AxesLabel->{"time in system","time in line"},Filling->Axis]

We begin by verifying that the integral of the density function over the entire domain is 1.

Integrate[f[x,y],{x,0,∞},{y,0,x}]

To compute the average or expected values of x and y over the domain, we find the first moments of the density function with respect to the x and y axes. The second integral is difficult to evaluate as x approaches ∞, so we substitute a large value of x for the upper bound.

Print[ "Expected time in system = ",Integrate[x f[x,y],{x,0,∞},{y,0,x}]," minutes."]
Print[ "Expected time waiting in line = ",Integrate[y f[x,y]//Simplify,{x,0,1000000},{y,0,x}]," minutes."]

Now we compute the probability that y is at least one-half as large as x. First we can draw the region over which we will integrate (the part for x < 5).

Plot[{x,(x/2)},{x,0,5},AxesLabel->{"time in system","time in line"},Filling->{1->{2}}]

The probability of this event is the integral of the density function over this region.

Print["The probability that the time in line is at least half as long as the time in the system is ",Integrate[f[x,y],{x,0,∞},{y,1/2x,x}]]

If we want to know the probability that a customer spends no more than six minutes at the service window itself, we represent that quantity x - y, and we integrate over the portion of the domain where x - y ≤ 6.

p1=Plot[{x,x-6},{x,6,20},Filling->{1->{{2},GrayLevel[0.8`]}}];
p2=Plot[x,{x,0,6},Filling->Axis];
Show[p1,p2,PlotRange->All,AxesOrigin->{0,0},AxesLabel->{"time in system","time in line"}]

Looking at the domain, you can see that it is easier here to integrate with respect to x first.

Print["The probability that the time spent at the service desk is no more than 6 minutes is ",Integrate[f[x,y],{y,0,∞},{x,y,y+6}]]

You Try It: Part I

Try your hand at computing the probability that the time at the service desk will be at least 8 minutes. Change the entries in red to the appropriate terms.

Print["The probability is ",Integrate[f[x,y],{y,0,∞},{x,y,y+6}]]

Compute the probability that the time in line is no more than 30% of the time in the system. Change the entries in red to the appropriate terms.

Print["The probability is ",Integrate[f[x,y],{x,0,∞},{y,1/2x,x}]]

Part II: Using Moments to Determine Acceptable Heights

Designing a Buoy

Suppose that you have a buoy in the form of a solid hemisphere of radius 40 centimeters and uniform density that is topped with a solid cylinder of the same material and radius 40 centimeters. This buoy floats with the center of the hemisphere above the water surface. Determine the maximum height that the cylinder may attain that permits the buoy to continue floating upright. To picture this, we plot the shape of the body.

Clear[x,y,z,u,v,t,r,θ]
hemisplot=RevolutionPlot3D[-(Sqrt[1600-(r^2)]),{r,0,40},{θ,0,2 π},AxesLabel->{x,y,z}];
cylpts=Table[{40 Cos[t],40 Sin[t],u},{t,0,2 π,(π/12)},{u,0,40}];
cylplot=ListPointPlot3D[cylpts];
topplot=RevolutionPlot3D[40,{r,0,40},{θ,0,2 π},AxesLabel->{x,y,z}];
flagpolepts=Table[{0,0,t},{t,35,90}];
flagpts=Table[{0,u,v},{u,0,15},{v,70,90}];
fp1=ListPointPlot3D[flagpolepts];
fp2=ListPointPlot3D[flagpts];
Show[hemisplot,cylplot,topplot,fp1,fp2,PlotRange->All]

You may assume that the flag on top is weightless. 

The method of moments can help us. It can be determined that the buoyancy force exerted by the body of water due to the amount of water displaced will be in a direction through the point that would be the midpoint of the sphere, as long as the object is upright or tilted so that no part of the cylinder is immersed. Hence, if the center of gravity of the solid is at that point or closer to to the bottom of the hemisphere, the object will stay upright. If the center of gravity is in the region of the cylinder, the object will tilt. 

We first determine the z-coordinate of the center of gravity of the hemisphere. Recall that the volume of a hemisphere with radius of r is 2/3 π (r^3).

volhemis=2/3 π (40^3);
hemis=-(Sqrt[1600-(x^2)-(y^2)]);
zbar=Integrate[z ,{x,-40,40},{y,-(Sqrt[1600-(x^2)]),(Sqrt[1600-(x^2)])},{z,hemis,0}]/volhemis;
Print["z-coordinate for center of gravity of hemisphere = ",zbar]

The z-coordinate of the center of mass of the cylinder is h/2, so if we equate the first moment of the hemisphere to the first moment of the cylinder, we can solve for h.  Since the radius is 40 here, the volume of the cylinder is 1600 π h.

solh=Solve[-zbar*volhemis==h/2*(1600π h),h]
Print["The height of the cylinder can be at most ",solh[[2,1,2]]," or ",N[solh[[2,1,2]]]," centimeters."]

In general, if the radius were R, the maximum height of the cylinder could be at most  ((R(Sqrt[2]))/2).

You Try It: Part II

Experiment with some of the 3-D plot functions presented in this section. First, just practice with a cylindrical plot. It is understood in these plots that z is a function of r and θ . Try some of your own functions. Change the terms in red.

z=r ((Cos[θ])^2);
newcyl=RevolutionPlot3D[z,{r,0,10},{θ,0,2 π},AxesLabel->{x,y,z}]

You can also create spherical plots. In these plots, it is assumed that ρ is a function of ϕ and θ. Try some of your own by changing the terms in red. Not everything will give satisfying results.

ρ=θ Cos[ϕ];
SphericalPlot3D[ρ,{ϕ,0,π},{θ,0,2 π},AxesLabel->{x,y,z}]

You may have noticed that to draw the flag we had to first generate points and then plot the points. The difference between the two was that the pole could be thought of as a line, whereas, the flag was more like a surface. Experiment with the following commands and draw something else. Remember to leave one parameter in the first expression and two parameters in the second expression.

newflagpolepts=Table[{0,0,t},{t,0,50,0.1}];
newflagpts=Table[{0,u,v},{u,0,20,.1},{v,35,50,.1}];
fp1=ListPointPlot3D[newflagpolepts];
fp2=ListPointPlot3D[newflagpts];
Show[fp1,fp2,AspectRatio->1,Boxed->False,Axes->False]