# Eleisha's Segment 22: Uncertainty of Derived Parameters

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).

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
b5 = bfit
#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