HW #3
Dr. Teresa Vasquez, EdD, MBA, MSDS (in progress)
September 09, 2026
Reflection

Based on my handwritten notes below and my tutoring session transcript (session_transcript.Rmd, included in this submission).

How this one went

Estimated time: ~4 hours.

Q1 and Q2 were the light part of this assignment. On Q1 I worked out the complement rule for P(X ≥ 2) cleanly on the first pass — no correction needed. On Q2 I did make a mistake (writing the Gamma mean as 1.5 instead of 6), but I caught it myself immediately by cross-checking a different way: each Exponential(0.5) has mean 2, and I’m summing three of them, so 2+2+2 = 6 had to be right, and 3/0.5 confirmed it.

Q3 was where the real friction was, but it wasn’t a math problem — it was a sourcing problem. The class weights (0.2/0.5/0.3) for the mixture model weren’t sitting in the worksheet; they were supposed to come from lecture. I went looking for the deck, and the course Drive folder stopped at “05 Poisson” — the mixture/continuous lecture just wasn’t there. I ended up photographing an expected-value slide that had the same three numbers (0.2/0.5/0.3) attached to a different problem, which would have been an easy wrong turn to take. It took finding and photographing the actual latent-class slide before I could trust the weights enough to build the hand calc on top of them. That gap is also where the overnight break in this session landed — I stopped mid-search and picked it back up rested the next day, which is exactly the right call when you catch yourself about to build a whole calculation on an unverified number.

Once the weights were locked in, the rest of Q3 went smoothly: I built the general mixture CDF formula, worked P(X ≤ 0) and P(X ≤ 1) by hand, and then caught a small R bug on my own — punif(0, min = 1, max = 1) came back 0, which didn’t match what I expected (0.5), and the fix was realizing the lower bound should have been min = -1, not min = 1.

Spoons it took: on a 12-spoon day, I’d put this one around 5-6 spoons — my own read of the transcript, not a measurement. Q1 and Q2 barely registered; nearly the whole cost was Q3’s slide-hunting detour and the overnight pause in the middle of it.

Growth rating for this session: 7/10. Evidence: I caught my own Q2 mean error through an independent cross-check instead of just accepting the first number I wrote down, and I caught the punif bound bug myself by noticing the output didn’t match the value I already knew it should be. The point I’m not claiming: locating the correct class weights took real searching through external material rather than something I could reason my way to from the homework worksheet alone.

Skills exercised this session: the Poisson pmf and the complement rule, the Binomial-to-Poisson approximation (matching \(n \cdot p = \lambda\)), the relationship between a sum of independent Exponentials and the Gamma distribution, R’s shape/rate parameterization for rexp()/dgamma(), simulating and plotting a kernel density estimate and an empirical CDF, finite mixture models and latent-class simulation, the law of total probability applied to a continuous mixture, and general R fluency (replicate(), sample(), rnorm()/rexp()/runif(), pnorm()/pexp()/punif(), and using mean() on a logical vector to turn a count into a simulated probability).

Things I’m doing well

  • I cross-check a computed value a second way before trusting it. On Q2, catching mean = 1.5 as wrong wasn’t from staring at the formula again — it was from noticing 2+2+2 = 6 had to be the answer since each Exponential(0.5) has mean 2, and that disagreed with what I’d written.
  • I debug R output by comparing it to a value I already trust. The punif bug in Q3 (min = 1 instead of min = -1) was caught because the function returned 0 when I already knew the answer had to be 0.5 — the mismatch is what sent me back to the arguments.
  • I stopped instead of guessing when a needed number wasn’t where it should be. Finding a slide with matching numbers on a different problem didn’t stop me from continuing to search for the actual right one.

Things to work on

  • Verifying that a source value is for this problem, not just that the numbers match. The expected-value slide with 0.2/0.5/0.3 attached to a different problem was a near-miss — I need a faster way to confirm a number belongs to the calculation I’m actually doing before I build on it.
  • R’s shape/rate convention for Gamma and Exponential. I still have to think through, each time, that dgamma() wants rate to match rexp(rate = ...) — this hasn’t become automatic yet.
  • Reading a formula’s inequality precisely before typing it out. The mixture CDF formula needed a fix from \(P(X = x)\) to \(P(X \le x)\) during setup — worth double-checking which one a problem is actually asking for before committing it to the page.

This report was generated using a custom Claude Code skill built by Dr. Teresa Vasquez. The skill takes the assignment PDF and a PDF of my own worked-out solutions — including notes on my reasoning, where I struggled, and what I learned — and digitizes that existing work into this formatted report. AI is not solving the problems; it is transcribing and formatting work I have already completed by hand. The original uploaded files (assignment PDF, handwritten work PDF, and chat transcript) are included in the submission ZIP, which is available to professors.

Question 1

The number of typographical errors on a printed page of a certain textbook follows a Poisson distribution with rate \(\lambda = 1.5\) errors per page.

  1. Simulate 5,000 draws from Poisson(\(\lambda = 1.5\)). Use the simulated draws to estimate the three probabilities below, and report your simulated estimates next to the exact values from the pmf: \(P\)(a page has no errors), \(P\)(a page has at least 2 errors), \(P\)(a page has exactly 3 errors).

  2. Simulate 5,000 draws from Binomial(2000, 0.00075) and, separately, 5,000 draws from Poisson(1.5). Overlay histograms of the two simulated distributions on one figure and comment on how similar they look.

Approach. For the exact side of part (a), I worked from the Poisson pmf, \(P(X=k) = e^{-\lambda}\lambda^k/k!\), plugging in \(\lambda = 1.5\). “No errors” is just \(k=0\), and since \(1.5^0 = 1\) and \(0! = 1\), that collapses down to \(e^{-1.5}\) with nothing left to simplify. “Exactly 3” is \(k=3\), plugged straight into the pmf. “At least 2” needed a different move: the complement of “at least 2” is the finite set \(\{0, 1\}\), so \(P(X \ge 2) = 1 - P(0) - P(1)\) — I don’t have to sum an infinite tail, I just subtract the two cases I can compute directly. For the simulated side, I drew 5,000 values from rpois() and used mean() on a logical vector (draws == 0, draws >= 2, draws == 3) — a mean of TRUE/FALSE values is exactly the fraction of draws meeting the condition, which is the simulated probability.

For part (b), the point being demonstrated is that a Binomial with a large \(n\) and a small \(p\), where \(n \cdot p\) still equals \(\lambda\), behaves like the matching Poisson. Here \(n \cdot p = 2000 \times 0.00075 = 1.5 = \lambda\), so I expected the two histograms to sit almost on top of each other.

\[P(X = k) = \frac{e^{-\lambda}\lambda^k}{k!}, \qquad \lambda = 1.5\] \[P(X \ge 2) = 1 - P(X = 0) - P(X = 1)\]

Exact values from the pmf: \(P(X=0) =\) 0.2231302, \(P(X \ge 2) =\) 0.4421746, \(P(X=3) =\) 0.1255107.

# exact P(X = 0): the Poisson pmf at k = 0. 1.5^0 = 1 and 0! = 1, so this
# collapses straight down to e^(-1.5) with nothing left over.
exact_p0 <- dpois(0, lambda = 1.5)

# exact P(X = 1): the pmf at k = 1, used below as one of the two pieces of
# the complement rule for P(X >= 2).
exact_p1 <- dpois(1, lambda = 1.5)

# exact P(X = 3): the pmf at k = 3, plugged straight in.
exact_p3 <- dpois(3, lambda = 1.5)

# exact P(X >= 2): "at least 2" is the complement of the finite set {0, 1},
# so instead of summing an infinite tail I subtract the two cases I can get
# directly from the pmf. ppois(1, 1.5) already equals P(0) + P(1), so
# 1 - ppois(1, 1.5) is the same thing as 1 - exact_p0 - exact_p1.
exact_p_ge2 <- 1 - ppois(1, lambda = 1.5)

# print all four exact values so they're visible next to the code
exact_p0
## [1] 0.2231302
exact_p1
## [1] 0.3346952
exact_p3
## [1] 0.1255107
exact_p_ge2
## [1] 0.4421746

Now the simulated side, and the Binomial-vs-Poisson overlay for part (b):

set.seed(1) # set seed at 1
m <- 5000   # number of simulations = 5000 or larger

# simulate m draws from Poisson(lambda = 1.5)
draws <- rpois(m, lambda = 1.5)

# simulated estimates: mean() on a logical vector turns a count into a
# fraction of the 5,000 draws, which IS the simulated probability
sim_p0 <- mean(draws == 0)   # simulated P(X = 0)
sim_p_ge2 <- mean(draws >= 2)  # simulated P(X >= 2)
sim_p3 <- mean(draws == 3)   # simulated P(X = 3)

sim_p0
## [1] 0.2302
sim_p_ge2
## [1] 0.4312
sim_p3
## [1] 0.1304
# reset the RNG start for part (b) -- resetting the seed does NOT draw
# anything itself, it just makes the next draws reproducible from a known
# starting point
set.seed(1)

# simulate m draws from Binomial(2000, 0.00075), and separately from
# Poisson(1.5) -- n*p = 2000 * 0.00075 = 1.5 = lambda, so these two should
# look nearly identical
binom_draws <- rbinom(m, size = 2000, prob = 0.00075)
poisson_draws <- rpois(m, lambda = 1.5)

par(bg = "#0a0a0a", col.axis = "#f7f1e6", col.lab = "#f7f1e6",
    col.main = "#f7f1e6", fg = "#f7f1e6")
hist(binom_draws, breaks = seq(-0.5, 12.5, by = 1),
     col = rgb(0, 0, 1, 0.4), border = "black",
     main = "Binomial(2000, 0.00075) vs Poisson(1.5)", xlab = "count")
hist(poisson_draws, breaks = seq(-0.5, 12.5, by = 1),
     col = rgb(1, 0, 0, 0.4), border = "black", add = TRUE)
legend("topright", c("Binomial(2000, 0.00075)", "Poisson(1.5)"),
       fill = c(rgb(0, 0, 1, 0.4), rgb(1, 0, 0, 0.4)), text.col = "#f7f1e6",
       bty = "n")

My simulated estimates (0.2302, 0.4312, 0.1304) land close to the exact pmf values above, well within what I’d expect from sampling noise at 5,000 draws. And for part (b): the two histograms sit almost exactly on top of each other, which is exactly what should happen when a Binomial’s \(n \cdot p\) matches the Poisson’s \(\lambda\) — the large-\(n\), small-\(p\) Binomial is well approximated by the matching Poisson.

↑ Back to top

Question 2

Confirm computationally that a sum of independent Exponential random variables follows a Gamma distribution.

  1. In R, simulate 5,000 values of the sum of 3 independent Exponential(0.5) draws. Overlay a histogram of your simulated sums with the theoretical Gamma(3, 0.5) density, evaluated over an appropriate grid of \(x\) values.

  2. Compare the sample mean and sample variance of your simulated sums to the theoretical Gamma(3, 0.5) mean and variance.

Approach. For part (a), the parameterization is the part to watch: R’s dgamma() takes a rate, not a scale, so Gamma(shape = 3, rate = 0.5) is the matching theoretical distribution for a sum of three rexp(rate = 0.5) draws — the rate has to agree on both sides. For part (b), the theoretical mean and variance of a Gamma with shape \(a\) and rate \(b\) are \(\text{mean} = a/b\) and \(\text{variance} = a/b^2\). My first pass at the mean gave 1.5, which was wrong — I caught it by checking the answer a different way: each individual Exponential(0.5) has mean \(1/0.5 = 2\), and I’m summing three of them, so the total mean has to be \(2+2+2=6\). That matched \(3/0.5 = 6\) once I redid the division correctly, not \(1.5\).

\[\text{mean} = \frac{a}{b} = \frac{3}{0.5}, \qquad \text{variance} = \frac{a}{b^2} = \frac{3}{0.5^2}, \qquad a = 3,\ b = 0.5\]

Theoretical Gamma(3, 0.5) mean = 6, variance = 12.

# Gamma(shape = a, rate = b) parameters matching the sum of 3
# Exponential(0.5) draws -- a = number of Exponentials summed, b = their
# shared rate.
a <- 3
b <- 0.5

# theoretical mean = a / b. My first pass at this gave 1.5 -- wrong. I
# caught it by cross-checking a different way: each Exp(0.5) alone has mean
# 1/0.5 = 2, and summing 3 of them gives 2+2+2 = 6, which is what 3/0.5
# actually equals (not 1.5).
theoretical_mean <- a / b

# theoretical variance = a / b^2 -- this one I got right on the first pass.
theoretical_var <- a / b^2

theoretical_mean
## [1] 6
theoretical_var
## [1] 12

Now the simulation, the histogram-vs-density overlay, and the sample mean/variance comparison:

set.seed(2) # set seed at 2
m <- 5000   # number of simulations = 5000 or larger

# simulate m sums, each the sum of 3 independent Exponential(0.5) draws
sums <- replicate(m, sum(rexp(3, rate = 0.5)))

par(bg = "#0a0a0a", col.axis = "#f7f1e6", col.lab = "#f7f1e6",
    col.main = "#f7f1e6", fg = "#f7f1e6")
hist(sums, breaks = 30, freq = FALSE,
     col = rgb(0, 1, 0.9, 0.5), border = "black",
     main = "Sum of 3 Exp(0.5) vs Gamma(3, 0.5)", xlab = "sum")
# dgamma with rate = 0.5 matches rexp(rate = 0.5) -- same rate on both sides
curve(dgamma(x, shape = 3, rate = 0.5), add = TRUE, col = "#FF00A8", lwd = 3)
legend("topright", c("Simulated sums", "Gamma(3, 0.5) density"),
       fill = c(rgb(0, 1, 0.9, 0.5), NA), border = c("black", NA),
       lty = c(NA, 1), lwd = c(NA, 3), col = c(NA, "#FF00A8"),
       text.col = "#f7f1e6", bty = "n")

# sample mean and variance of the simulated sums, to compare against the
# theoretical values computed above
sample_mean <- mean(sums)
sample_var <- var(sums)
sample_mean
## [1] 6.005578
sample_var
## [1] 12.18892

The simulated sums line up closely with the Gamma(3, 0.5) curve. My sample mean (6.006) and sample variance (12.189) both land close to the theoretical values of 6 and 12 — well within what 5,000 simulated draws should produce, confirming that a sum of 3 independent Exponential(0.5) draws really does follow a Gamma(3, 0.5) distribution.

↑ Back to top

Question 3

This question generates a pdf (or CDF) picture by generating data directly from the finite mixture model below.

  1. In R, simulate 5,000 draws from the finite mixture distribution and plot the pdf and CDF.

  2. Estimate \(P(X \le 0)\) and \(P(X \le 1)\) from your simulation, and compare them to your hand calculations using pnorm(), pexp(), and punif().

The class-conditional distributions and their probabilities, from lecture:

Class Distribution of X | Class Probability
0 N(0, 2) 0.2
1 Exp(3) 0.5
2 Uniform(-1, 1) 0.3

Approach. The weights in that table (0.2/0.5/0.3) weren’t in the worksheet — they came from lecture, and finding the right slide took real work: the course Drive’s deck folder stopped at “05 Poisson,” and I initially photographed an expected-value slide that happened to have the same three numbers attached to a completely different problem. Once I found and photographed the actual latent-class mixture slide, the weights were confirmed (they sum to 1, a valid pdf) and I could trust them.

For simulating \(X\), the mechanism is: first draw a latent class (0, 1, or 2) with those probabilities, then draw \(X\) from that class’s own distribution. One thing to watch: R’s rnorm()/pnorm() take a standard deviation, not a variance, so \(N(0,2)\) meaning variance 2 needs sd = sqrt(2).

For part (b), because \(X\) comes from a mixture, \(P(X \le x)\) is a weighted average across the three classes — I originally wrote the formula with \(P(X = x)\) instead of \(P(X \le x)\) and had to fix that before plugging in numbers:

\[P(X \le x) = w_0 \cdot P(X \le x \mid \text{class } 0) + w_1 \cdot P(X \le x \mid \text{class } 1) + w_2 \cdot P(X \le x \mid \text{class } 2)\]

Class 0’s term uses pnorm, class 1’s uses pexp, class 2’s uses punif. I worked this out once for \(x=0\) and once for \(x=1\):

\(x\) Class 0 term Class 1 term Class 2 term \(P(X \le x)\)
0 0.5 0 0.5 (0.2)(0.5) + (0.5)(0) + (0.3)(0.5) = 0.25
1 0.7602 0.9502 1 (0.2)(0.7602) + (0.5)(0.9502) + (0.3)(1) = 0.927

Hand calculation: \(P(X \le 0) =\) 0.25, \(P(X \le 1) =\) 0.927.

# class weights, confirmed from the lecture's latent-class mixture slide
# (sums to 1, so it's a valid pdf over the 3 classes)
w0 <- 0.2
w1 <- 0.5
w2 <- 0.3

# class-conditional P(X <= x) terms for x = 0. N(0,2) means VARIANCE 2, so
# sd = sqrt(2), not sd = 2. punif's lower bound has to be -1 -- I first
# wrote punif(0, min = 1, max = 1) here by mistake, which returned 0
# instead of the 0.5 I already knew was right, and that mismatch is what
# sent me back to fix min = -1.
term0_x0 <- pnorm(0, mean = 0, sd = sqrt(2))  # class 0: N(0,2)
term1_x0 <- pexp(0, rate = 3)                 # class 1: Exp(3)
term2_x0 <- punif(0, min = -1, max = 1)       # class 2: Uniform(-1,1)

# same three class-conditional terms, now at x = 1
term0_x1 <- pnorm(1, mean = 0, sd = sqrt(2))
term1_x1 <- pexp(1, rate = 3)
term2_x1 <- punif(1, min = -1, max = 1)

# P(X <= x) is the weight-average of the three class terms -- this is the
# law of total probability applied to a continuous mixture
p_le_0 <- w0 * term0_x0 + w1 * term1_x0 + w2 * term2_x0
p_le_1 <- w0 * term0_x1 + w1 * term1_x1 + w2 * term2_x1

term0_x0
## [1] 0.5
term1_x0
## [1] 0
term2_x0
## [1] 0.5
term0_x1
## [1] 0.7602499
term1_x1
## [1] 0.9502129
term2_x1
## [1] 1
p_le_0
## [1] 0.25
p_le_1
## [1] 0.9271565

Now the simulation itself, the pdf/CDF plots, and the simulated estimates to compare against the hand calculation above:

set.seed(3) # set seed at 3
m <- 5000   # number of simulations = 5000 or larger

# draw the latent class (0, 1, or 2) for each of the m draws, according to
# the class probabilities confirmed above
classes <- sample(0:2, size = m, replace = TRUE, prob = c(w0, w1, w2))

# given the latent class drawn above, simulate X from the corresponding
# class-conditional distribution
x <- numeric(m) # empty container for 5000 numbers
x[classes == 0] <- rnorm(sum(classes == 0), mean = 0, sd = sqrt(2)) # N(0,2)
x[classes == 1] <- rexp(sum(classes == 1), rate = 3)                # Exp(3)
x[classes == 2] <- runif(sum(classes == 2), min = -1, max = 1)      # Uniform(-1,1)

par(bg = "#0a0a0a", col.axis = "#f7f1e6", col.lab = "#f7f1e6",
    col.main = "#f7f1e6", fg = "#f7f1e6", mfrow = c(1, 2))

# left panel: kernel density estimate of the simulated draws -- approximates
# the true mixture pdf, 0.2*dnorm(x,0,sqrt(2)) + 0.5*dexp(x,3) + 0.3*dunif(x,-1,1)
plot(density(x), col = "#00FFE1", lwd = 3,
     main = "PDF (density) of simulated X", xlab = "x", ylab = "density")
polygon(density(x), col = rgb(0, 1, 0.9, 0.3), border = NA)

# right panel: empirical CDF of the simulated draws -- approximates
# F(x) = w0*P(N(0,2)<=x) + w1*P(Exp(3)<=x) + w2*P(Unif(-1,1)<=x)
plot(ecdf(x), col = "#FF00A8", lwd = 2, pch = NA,
     main = "CDF (ecdf) of simulated X", xlab = "x", ylab = "F(x)")

# simulated estimates of P(X <= 0) and P(X <= 1), to compare against the
# hand-calculated values above -- mean() on a logical vector again turns a
# count into a probability
sim_p_le_0 <- mean(x <= 0)
sim_p_le_1 <- mean(x <= 1)
sim_p_le_0
## [1] 0.2558
sim_p_le_1
## [1] 0.9258

My simulated estimates — \(P(X \le 0) \approx\) 0.256 and \(P(X \le 1) \approx\) 0.926 — land close to the hand calculations of 0.25 and 0.927, confirming the mixture CDF formula and the class weights I finally tracked down.

↑ Back to top