\[\newcommand{\arr}[1]{\underline{\underline{#1}}}\]

\[\newcommand{\vec}[1]{\underline{#1}}\]

\[\require{mhchem}\]

Non-linear regression#

Recap on linear regression#

Last class, we talked about how we could turn linear regression into a linear algebra problem

You can calculate this yourself
You can also do linear regression using statistical packages like statsmodel

Today we will discuss two ways of solving non-linear regression problems

Turn a non-linear problem into a linear one and solve
Non-linear curve fitting

Turn a non-linear regression problem into a linear regression problem#

Rate constants and reaction orders are determined by using models that are fit to experimental data
A common case is to monitor concentration vs. time in a constant volume, batch reactor
We consider the disappearance of $A$
From the mole balance we know:

(507)#\[\begin{equation} \frac{dN_A}{dt} = r_A V \end{equation}\]

Let us assume the rate law is of the form: $r_A = k C_A^\alpha$ and a constant volume so that:

(508)#\[\begin{equation} \frac{dC_A}{dt} = -k C_A^\alpha \end{equation}\]

Let us be loose with mathematics, rearrange the equation, and take the log of both sides.
- By loose I mean we take logs of quantities that are not dimensionless

(509)#\[\begin{equation} \ln(-\frac{dC_A}{dt}) = \ln{k} + \alpha \ln C_A \end{equation}\]

This suggests that if we could numerically compute $\frac{dC_A}{dt}$ from our data of $C_A(t)$ then a plot of the log of the negative derivative vs the log of concentration would have
- an intercept equal to the log of the rate constant, $k$
- and a slope equal to the reaction order $\alpha$
Given the following data, determine the reaction order in A and the rate constant with 95% confidence intervals.

time (min)	C\_A (mol/L)
0	0.0500
50	0.0380
100	0.0306
150	0.0256
200	0.0222
250	0.0195
300	0.0174

We can get the derivatives by first fitting a spline through the data. The spline is essentially just a smoothing function
We will use the splev function to numerically compute derivatives from spline fits of the function.
This works best when the $x$ points are evenly spaced, and they should be monotically increasing or decreasing

import numpy as np
import matplotlib.pyplot as plt

data=np.array([[0,0.05],
               [50,.038],
               [100,.0306],
               [150,.0256],
               [200,.0222],
               [250,.0195],
               [300,.0174]])

plt.plot(data[:,0],data[:,1],'o')
plt.xlabel('Time [min]')
plt.ylabel('Conc. [mol/L]')

Text(0, 0.5, 'Conc. [mol/L]')

../../_images/26-bonus-nonlinear-regression_5_1.png

So, we need to convert the list of numbers to a numpy array so we can do the analysis.

import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate

# data will be a 2d list, which we convert to an array here
data = np.array(data)
t = data[:, 0]   # column 0
Ca = data[:, 1]  # column 1

# calculate a spline through the data
tck = interpolate.splrep(t, Ca)

t_eval = np.linspace(0,300)
Ca_spline = interpolate.splev(t_eval, tck)
plt.plot(data[:,0],data[:,1],'o', label='Exp Data')
plt.plot(t_eval, Ca_spline,'--.',label='Spline Fit')
plt.xlabel('Time [min]')
plt.ylabel('Conc. [mol/L]')
plt.legend()

<matplotlib.legend.Legend at 0x7f001ee04cd0>

../../_images/26-bonus-nonlinear-regression_7_1.png

import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
import statsmodels.api as sm

# data will be a 2d list, which we convert to an array here
data = np.array(data)
t = data[:, 0]   # column 0
Ca = data[:, 1]  # column 1

# calculate numerical derivatives
tck = interpolate.splrep(t, Ca)
dCadt = interpolate.splev(t, tck, der=1)

# do the transformation
x = np.log(Ca)
y = np.log(-dCadt)

# setup and do the regression
# column of ones and x:  y = b + mx
X = np.column_stack([x, x**0])

mod = sm.OLS(y, X)
res = mod.fit()

intercept = res.params[1]
alpha = res.params[0]

confidence_intervals = res.conf_int(0.05)
intercept_error = confidence_intervals[1]
alpha_error = confidence_intervals[0]


print('alpha = {0}, conf interval {1}'.format(alpha, alpha_error))
print('k = {0}, conf interval {1}'.format(np.exp(intercept), 
                                          np.exp(intercept_error)))


# always visually inspect the fit
plt.plot(x, y,'o')
plt.plot(x, res.predict(X))
plt.xlabel('$\ln(C_A)$')
plt.ylabel('$\ln(-dC_A/dt)$')
plt.show()

alpha = 2.0354816446001145, conf interval [1.92418422 2.14677907]
k = 0.1402128334966662, conf interval [0.09372748 0.20975319]

../../_images/26-bonus-nonlinear-regression_8_1.png

res.summary()

/home/runner/micromamba-root/envs/buildenv/lib/python3.10/site-packages/statsmodels/stats/stattools.py:74: ValueWarning: omni_normtest is not valid with less than 8 observations; 7 samples were given.
  warn("omni_normtest is not valid with less than 8 observations; %i "

OLS Regression Results
Dep. Variable:	y	R-squared:	0.998
Model:	OLS	Adj. R-squared:	0.997
Method:	Least Squares	F-statistic:	2210.
Date:	Wed, 19 Oct 2022	Prob (F-statistic):	8.22e-08
Time:	19:29:39	Log-Likelihood:	13.785
No. Observations:	7	AIC:	-23.57
Df Residuals:	5	BIC:	-23.68
Df Model:	1
Covariance Type:	nonrobust

	coef	std err	t	P>\|t\|	[0.025	0.975]
x1	2.0355	0.043	47.013	0.000	1.924	2.147
const	-1.9646	0.157	-12.539	0.000	-2.367	-1.562

Omnibus:	nan	Durbin-Watson:	2.377
Prob(Omnibus):	nan	Jarque-Bera (JB):	0.735
Skew:	-0.181	Prob(JB):	0.692
Kurtosis:	1.454	Cond. No.	40.4

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

You can see there is a reasonably large range of values for the rate constant and reaction order (although the confidence interval does not contain zero)
The fit looks ok, but you can see the errors are not exactly random
- There seems to be systematic trends in a sigmoidal shape of the data
- That suggests small inadequacy in the model
Let us examine some methods of evaluating the quality of fit
First we examine the residuals, or the errors between the data and the model.
In a good fit, these will be randomly distributed
In a less good fit, there will be trends

residuals = y - res.predict(X)

# always visually inspect the fit
plt.plot(x, residuals, 'o-')
plt.xlabel('$\ln(C_A)$')
plt.ylabel('residuals')
plt.show()

You can see there are trends in this data
- That means the model may not be complete
There is uncertainty in the data
- In each concentration measurement there is uncertainty in the time and value of concentration
- You need more data to reduce the uncertainty
- You may also need better data to reduce the uncertainty
Derivatives tend to magnify errors in data
- The method we used to fit the data contributed to the uncertainty
We also nonlinearly transformed the errors by taking logs and exp of the data and results, which may have skewed the confidence limits