Class 14

Between Class Sessions - Prep for Day 14

Please spend around 2 hours working between class sessions, focusing on the tasks below. Use any extra time to complete KnewtonAlta assignments and/or work on Project tasks.

Pick something from your prep today that you can share with your group in class. It might be something new that you learned. It might be questions you have that are still unanswered, or a question along with what helped you eventually answer it. It might be something tricky that you solved. It might be a review topic that helped you remember something. It might be a conversation you had with AI that was helpful. You will have a chance to share this with your peers during class. Come ready to articulate your thinking and questions.

Preparation

(1) Probability Exercise (Tropical Storms Example)

Read the definitions and complete the exercises below.

Definition: Factorial

The product of positive integers from 1 to \(n\) is called \(n\) factorial and is written \(n\)!. We define 0! = 1.

Examples: \(3! = (3)(2)(1) = \prod_{i=1}^{3} i = 6\)

Either of the following commands will calculate \(3!\) in R.

prod(1:3)
factorial(3)

\(n! = (n)(n-1)(n-2) ... (2)(1) = \prod_{k=1}^{n} k\)

Definition: Poisson Random Variable

A Poisson random variable (or Poisson distribution) is a model for the number of events that occur during a fixed time interval. See Poisson distribution on Wikipedia.

The probability mass function (pmf) for the Poisson distribution is \(p(x; \lambda) = \frac{\lambda^x}{x!}e^{-\lambda}\) when \(x = 0, 1, 2, 3, ...\) and zero otherwise where \(\lambda > 0\).

Examples:

Let \(X\) be the random variable that counts the number of hurricanes in Florida in a given year.
Let \(Y\) be the random variable that counts the number of hits on BYU-Idaho’s website in a given minute.
Let \(W\) be the random variable that counts the number of cars that go through an intersection in a given hour. The random variables \(X\), \(Y\), and \(W\) are Poisson random variables.

The number of Florida hurricanes (2000–present) can be thought of as a Poisson random variable, which means we can model the probabilty of having \(x\) tropical storms in a year using \[p(x; \lambda) = \frac{\lambda^x}{x!}e^{-\lambda}\text{ with }x = 0, 1, 2, 3, ...\text{ where }\lambda > 0\] The parameter \(\lambda\) is a value we would like to determine. The code chunk below defines this probability mass function in R, using a default value of \(\lambda = 2\). Start by exploring the function \(p(x; \lambda)\) and how the parameter \(\lambda\) affects the plot.

p <- function(x,lambda=2){
# x must be a whole number
  (lambda^x/factorial(x))*exp(-lambda)
}

x <- seq(0,10,1) # x must be a whole number
plot(x,p(x))
plot(x,p(x,lambda=5))
plot(x,p(x,lambda=2.4)) # lambda just needs to be positive.

Assuming the parameter \(\lambda = 2\), calculate the following and check your answers.

The probability of 4 Florida tropical storms this year.
The probability of 5 Florida tropical storms this year.
The probability of 2 Florida tropical storms this year.
The probability of 3 or more Florida tropical storms this year.
The probability of more than 10 Florida tropical storms this year.
The probability of 2 or less Florida tropical storms this year.

Answers

Let \(X\) be the random variable that counts the number of hurricanes in Florida this year.

Use the following code to graph the pmf of \(X\).

x <- 0:30

par(mfrow=c(1,1),mar=c(2.5,2.5,0.25,0.25))
barplot(p(x),ylim=c(0,0.3),width=rep(1,length(x)),space=1)

\(P(X = 4) = p(4;2) = 0.09022352\)
\(P(X = 5) = p(5;2) = 0.03608941\)
\(P(X = 2) = p(2;2) = 0.2706706\)
\(P(X \geq 3) = \sum_{i=3}^{\infty} P(X=i) = 1 - \sum_{i=0}^{2} P(X=i) = 1 - (P(X=0) + P(X=1) + P(X=2)) = 0.3233236\)
\(P(X > 10) = \sum_{i=11}^{\infty} P(X=i) = 1 - \sum_{i=0}^{10} P(X=i) = 0.000008308\)
\(P(X \leq 2) = P(X=0) + P(X=1) + P(X=2) = 1 - P(X \geq 3) = 0.6766764\)

p(4)
p(5)
p(2)
1 - (p(0) + p(1) + p(2))
1 - sum(p(0:10))
sum(p(0:2))
    1-(1 - (p(0) + p(1) + p(2)))

Now answer the following and come with answers to compare with your group in class.

The probability of less than 8 Florida tropical storms this year using \(\lambda = 3\).
The probability of less than 8 Florida tropical storms this year using \(\lambda = 6\).
The probability of at least 8 Florida tropical storms this year using \(\lambda = 6\).
The probability of more than 12 Florida tropical storms this year using \(\lambda = 5\).

Regular Reminders

Skill Practice (KA Homework)

Complete 3 – Probability Basics assignment in Knewton Alta

Applied Practice (Project Work)

Complete and Submit Project 1 Task 3

During Class

Brain Gains

Let \(h(t) = \sqrt{3t+1}\).

Compute \(h(2)\)
Solve \(h(t) = 2\) for \(t\)

Let \(X_1\), \(X_2\), and \(X_3\) be independent random variables which each return the value of a single toss of one fair six-sided die.

Compute \(P(X_1 = 1, X_2 = 4, \text{ and } X_3 = 5)\)
Compute \(P(X_1 \leq 3, X_2\leq 3,\text{ and }X_3\leq 3)\)

Recall that the probability mass function (pmf) for a Poisson distribution is \[p(x; \lambda) = \frac{\lambda^x}{x!}e^{-\lambda}\text{ with }x = 0, 1, 2, 3, ...\text{ where }\lambda > 0\] Suppose that \(Y\) is a random variable with a Poisson distribution with \(\lambda = 1\).

Compute \(P(Y=1)\)
Compute \(P(Y=0)\)
Compute \(P(Y>1)\)

You may use the code below to help you code in the Poisson distribution.

p <- function(x,lambda=2){
# x must be a whole number
  (lambda^x/factorial(x))*exp(-lambda)
}

Group Meeting

Start by giving each person a moment to share what they chose to prepare for class. Help each other address any questions. When each person has had a chance to share, move on the other activities.

Probability Exercise (Tropical Storms Example from Prep)

Between class we worked with a Poisson distribution to model the number of tropical storms in Florida in a given year. The probability mass function for a Poisson distribution is given by \(p(x; \lambda) = \frac{\lambda^x}{x!}e^{-\lambda}\) with \(x = 0, 1, 2, 3, ...\) where \(\lambda > 0\). You can use the following code to use this function in R.

p <- function(x,lambda=2){
# x must be a whole number
  (lambda^x/factorial(x))*exp(-lambda)
}

Compare your answers that you got for the last 4 questions of yesterday’s prep work. These 4 problems asked you to compute each of the things below.

The probability of less than 8 Florida tropical storms this year using \(\lambda = 3\).
The probability of less than 8 Florida tropical storms this year using \(\lambda = 6\).
The probability of at least 8 Florida tropical storms this year using \(\lambda = 6\).
The probability of more than 12 Florida tropical storms this year using \(\lambda = 5\).

Answers

\[ \begin{aligned} P(X < 8 \mid \lambda = 3) &= \sum_{i=0}^{7} p(i;3) &&\approx 0.9881 \\ P(X < 8 \mid \lambda = 6) &= \sum_{i=0}^{7} p(i;6) &&\approx 0.7440 \\ P(X \geq 8 \mid \lambda = 6) &= 1 - P(X < 8 \mid \lambda = 6) &&\approx 0.2560 \\ P(X > 12 \mid \lambda = 5) &= 1 - \sum_{i=0}^{12} p(i;5) &&\approx 0.0020 \end{aligned} \]

sum(p(0:7, 3))
sum(p(0:7, 6))
1 - sum(p(0:7, 6))
1 - sum(p(0:12, 5))

Activity - Florida Tropical Storms Revisited

Read through the following work, execute the code chunks, and then complete the additional calculations and discuss the questions shown at the end.

Let \(X\) be the random variable that counts the number of tropical storms in Florida this year. Assuming the parameter \(\lambda = 2\), we calculate the following:

#Define the Poisson distribution with a default value for lambda
p <- function(x,lambda=2){
# x must be a whole number
  (lambda^x/factorial(x))*exp(-lambda)
}

#The probability of 4 Florida tropical storms this year (lambda = 2 is assumed)
p(4) 

#The probability of 4 Florida tropical storms this year (lambda = 2 is assigned - should match above. )
p(4,2) 

#The probability of 5 Florida tropical storms this year (lambda = 2 is assumed)
p(5)

#The probability of 2 Florida tropical storms this year (lambda = 2 is assumed)
p(2)

#The probability of $x$ Florida tropical storms this year (lambda = 2 is assumed) for each x from 0 to 10
p(0:10)

The value for \(\lambda\) does not have to be 2. We can instead change it to 5, or some other value, and repeat the computations above.

#The probability of 4 Florida tropical storms this year, using lambda = 5.
p(4,5) 

#The probability of $x$ Florida tropical storms this year, using lambda = 5, for each x from 0 to 10
p(0:10,5)

#The probability of $x$ Florida tropical storms this year, using lambda = 1, for each x from 0 to 10
p(0:10,1)

Additional Calculations and Questions

Compute the probability of 7 Florida tropical storms in a year assuming \(\lambda = 8\)
Compute the probability of 7 Florida tropical storms in a year assuming \(\lambda = 4\)
If we actually observed of 7 Florida tropical storms in a year, would you be more likely to say that \(\lambda =8\) or \(\lambda = 4\)? Why?
Compute the probability of 7 Florida tropical storms in a year assuming \(\lambda = 6.3\). (Wait, can we use decimals for \(\lambda\)?)
Compute the probability of 7 Florida tropical storms in a year using various values for \(\lambda\). Then as a group decide what value you think is the best value to assume for \(\lambda\) if we actually did see 7 tropical storms in a year.
Construct a plot that has \(\lambda\) on the horizontal(\(x\)) axis, and on the vertical(\(y\)) axis we place the probability of 7 Florida tropical storms in a year assuming that value for \(\lambda\). Describe the shape of this plot, and how can we use it to find the “best value” for \(\lambda\).
How does your choice of “best” \(\lambda\) change if we know there will be 4 tropical storms (rather than 7 storms)?
How does your choice for \(\lambda\) change if we know there will be 5 tropical storms (rather than 7 storms)?

Answers

Computed probabilities of 7 storms under several values of \(\lambda\):

\[ \begin{aligned} p(7; \lambda = 8) &\approx 0.1396 \\ p(7; \lambda = 4) &\approx 0.0595 \\ p(7; \lambda = 6.3) &\approx 0.1435 \\ p(7; \lambda = 7) &\approx 0.1490 \end{aligned} \]

Observing 7 storms is more likely under \(\lambda = 8\) than \(\lambda = 4\) (the larger probability), so \(\lambda = 8\) is the more plausible choice. (Discuss with your group why.)
Yes, \(\lambda\) may be any positive real number; only \(x\) must be a non-negative integer.
The probability is maximized at \(\lambda = 7\). (Discuss with your group why \(\lambda = x\) maximizes \(p(x;\lambda)\) here.)
The plot of \(p(7; \lambda)\) vs. \(\lambda\) is a hump that rises, peaks at \(\lambda = 7\), then decays. The “best value” for \(\lambda\) is the location of the peak — i.e., the value that maximizes the probability of the observed outcome. (Discuss with your group.)
If the observation is 4 storms, the peak shifts to \(\lambda = 4\).
If the observation is 5 storms, the peak shifts to \(\lambda = 5\). The pattern:

\[ \lambda_\text{best} = x_\text{observed} \]

p(7, 8)
p(7, 4)
p(7, 6.3)
p(7, 7)
lambda <- seq(0, 20, 0.1)
plot(lambda, p(7, lambda), type = "l")

Activity - Florida Tropical Storms Extended

Let \(X_1\) be the random variable that counts the number of tropical storms in Florida this year, \(X_2\) be the random variable that counts the number of tropical storms in Florida next year, and \(X_3\) be the random variable that counts the number of tropical storms in Florida the year after. Assuming \(X_1\), \(X_2\), and \(X_3\) are independent Poisson random variables each with the the same parameter \(\lambda = 2\). This means we have \[ \begin{aligned} p_{3yr}(x_1, x_2, x_3; \lambda) &= \left(\frac{\lambda^{x_1}}{x_1!}e^{-\lambda}\right) \left(\frac{\lambda^{x_2}}{x_2!}e^{-\lambda}\right) \left(\frac{\lambda^{x_3}}{x_3!}e^{-\lambda}\right) \\ &= \frac{\lambda^{x_1 + x_2 + x_3}}{x_1!x_2!x_3!}e^{-3\lambda} \end{aligned} \] with \(x_1 = 0, 1, 2, 3, ...\), \(x_2 = 0, 1, 2, 3, ...\), and \(x_3 = 0, 1, 2, 3, ...\), (and zero otherwise) where \(\lambda > 0.\) The following code chunk provides two versions to input this function. The first requires you input a vector \(x\) consisting of the number of tropical storms that occurred each year (and can be used if you want to look at a period of 3 years, or even 7 years, by just changing the number of inputs to the vector). The second version will only work with a 3 year period.

p3v1 <- function(x,lambda=2){
# each element of x must be a whole number
  prod((lambda^x/factorial(x))*exp(-lambda))
}

p3v2 <- function(x1,x2,x3,lambda=2){
# x1, x2, and x3 must be whole numbers
  (lambda^(x1+x2+x3)/(factorial(x1)*factorial(x2)*factorial(x3)))*exp(-3*lambda)
}

We can now calculate the following probabilities:

#The probability of 4 Florida tropical storms this year, 4 Florida tropical storms next year, and 8 Florida tropical storms the year after (using $\lambda = 2$ as assumed).
p3v1(c(4,4,8))
#The same probability as above, using the other version
p3v2(4,4,8)
#The same probabilty as above by just multiplying probabilities of independent events together. 
p(4)*p(4)*p(8)

#The probability of 2 Florida tropical storms this year, 5 Florida tropical storms next year, and 3 Florida tropical storms the year after.
p3v1(c(2,5,3))

Again, we can change \(\lambda\) to another value, and see how these probabilities change.

#The probability of 4 Florida tropical storms this year, 4 Florida tropical storms next year, and 8 Florida tropical storms the year after, assuming the parameter $\lambda = 5$.
p3v1(c(4,4,8),5)

#Repeat the above, but with lambda = 1, and then with lambda = 10
p3v1(c(4,4,8),1)
p3v1(c(4,4,8),10)

p3v2(4,4,8,1)
p3v2(4,4,8,10)

If we know there will be 4 Florida tropical storms this year, 4 Florida tropical storms next year, and 8 Florida tropical storms the year after, then our calculations above can help us determine a good value for the parameter \(\lambda\). The value \(\lambda=10\) is a better option that \(\lambda = 1\) because \(p3(4,4,8,10) = 4.029234e-05 > 2.143747e-09 = p3(4,4,8,1)\).

Try various other values for \(\lambda\) and calculate the probability of 4 Florida tropical storms this year, 4 Florida tropical storms next year, and 8 Florida tropical storms the year after.
What criteria would you use to decide “what is the best value for the parameter \(\lambda\)?”
Construct a plot (use p3v2) that has \(\lambda\) on the horizontal(\(x\)) axis, and on the vertical(\(y\)) axis we place the probability of 4, then 4, then 8 Florida tropical storms for the next three years assuming that value for \(\lambda\). Describe the shape of this plot, and how can we use it to find the “best value” for \(\lambda\).
In the plot you made, update your code to use p3v1 instead of p3v2. If you notice an error, that’s OK. The problem is that R does not know how to use two vectors as an input to \(p\) without more information. The sapply function deals with this and shows up in the prep tomorrow.
Use the data from the List of Florida hurricanes (2000–present) for the years 2000-2005 (6 years) to estimate a good value for lambda (you’ll need to use p3v1 now, as there are 6 data points instead of 3).

Answers

Reference probabilities for the observation \((x_1, x_2, x_3) = (4, 4, 8)\):

\[ \begin{aligned} p_{3yr}(4,4,8;\,\lambda = 2) &\approx 6.99\times 10^{-6} \\ p_{3yr}(4,4,8;\,\lambda = 5) &\approx 2.01\times 10^{-3} \\ p_{3yr}(4,4,8;\,\lambda = 10) &\approx 4.03\times 10^{-5} \end{aligned} \]

A good criterion is maximum likelihood: pick the \(\lambda\) that makes the observed outcome most probable. (Discuss with your group why this is reasonable.)
The plot of \(p_{3yr}(4, 4, 8;\,\lambda)\) vs. \(\lambda\) rises, peaks, then decays — same hump shape as the single-year case. The peak occurs at the sample mean:

\[ \begin{aligned} \lambda_\text{best} &= \frac{4 + 4 + 8}{3} \\ &= \frac{16}{3} \\ &\approx 5.33 \end{aligned} \]

with \(p_{3yr}(4, 4, 8;\, \lambda_\text{best}) \approx 2.08\times 10^{-3}\).

Using p3v1 directly with a vector of \(\lambda\) values fails because p3v1 expects a single \(\lambda\), not a vector — sapply will fix this (shows up in tomorrow’s prep).
For 2000–2005 data (using p3v1 with 6 inputs), the best \(\lambda\) is again the sample mean of the 6 observed storm counts. (Look up the actual counts and compute the mean as a group.)

lambda <- seq(0, 20, 0.1)
plot(lambda, sapply(lambda, function(L) p3v2(4, 4, 8, L)), type = "l")
abline(v = (4 + 4 + 8) / 3, lty = 2)

Source: Class.14 on byuimath.com