|
|
@@ -0,0 +1,46 @@
|
|
|
+#!/usr/bin/env python3
|
|
|
+""" LPC stats HW2, Problem 2
|
|
|
+ Author: Caleb Fangmeier
|
|
|
+ Created: Sep. 27, 2017
|
|
|
+"""
|
|
|
+from scipy.stats import beta, poisson
|
|
|
+from scipy.optimize import root
|
|
|
+
|
|
|
+# Declare the experimental results
|
|
|
+N = 25
|
|
|
+M = 353
|
|
|
+κ = 37.6
|
|
|
+x = 1/(1+κ)
|
|
|
+
|
|
|
+p = 0.683
|
|
|
+
|
|
|
+
|
|
|
+def DL(z):
|
|
|
+ """
|
|
|
+ This is just the DL function as defined in the problem statement. It is the
|
|
|
+ cumulative distribution function for p(D|s).
|
|
|
+ """
|
|
|
+ num = 0
|
|
|
+ den = 0
|
|
|
+ for r in range(N+1):
|
|
|
+ b = beta.pdf(x, r+1, M)
|
|
|
+ num += b*poisson.cdf(N-r+1, z)
|
|
|
+ den += b
|
|
|
+ return num/den
|
|
|
+
|
|
|
+
|
|
|
+def DR(z):
|
|
|
+ """
|
|
|
+ Same as DL, but for the opposite side.
|
|
|
+ """
|
|
|
+ return 1 - DL(z)
|
|
|
+
|
|
|
+
|
|
|
+# Scipy's root function finds a zero crossing for the given function.
|
|
|
+# The second argument is just the starting point for the algorithm.
|
|
|
+root_DL = root(lambda z: DL(z) - (1-p)/2, 10).x[0]
|
|
|
+root_DR = root(lambda z: DR(z) - (1-p)/2, 10).x[0]
|
|
|
+print(f"The interval is: [{root_DR:.2f}, {root_DL:.2f}]")
|
|
|
+
|
|
|
+# Output
|
|
|
+# > The interval is: [12.39, 22.78]
|
