Difference between revisions of "Eleisha's Segment 22: Uncertainty of Derived Parameters"

Jump to navigation Jump to search

To Compute:

1. In lecture slide 3, suppose (for some perverse reason) we were interested in a quantity $\displaystyle f = b_3/b_5$ instead of $\displaystyle f = b_3b_5$ . Calculate a numerical estimate of this new $\displaystyle f$ and its standard error.

The variance of $\displaystyle f$ can be calculated as:

$\displaystyle \text{Variance} = \bigtriangledown f \Sigma \bigtriangledown f^T$

In this case:

$\displaystyle \bigtriangledown f = (0, 0, \frac{1}{b_5}, 0, \frac{-b_3}{b_5^2} )$

So,

$\displaystyle \text{Variance} = \bigtriangledown f \Sigma \bigtriangledown f^T = \frac{\Sigma_{33}}{b_5} - \frac{2\Sigma_{35} b_3}{b_5^3} + \frac{\Sigma_{55}b_3^2}{b_5^4}$

Using the given estimates of $\displaystyle b_3, b_5$ and the covariance matrix $\displaystyle \Sigma$ this can be calculated to be approximately equal to: 0.0145.

Therefore the standard error of $\displaystyle f$ is approximately: 0.1204 (take the square root of the variance).

As for $\displaystyle f$ , $\displaystyle f = \frac{b_3}{b_5} \approx 0.44377$

In summary,

$\displaystyle f = \frac{b_3}{b_5} \approx 0.44377 \pm 0.1204$ .

2. Same set up, but plot a histogram of the distribution of $\displaystyle f$ by sampling from its posterior distribution (using Python, MATLAB, or any other platform).

Below is a histogram of the distribution of $\displaystyle f$ by sampling from its posterior distribution

Histogram:

Below is the python script that was used to generate the histogram and perform all of the calculations for question 1.

import numpy as np
import matplotlib.pyplot as plt

bfit = [1.1235, 1.5210, 0.6582, 3.2654, 1.4832]
b3 = bfit[2]
b5 = bfit[4]
#print b3, b5

covar = np.loadtxt("covar.txt") #Import the covariance data from the slides (I put it in a text file)
#print covar

cov_33 = covar[2,2] #Save the needed covariance elements as variables
cov_35 = covar[2,4]
cov_55 = covar[4,4]
#print cov_33, cov_35, cov_55

func_covar = (cov_33/b5**2) - (2*cov_35*b3)/(b5**3) + (cov_55*b3**2)/(b5**4) #

print "Variance of f: " + str(func_covar)
print "Standard error of f: " + str(np.sqrt(func_covar))
f_estimate = (b3/b5) #Get estimate of new quantity f
print "Estimate of f: " + str(f_estimate)

#These next few lines sample f from the its posterior distribution
bees = np.random.multivariate_normal(bfit, covar, 10000)
#print bees
generated_fs = bees[:,2]/bees[:,4]
#print generated_fs

#Plot the sample in a histogram
plt.hist(generated_fs, 30)
plt.xlabel("b3/b5")
plt.ylabel("Frequency")
plt.savefig("Eleisha_HW22_Figure.png")
plt.show()


Sample output

Variance of f: 0.0145062470173
Standard error of f: 0.120441882322
Estimate of f: 0.443770226537


To Think About:

1. Lecture slide 2 asserts that a function of normally distributed RVs is not, in general, normal. Consider the product of two independent normals. Is it normal? No! But isn't the product of two normal distribution functions (Gaussians) itself Gaussian? So, what is going on?

2. Can you invent a function of a single normal N(0,1) random variable whose distribution has two separate peaks (maxima)? How about three? How about ten?

Back To: Eleisha Jackson