Is Error Propogation using partial derivatives exact or an approximation? If it is an approximation, when and why is it valid?
If you have a constant error in your independent variable, you almost always require weighted least-squares instead of ordinary least-squares for non-linear regression. Why?
Is attenuation possible in non-linear weighted least-squares?
If you have a measurement error of $\sigma_{\lambda} = 2.0$ and you are using the equation:
$$ y = \lambda e^{-\lambda x}$$
what is your uncertainty in $y$? Assume $x$ is exact.
What is the equation for standard error in residuals with constant uncertainty in the independent variable?
What is the equation for standard error in residuals with variable uncertainty in the dependent variable?
Give two examples of when to use weighted least-squares
The Arrhenius Equation allows you to find a reaction rate at a temperature. It is:
$$ k = A e^{-\frac{E_a}{RT}} $$where $k$ is the reaction rate, $A$ is an empirical constant, $E_a$ is the activation energy, $R$ is the Universal gas constant and $T$ is temperature. In this homework problem, you will determine $A$ and $E_a$ and their confidence intervals.
The chemical reaction that you're observing is the unimolecular dissolution of a 125 Dalton compound in water. This reaction occurs over the course of days inside a large tank which has a controlled temperature. You take an approximately 250mL sample every four hours from the tank and can compute a concentration of the compound by drying the water off and weighing the dry mass. The measurement process takes about 10 minutes and so you may estimate your measurement error in time as 5 minutes, since the reaction continues during the drying process and taking the sample requires some time. Your volume measurements have a measurement error of 0.25 mL and your balance has a measurement error of 0.5 g. You know the temperature exactly.
Knowing the information above, you may fit your concentration data vs time to an exponential decay:
$$C(t) = C_0 e^{-kt}$$and thus compute a reaction rate. Using reaction rates at four temperatures, you may estimate the activation energy and pre-exponential empricial factor $A$ from the Arrhenius equation. Carefully propogate error at all steps and ensure your final confidence interval includes error from the measurements and the regression. Report the reaction rate in mol / s$\cdot$L and activation energy in kJ. Also, plot concentration vs time for one experiment with error bars on the data and your best fit line. Also plot the best-fit Arrhenius equation and your data with error bars. As a hint, here's the basic workflow of this problem:
Hints: you'll need to use basinhopping to do the minimization. Also, you should have variable measurement error for concentration due to error propogation.
In [2]:
%matplotlib inline
import numpy as np
import scipy
import matplotlib.pyplot as plt
import numpy as np
volume_25 = np.array([276.1, 240.5, 277.6, 260.3, 296.0, 237.1, 247.5, 263.2, 232.2, 264.1, 236.0, 262.2, 293.4, 268.6, 222.7, 280.5, 239.3, 267.6, 242.6, 225.5, 278.0])
mass_25 = np.array([ 13.8, 10.1, 12.8, 12.3, 12.8, 9.7, 9.2, 10.3, 8.2, 8.1, 7.8, 8.1, 8.4, 7.6, 5.3, 7.4, 5.2, 5.4, 5.2, 4.2, 4.9])
time = np.array([0, 480, 960, 1440, 1920, 2400, 2880, 3360, 3840, 4320, 4800, 5280, 5760, 6240, 6720, 7200, 7680, 8160, 8640, 9120, 9600])
Let's first compute the measurement error in concentration due to our error in volume and mass:
$$c = f(m,v) = \frac{m}{MW v}$$Notice we do not need to consider molecular weight, since it is known exactly.
$$\sigma_{c}^2 = \left(\frac{\partial f(m, v)}{\partial v}\right)^2 \sigma_v^2 + \left(\frac{\partial f(m, v)}{\partial m}\right)^2 \sigma_m^2$$$$\sigma_c^2 = \frac{m^2}{MW^2 v^4}\sigma_v^2 + \frac{1}{MW^2 v^2}\sigma_m^2$$
In [3]:
#Setup constants and divens
MW = 125.0
terr = 5.0
merr = 0.5
verr = 0.25
#compute concentration and error in concentration
c_25 = mass_25 / volume_25 * (1000 / MW)
c_25_err = np.sqrt( 1.0 / (volume_25)**2 * merr**2 + mass_25**2 / volume_25**4 * verr**2 ) * 1000.0 / MW
print c_25_err, c_25
#plot em
plt.errorbar(time, c_25, xerr=terr, yerr=c_25_err, capthick=3)
plt.ylabel('Concentration [mol/L]')
plt.xlabel('Time [minutes]')
plt.show()
For the analysis, we have two adjustable parameters: $C_0$ and $k$. This requires weighted-least squares in the non-linear regression setting. First, let's compute our weights. The equation is:
$$W_i = \frac{1}{\sigma^2_{c_i} + \sigma_{f(t_i)}^2}$$$$\sigma_{f(x_i)}^2 = \left(\frac{\partial C_0e^{-kt}}{\partial t}\right)^2 \sigma_{t}^2 = C_0^2k^2e^{-2kt} \sigma_t^2$$We'll just use $\hat{C}_0$ as our first concentration measurement. However, we cannot forget it is an adjustable parameter in our error analysis.
In [4]:
import scipy.stats
#This our weighted SSR following the derivation above.
def WSSR(k, t, c, terr, cerr):
w = 1.0 / (cerr**2 + terr**2 * (k * np.exp(-k * t))**2)
chat = c[0] * np.exp(-k * t)
return np.sum(w*(chat - c)**2)
#WE need to bound this and use basin_hopping because it is a non-linear problem
mini_args = {'bounds': [(10**-10, 10)]}
result = scipy.optimize.basinhopping(lambda s: WSSR(s, time, c_25, terr, c_25_err), x0=0.1, minimizer_kwargs=mini_args, niter=1000)
print result
#now get the predicted chats and compare
chat = c_25[0] * np.exp(-result.x * time)
plt.errorbar(time, c_25, xerr=terr, yerr=c_25_err, capthick=3)
plt.plot(time, chat)
plt.show()
We need to find the standard error in the fit rate. We start by building our $\mathbf{F}$ matrix:
$$F_{ij} = \left. \frac{\partial f(\beta, x)}{\partial \beta_j}\right|_{x=x_i}$$$$F_{i1} = \left. \frac{\partial (C_0e^{-kt})}{\partial k}\right|_{t=t_i} = -t_iC_0e^{-kt_i}$$$$F_{i2} = \left. \frac{\partial (C_0e^{-kt})}{\partial C_0}\right|_{t=t_i} = e^{-kt_i}$$Recall, that even though we just used our first concentration point as our estimate for $\hat{C}_0$, we have to propogate error in that measurement. Next we need the standard error in the residuals:
$$S_{\epsilon}^2 = \frac{\sigma_{\epsilon}^2 + \frac{1}{N}\sum_i^N \sigma_{f(x_i)}^2 + \frac{1}{N} \sum_i^N \sigma_c^2}{N - 2}$$where $\sigma_\epsilon^2 = \frac{1}{N}\sum_i^N (\hat{y} - y)^2$
The $-2$ is OK to omit, but accounts for our degrees of freedom. We need to find the variance in $f(x_i)$ term, which we already found for the weights. It is:
$$\sigma_{f(x_i)}^2 = \left(\frac{\partial C_0e^{-kt}}{\partial t}\right)^2 \sigma_{t}^2 = C_0^2k^2e^{-2kt} \sigma_t^2$$Now, combining the $\mathbf{F}$ matrix with the standard error in the residuals:
$$S_{\beta}^2 = S_{\epsilon}^2(\mathbf{F}^T\mathbf{F})^{-1}$$
In [5]:
import numpy.linalg as linalg
#Get the best fit
khat = result.x
#Now we computer error in residuals
s2_e = (np.sum((chat - c_25)**2) + np.sum(c_25_err**2) + np.sum((c_25[0] * khat * np.exp(-khat * time) * terr)**2)) / (len(time) - 2)
print s2_e
#Stack our f-matrix.
F_mat = np.column_stack( (-time * c_25[0] * np.exp(-khat * time), np.exp(-khat * time) ) )
#and now compute uncertainty
s2_beta = s2_e * linalg.inv(F_mat.transpose().dot(F_mat))
print khat, np.sqrt(s2_beta[0, 0])
khat_25 = khat
kerr_25 = np.sqrt(s2_beta[0,0])
Now, let's turn this into a function to make it easier
In [7]:
#Givens
MW = 125.0
terr = 5.0
merr = 0.5
verr = 0.25
#The SSR function to minimize
def WSSR(k, t, c, terr, cerr):
w = 1.0 / (cerr**2 + terr**2 * (k * np.exp(-k * t))**2)
chat = c[0] * np.exp(-k * t)
return np.sum(w*(chat - c)**2)
def regress_data(time, volume, mass):
#compute c and error
c = mass / volume * (1000 / MW)
c_err = np.sqrt( 1.0 / (volume)**2 * merr**2 + mass**2 / volume**4 * verr**2 ) * 1000.0 / MW
#compute the best-fit
mini_args = {'bounds': [(10**-10, 10)]}
result = scipy.optimize.basinhopping(lambda s: WSSR(s, time, c, terr, c_err), x0=0.1, minimizer_kwargs=mini_args, niter=1000)
#plot it and compute chats
chat = c[0] * np.exp(-result.x * time)
plt.plot(time, chat)
plt.errorbar(time, c, xerr=terr, yerr=c_err, capthick=3)
plt.ylabel('Concentration [mol/L]')
plt.xlabel('Time [minutes]')
plt.show()
#uncertainty analysis
khat = result.x
s2_e = (np.sum((chat - c)**2) + np.sum(c_err**2) + np.sum((c[0] * khat * np.exp(-khat * time) * terr)**2)) / (len(time) - 2)
F_mat = np.column_stack( (-time * c[0] * np.exp(-khat * time), np.exp(-khat * time) ) )
#print F_mat
s2_beta = s2_e * linalg.inv(F_mat.transpose().dot(F_mat))
return khat, np.sqrt(s2_beta[0, 0])
In [8]:
volume_30 = np.array([243.7, 262.7, 247.4, 224.4, 225.0, 298.5, 233.0, 285.7, 227.5, 243.3, 240.3, 276.1, 276.4, 228.0, 245.2, 236.0, 261.1, 243.3, 275.9, 239.5, 208.0])
mass_30 = np.array([ 12.1, 12.4, 9.0, 8.5, 7.3, 8.3, 4.8, 5.3, 3.3, 3.6, 2.6, 2.6, 2.4, 1.5, 2.4, 1.2, 2.3, 2.0, 2.2, 0.9, 0.4])
In [9]:
khat_30, kerr_30 = regress_data(time, volume_30, mass_30)
In [10]:
volume_10 = np.array([254.1, 280.5, 265.7, 233.8, 258.8, 239.6, 279.0, 194.8, 285.9, 278.6, 281.0, 302.7, 265.3, 209.9, 278.0, 213.2, 243.3, 274.4, 238.9, 259.3, 251.6])
mass_10 = np.array([ 12.8, 13.6, 12.5, 11.4, 13.4, 12.5, 12.7, 10.2, 14.0, 13.7, 13.1, 14.2, 12.8, 9.2, 13.0, 8.9, 11.7, 12.8, 10.6, 12.2, 11.1])
In [11]:
khat_m10, kerr_m10 = regress_data(time, volume_10, mass_10)
In [12]:
volume_8 = np.array([249.3, 251.0, 202.6, 291.9, 251.6, 214.1, 272.5, 257.8, 274.7, 261.7, 227.6, 260.9, 295.4, 277.4, 237.2, 260.8, 243.0, 249.6, 258.0, 260.6, 272.4])
mass_8 = np.array([ 12.8, 12.3, 9.6, 14.2, 11.8, 10.1, 12.0, 11.0, 11.4, 11.1, 9.2, 11.0, 11.4, 11.2, 9.6, 10.1, 10.2, 9.2, 9.2, 8.8, 10.5])
In [13]:
khat_8, kerr_8 = regress_data(time, volume_8, mass_8)
Now combine the results. Remember, that the temperature must be in an absolute scale
In [16]:
#Put together our data
k = np.array([khat_m10, khat_8, khat_25, khat_30])
k = np.reshape(k, (4,))
kerr = np.array([kerr_m10, kerr_8, kerr_25, kerr_30])
temp = np.array([-10, 8, 25, 30]) + 273.15
#Convert to kJ/ mol
R = 8.314 / 1000.0
#our new WSSR. Notice our weights do not have partial derivatives, since we have no dependent variable uncertainty.
def WSSR(A, E_a, temp, k, kerr):
w = 1.0 / kerr**2
khat = A * np.exp(-E_a / (R * temp) )
ssr = np.sum( w * (khat - k)**2)
return ssr
#Again, ust use bounded optimization
mini_args = {'bounds': [(10**-10, 10**10), (10, 1000000)]}
result = scipy.optimize.basinhopping(lambda x: WSSR(x[0], x[1], temp, k, kerr), x0=[1,1], minimizer_kwargs=mini_args, niter=10000)
print result
#Get estimated values for the fit
A = result.x[0]
E_a = result.x[1]
khat = A * np.exp(-E_a / (R * temp) )
#plot it
plt.errorbar(temp, k, fmt='o', yerr=kerr)
plt.plot(temp, khat)
plt.show()
In [17]:
#The math is not shown here, but it's very similar to above (same equation).
#Nootice we don't have partials, just the independent variable uncertainty.
s2_e = (np.sum((khat - k)**2) + np.sum(kerr**2)) / (len(k) - 2)#subtract 2 because we have two adjustable parameters
#stuck up the F-matrix
F_overall_mat = np.column_stack( (np.exp(-E_a / (temp * R)), -A / (temp * R) * np.exp(-E_a / (temp * R)) ) )
print F_overall_mat
print s2_e
s2_beta = s2_e * linalg.inv(F_overall_mat.transpose().dot(F_overall_mat))
#print results
print A, np.sqrt(s2_beta[0,0]) * scipy.stats.t.ppf(0.975, len(k) - 2)
print E_a, np.sqrt(s2_beta[1,1]) * scipy.stats.t.ppf(0.975, len(k) - 2)
So our final answer is $$E_a = 52 \pm 251 \textrm{kJ/mol}$$
$$A = 200,000 \pm 10,000,000 \frac{1}{\textrm{minutes}}$$For grading, the $E_a$ is all that matters. The uncertainty value is very sensitive to the value of $E_a$ and $A$, so be easy on them. the most important part is it's done correctly, the error bars are correct on the plot, and that they have near the correct $E_a$ value
In [ ]: