Class 36
Between Class Sessions - Prep for Day 36
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) Reading - Rectangles, Targets, and Sums - Sections 4.1 - 4.3
For Unit 3, the reading as well as many practice exercises appear in the document Rectangles, Targets, and Sums
Read sections 4.1-4.3.
Complete a few exercises, and come ready to share what you did with your group. Remember it’s OK if you don’t have the correct answers. The point is to have made a good faith attempt, and then come ready to share what you did.
Using Mathematica, compute each of the following.
- The area under \(g(x) = 2x\) for \(0\leq x\leq 4\), so compute \(\int_0^4 2x dx\).
- The area under \(g(x) = \sqrt{9-x^2}\) for \(-3\leq x\leq 3\), so compute \(\int_{-3}^{3} \sqrt{9-x^2} dx\).
- The centroid of the region under \(g(x) = \sqrt{9-x^2}\) for \(-3\leq x\leq 3\), so compute \[\dfrac{\int_{-3}^{3} x\sqrt{9-x^2} dx}{\int_{-3}^{3} \sqrt{9-x^2} dx}\]
Solutions
Come with your attempts, ready to share with each other. The code below will get you the first bullet point.
Integrate[2 x, {x, 0, 4}](2) Simulation and Probability
- Use the following code to plot the function \(f_1(x) = \frac{1}{\sqrt{0.02\pi}}e^{-\frac{1}{2}\left(\frac{x-4}{0.1}\right)^2}\). That is the Normal probability model (or distribution) with the parameters \(\mu = 4\) and \(\sigma = 0.1\).
f1 <- function(x,mu=0,sigma=1){(1/sqrt(2*pi*sigma^2))*exp(-0.5*((x-mu)/sigma)^2)}
x <- seq(3,5,0.01)
y <- f1(x,4,0.1)
plot(x,y,type='l')- Use the following code to simulate a sample of 25000 random measurements from the Normal probability model with the parameters \(\mu = 4\) and \(\sigma = 0.1\) and plot the density histogram of this sample. This code uses the set.seed() command in R to set the seed so our simulated sample is reproducible (we will all have the exact same random sample of 25000 measurements. We will use the seed 123.
set.seed(123)
tmp <- rnorm(25000, mean=4, sd=0.1)
hist(tmp,probability = TRUE)- Compare your plot of f1 and the density histogram of our sample. Use the following code to plot the sample and f1 together.
hist(tmp,probability = TRUE)
lines(x,y,col=5)- Use the following code to calculate the probability that a measurement in our sample will be less than 3.8.
#Count the number of measurements in our sample less than 3.8.
x <- length(which(tmp < 3.8))
x
#Calculate the probability as the ratio of number of measurements with the characteristic (being less than 3.8) compared to the total number of measurements in our sample.
p <- x/length(tmp)
pAdjust the code above to calculate the probability that a measurement in our sample will be more than 4.7.
Adjust the code again to calculate another probability of your choice using our sample.
Regular Reminders
Skill Practice (KA Homework)
- Begin 3 – Riemann Sums & Definite Integrals
Applied Practice (Project Work)
- Continue working on Project 3 Task 1.
- You should be getting close to having the task complete.
During Class
Code chunks from Rectangles, Targets, and Sums
#Shades a target diagram for a probability mass function.
#Inputs:
# x - a vector of data points
# p - a corresponding vector of probabilities or frequencies
#All widths are 1 unit wide.
draw_pmf <- function(x,p){
xs <- c(rbind(x-1/2,x-1/2,x+1/2,x+1/2))
px <- c(rbind(0,p,p,0))
par(mar=c(2.5,2.5,0.25,0.25))
plot.new()
plot(xs,px,type="l")
polygon(xs,px,col="gray")
}
#Shades a target diagram (shades area under) for a function f from a to b.
#Inputs:
# f - a function f(x)
# a - left end of the target
# b - right end of the target
# num_points - how many point are sent into f for plotting.
draw_target <- function(f,a,b,num_points=100){
x <- c(a,seq(a,b,(b-a)/num_points),b,a)
y <- c(0,f(seq(a,b,(b-a)/num_points)),0,0)
par(mar=c(2.5,2.5,0.25,0.25))
plot(x,y,type = "l")
polygon(x,y,col="gray")
}
#Draws rectangles over the top of a given function.
#The midpoint of top of each rectangle passes through the function.
# f - a function f(x)
# a - left end of graph
# b - right end of graph
# num_rectangles - how many rectangles to plot.
# method - One of "left", "right", or "mid". Defaults to mid.
draw_rect_approx <- function(f,a,b,num_rectangles, method = "mid"){
n <- num_rectangles
dx <- (b-a)/n
x <- c(a,seq(a,b,dx/100),b,a)
y <- c(0,f(seq(a,b,dx/100)),0,0)
par(mar=c(2.5,2.5,0.25,0.25))
plot(x,y,type = "l")
if(method == "left"){
xi <- seq(a+0*dx/2,b-dx/2,dx)
lines(xi,f(xi),type = "h")
lines(xi,f(xi),type = "s")
lines(c(xi[n],xi[n]+dx),f(c(xi[n],xi[n])),type = "l")
lines(c(xi[n],xi[n]+dx),f(c(xi[n],xi[n])),type = "h")
}
else if(method == "right"){
xi <- seq(a+dx,b+dx/2,dx)
lines(xi-dx,f(xi),type = "h")
lines(xi-dx,f(xi),type = "s")
lines(c(xi[n]-dx,xi[n]),f(c(xi[n],xi[n])),type = "l")
lines(c(xi[n]-dx,xi[n]),f(c(xi[n],xi[n])),type = "h")
}
else{#Use midpoint
xi <- seq(a+dx/2,b,dx)
lines(xi-dx/2,f(xi),type = "h")
lines(xi-dx/2,f(xi),type = "s")
lines(c(xi[n]-dx/2,xi[n]+dx/2),f(c(xi[n],xi[n])),type = "l")
lines(c(xi[n]-dx/2,xi[n]+dx/2),f(c(xi[n],xi[n])),type = "h")
}
}Brain Gains
- Use the two properties of probability density functions to explain why \(f\) is not a PDF. \(f(x) = \begin{cases} 2 & \quad 0 \leq x \leq 6 \\ \\ 0 & \text{otherwise}. \end{cases}\)
f <- function(x){2 + 0*x}
draw_target(f,0,6)Find the value \(k\) so that this function is a PDF. \(f(x) = \begin{cases} 2k & \quad 0 \leq x \leq 6 \\ \\ 0 & \text{otherwise}. \end{cases}\)
Use the following code to plot the function \(f_0(x) = \frac{1}{3}\) for \(-1 < x < 2\) and 0 otherwise. This is the uniform probability model (or distribution) with parameters \(a = -1\) and \(b = 2\).
f0 <- function(x,a=0,b=1){(1/(b-a))+0*x}
x <- seq(-1,2,0.1)
y <- f0(x,-1,2)
plot(x,y,type='l',xlim=c(-2,3),ylim=c(0,0.5))Then simulate a sample of 25000 random measurements from this distribution (use the seed 123) and calculate the probability that a measurement in our sample will be less than 1.2.
set.seed(123)
tmp <- runif(25000,a,b)Answers
The function \(f\) is nonnegative (or \(f(x) \geq 0\)) but the area under the function \(f\) is \(A = 2(6) = 12\) which is not 1. So \(f\) is not a probability density function.
The function \(f\) is nonnegative (or \(f(x) \geq 0\)) as long as \(k > 0\) (and we get to pick \(k\)). If \(k = \frac{1}{12}\) then the area under the function \(f\) is \(A = (\frac{2}{12})(6) = 1\). Since \(f\) satisfies the two properties of a pdf it is a pdf.
The total number of measurements in our sample less than or equal to 1.2 is 18327. Comparing this to the number of measurements in our sample we find the probability is 18327/25000 = 0.73308.
set.seed(123)
tmp <- runif(25000,-1,2)
x <- length(which(tmp <= 1.2))
p <- x/length(tmp)
p- Note if we let \(X\) be the random variable which records the \(x\)-coordinate of dropping a dart on the target defined by \(f_0\) we could calculate this same probability \(P(X \leq 1.2)\) from the model. \[ P(X \leq 1.2) = \frac{(1/3)*(1.2+1)}{1} = \frac{11}{15} \approx 0.733333 \]
Key Ideas (Targets & Continuous Random Variables)
The cumulative distribution function (CDF) of a random variable
The cumulative distribution function of a random variable \(X\) is the function \(F(x) = P(X \leq x)\).
The probability density function (PDF) of a continuous random variable
The probability density function, \(f(x)\) of a continuous random variable \(X\) with cumulative distribution function \(F(x)\), is the derivative of \(F(x)\).
- Every probability density function is non-negative, in other words \(f(x) \geq 0\).
- The total area under a probability density function always equals 1.
Any function that is non-negative with a total area of 1 can be interpreted as the probability density function of some random variable.
If a function \(g(x)\) is non-negative and has a finite total area, we can normalize the function (dividing by the area), like we have been doing with our target functions, to make a PDF.
To normalize means to multiply by a factor to make some quantity a desired value. Examples from this semester:
- Compare the intensity of a light bulb to the original intensity. Thus the normalized intensity will be 1 (as a proportion) or 100 (as a percent) when t=0 hours. The “normalization” factor here is \(\frac{1}{\text{original intensity}}\).
- Scale a target so that it has area 1. The “normalization” factor here is \(\frac{1}{\text{total area of original target}}\).
Riemann Sums & Definite Integrals
Riemann sum for \(f\)
Let \(f(x)\) be defined on \(a \leq x \leq b\). Let \(n\) be a positive integer, and divide the interval \([ a, b ]\) into \(n\) subintervals of equal width. From each subinterval choose a point \(x_i\). We call
\(\sum_{i=1}^{n} f(x_i)\Delta x\)
a Riemann sum for \(f\).
Example 1
Let \(g(x) = 4 - x^2\) for \(-1 \leq x \leq 2\) and 0 otherwise.
- Using a Riemann Sum with \(n=10\) and \(x_i\) as the midpoints, approximate the area between \(g\) and the \(x\)-axis.
g <- function(x){4-x^2}
a <- -1
b <- 2
n <- 10
dx <- (b-a)/n
draw_rect_approx(g,a,b,n)
#Start at half of dx to the right of a, and then step by dx.
xi <- seq(a+dx/2,b,dx)
points(xi,g(xi),pch=16,col=2)
segments(xi,rep(0,length(xi)),xi,g(xi),col=2)
Ai <- g(xi)*dx
sum(Ai)- Using a Riemann Sum with \(n=10\) and \(x_i\) as the left end points, approximate the area between \(g\) and the \(x\)-axis.
draw_rect_approx(g,a,b,n,method='left')
#Start a, and then step by dx.
xi <- seq(a,b-dx,dx)
points(xi,g(xi),pch=16,col=2)
segments(xi,rep(0,length(xi)),xi,g(xi),col=2)
Ai <- g(xi)*dx
sum(Ai)- Using a Riemann Sum with \(n=10\) and \(x_i\) as the right end points, approximate the area between \(g\) and the \(x\)-axis.
draw_rect_approx(g,a,b,n,method='right')
#Start a plus dx, and then step by dx.
xi <- seq(a+dx,b,dx)
points(xi,g(xi),pch=16,col=2)
segments(xi,rep(0,length(xi)),xi,g(xi),col=2)
Ai <- g(xi)*dx
sum(Ai)Notice each of these approximations are different. What happens when we increase \(n\)?
Complete the table with \(n = 10, 75, 100, 1000, 50000, 750000,\) and \(4000000\).
g <- function(x){4-x^2}
a <- -1
b <- 2
n <- 10
dx <- (b-a)/n
#Mid: Start at half of dx to the right of a, and then step by dx.
xi.m <- seq(a+dx/2,b,dx)
#Left: Start a, and then step by dx.
xi.L <- seq(a,b-dx,dx)
#Start a plus dx, and then step by dx.
xi.R <- seq(a+dx,b,dx)
Ai.m <- g(xi.m)*dx
Ai.L <- g(xi.L)*dx
Ai.R <- g(xi.R)*dx
sum(Ai.m)
sum(Ai.L)
sum(Ai.R)Using the information from our table, what do you think is the exact area?
It appears that \(\lim_{n \rightarrow \infty} \sum_{i=1}^n g(x_i)\Delta x = 9\).
Definite Integral
For a function \(f(x)\) defined on \(a \leq x \leq b\), the definite integral of \(f\) from \(a\) to \(b\) is
\(\int_a^b f(x) dx = \lim_{n\to \infty }\sum_{i=1}^nf(x_i)\Delta x\),
provided the limit exists. If the limit exists, we say that \(f\) is integrable on \([ a, b ]\).
- \(f(x)\) is the integrand
- \(x\) is the variable of integration
- \(a\) is the lower bound (or lower limit of integration)
- \(b\) is the upper bound (or upper limit of integration)
- Compute \(\int_{-1}^{2} 4 - x^2 dx\) using Mathematica. (Remember to use shift + enter to run the command in Mathematica.)
Integrate[4-x^2,{x,-1,2}]- Is \(g\) the PDF of some random variable? If so, explain. If not, then find a value \(k\) so that \(f(x) = k g(x)\) is the PDF of some random variable.
Solution
The function \(g\) is nonnegative, but the area under \(g\) and above the \(x\)-axis is 9. As such, the function \(g\) is NOT the PDF of some random variable. Letting \(k=\frac{1}{9}\) gives the function \(f(x) = \frac{1}{9}(4 - x^2)\) for \(-1 \leq x \leq 2\) and 0 otherwise. The function \(f\) is the PDF of some random variable.
Group Meeting
Practice with Riemann Sums
- Use a Riemann sum with \(n=25\) to approximate \(E[X]\) for the target function \(f(x) = \frac{1}{2}(x-1)\) for \(1\leq x\leq 3\).
Solution
We can use the formula \(\frac{\sum{x_iA_i}}{\sum{A_i}} = \frac{\sum{x_if(x_i)dx}}{\sum{f(x_i)dx}}\) from Section 1.2 to compute the expected value, which is done below.
f <- function(x){(1/2)*(x-1)}
a <- 1
b <- 3
n <- 25
dx <- (b-a)/n
draw_rect_approx(f,a,b,n)
xi <- seq(a+dx/2,b,dx)
Ai <- f(xi)*dx
sum(xi*Ai)/sum(Ai)The limit of this Riemann Sum as \(n\) goes to infinite is a definite integral. Write down the definite integral that is equal to \(\ds \lim_{n \rightarrow \infty} \sum_{i=1}^{n} \frac{1}{2}(x-1)x \Delta x\).
Approximate the definite integral above by increasing the value of \(n\) from 25 to 50, then 100, etc., until you have a good estimate for limit.
Compute exactly the limit of this Riemann sum (the definite integral) using Mathematica. (Remember to use shift + enter to run the command in Mathematica.)
Integrate[(1/2)*(x-1)*x,{x,1,3}]Practice with identifying PDFs
Work on the chalkboard as you complete these problems, and pass the chalk as you finish each one. Leave up your work so that you can compare your answers with your neighbors. If you notice an answer differs from a neighboring group’s, then have a discussion with them.
- Use the two properties of probability density functions to explain why \(f\) is not a PDF. \(f(x) = \begin{cases} x & \quad 0 \leq x \leq 3 \\ \\ 0 & \text{otherwise}. \end{cases}\)
f <- function(x){x}
draw_target(f,0,3)Find the value \(k\) so that this function is a PDF. \(f(x) = \begin{cases} kx & \quad 0 \leq x \leq 3 \\ \\ 0 & \text{otherwise}. \end{cases}\)
Use the two properties of probability density functions to explain why \(f\) is not a PDF. \(f(x) = \begin{cases} -\frac{1}{2} & \quad 0 \leq x \leq 2 \\ \\ 0 & \text{otherwise}. \end{cases}\)
f <- function(x){-0.5+0*x}
draw_target(f,0,2)- Use the two properties of probability density functions to explain why \(f\) is not a PDF. \(f(x) = \begin{cases} \frac{1}{4}(x-2) & \quad 1 \leq x \leq 5 \\ \\ 0 & \text{otherwise}. \end{cases}\)
f <- function(x){0.25*(x-2)}
draw_target(f,1,5)- Use the two properties of probability density functions to determine whether or not \(f\) is a PDF. \(f(x) = \begin{cases} \frac{1}{8}(x+1) & \quad -1 \leq x \leq 3 \\ \\ 0 & \text{otherwise}. \end{cases}\)
f <- function(x){(1/8)*(x+1)}
draw_target(f,-1,3)- Use the two properties of probability density functions to determine whether or not \(f\) is a PDF. \(f(x) = \begin{cases} 5-x & \quad 2 \leq x \leq 5 \\ \\ 0 & \text{otherwise}. \end{cases}\)
f <- function(x){5-x}
draw_target(f,2,5)Consider the target with top defined by \(f\) (given below). When a dart falls on the point \(( x , y )\), we’ll record just the \(x\)-coordinate and let \(X\) represent this random variable. \(f(x) = \begin{cases} kx - k & \quad 1 \leq x \leq 3 \\ \\ 0 & \quad \text{otherwise}. \end{cases}\)
- Select \(k\) so \(f\) is the PDF of \(X\).
f <- function(x){x-1}
draw_target(f,1,3)Source: Class.36 on byuimath.com