Hints for Lab 9: Bayesian Analysis

For the class on Monday, February 12th

See also

Go back to Lab 9

A. Comparing Samples

Hints for Part A

• The probability density function for the exponential distribution is \(p(x; \lambda) = \lambda e^{-\lambda x}\)

• It might be easier to carry out the calculation by taking a logarithm first.

B. Coin Flips

Hints for Part B

• The likelihood function for coin flips is \(p^K (1-p)^{N-K} \binom{N}{K}\), where \(p\) is the probability the coin shows heads, \(N\) is total number of flips, and \(K\) is the number of heads you get. See Lab 3 Demo for example code.

• Here’s a full example to calculate the confidence interval for an unnormalized pdf, \(p(x) \propto x(1-x)\), defined for \(x \in [0,1]\). Note that this pdf peaks at \(x=0.5\).

  from scipy.integrate import trapezoid, cumulative_trapezoid
  
  def unnormalized_pdf(x):
      return x * (1-x)
  
  x = np.linspace(0, 1, 1001)  # should fill the integration interval
  y = unnormalized_pdf(x)
  total_integral = trapezoid(y, x)
  
  peak = 0.5  # you'll need to find this value for your unnormalized_pdf
  
  # calculate the right bound
  x_from_peak = np.linspace(peak, 1, 1001)
  y_from_peak = unnormalized_pdf(x_from_peak)
  integral_from_peak = cumulative_trapezoid(y_from_peak, x=x_from_peak, initial=0) / total_integral
  print("right interval bound:", x_from_peak[np.searchsorted(integral_from_peak, 0.95/2)])
  
  # calculate the left bound
  x_from_peak = np.linspace(peak, 0, 1001)
  y_from_peak = unnormalized_pdf(x_from_peak)
  integral_from_peak = cumulative_trapezoid(y_from_peak, x=x_from_peak, initial=0) / total_integral
  print(" left interval bound:", x_from_peak[::-1][np.searchsorted(integral_from_peak[::-1], -0.95/2)])