<jupyter_start><jupyter_text>Model fittingA lot of chemical analysis relies on our ability to fit a "line of best fit" to some data. For example, considering the determination of the activation energy, $E_a$, of some reaction from an Arrhenius relationship, $$ k(T) = A \exp{\frac{-E_a}{RT}}, $$where, $k(T)$ is the rate constant at time $T$, $A$ is the pre-exponential factor, $R$ is the ideal gas constant and $T$ is the temperature. Typically when we are performing our analysis, we will rearrange the above equation to give a straight line, $$ \ln{k(T)} = \ln{A} - \frac{E_a}{RT}. $$Let's looks at some experimental data for the decomposition of 2 moles of NO2 to 2 moles of NO and one mole of O2.The data for this is available [in this file](https://raw.githubusercontent.com/pythoninchemistry/ch40208/master/CH40208/working_with_data/arrhenius.txt).<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
T, k = np.loadtxt('arrhenius.txt', unpack=True)
plt.plot(T, k, 'o--')
plt.xlabel('Temperature/K')
plt.ylabel('Rate constant/M$^{-1}$s$^{-1}$')
plt.show()<jupyter_output><empty_output><jupyter_text>Alternatively, we plot this as a straight line.<jupyter_code>plt.plot(1/T, np.log(k), 'o--')
plt.xlabel('1/Temperature/K')
plt.ylabel('ln(Rate constant/M$^{-1}$s$^{-1}$)')
plt.show()<jupyter_output><empty_output><jupyter_text>Using Python, we can fit this data in either configuration, either by performing a [linear regression](https://en.wikipedia.org/wiki/Linear_regression) or by fitting a line of best fit. Linear regressionLinear regression is a relatively simple mathematical process, that is implemented in the `scipy.stats` module.<jupyter_code>from scipy.stats import linregress<jupyter_output><empty_output><jupyter_text>If we investigate the [documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.linregress.html) for the `linregress` function, we will see that we pass two arrays; `x` and `y`. These corresponding exactly to the *x* and *y*-axes of the plot above.<jupyter_code>result = linregress(1/T, np.log(k))
result<jupyter_output><empty_output><jupyter_text>We can see that the object that is returned contains a couple of pieces of information. The two parts that we are most interested in are the `slope` and the `intercept`, which are attributes of the `result` object and can be obtained as follows.<jupyter_code>print(result.slope)
print(result.intercept)<jupyter_output>15.161037389899018<jupyter_text>Therefore, we can evaluate the activation energy and the pre-exponential factor as follows.<jupyter_code>from scipy.constants import R
activation_energy = - R * result.slope
preexponential_factor = np.exp(result.intercept)
print(f"A = {preexponential_factor}; E_a = {activation_energy}")<jupyter_output>A = 3840209.1525825486; E_a = 112323.03523328649<jupyter_text>Finally, we can plot our linear regression line through the data as shown below.<jupyter_code>x = np.linspace(1/660, 1/590, 1000)

plt.plot(1/T, np.log(k), 'o')
plt.plot(x, result.slope * x + result.intercept, '-')
plt.xlabel('1/Temperature/K')
plt.ylabel('ln(Rate constant/M$^{-1}$s$^{-1}$)')
plt.show()<jupyter_output><empty_output><jupyter_text>Non-linear regressionOf course, our equation cannot always be modified to offer a straight line formulation or we may not want to as the logarithmic scaling of the *y*-axis may skew our understanding of the intercept. Therefore, necessary to fit a non-linear relationship. This can also be achieved with Python by minimising the difference between our model (the equation with given parameters) and the data. The value that we will minimize is known as the $\chi^2$-value, which has the following definitions.$$ \chi^2 = \sum\big\{[y - f(x)] ^ 2\big\} \;\;\;\text{or}\;\;\; \chi^2 = \sum\Bigg\{\frac{[y - f(x)] ^ 2}{y_{\text{err}}^2}\Bigg\},$$depending on whether or not uncertainty in the data is known. Consider the following data, which describe the decomposition of hydrogen peroxide in the presence of excess cerous ion CeIII,2H2O2 + 2CeIII &rarr; 2H2O + CeIV + O2.| Time/s | 2 | 4 | 6 | 8 | 10 | 12 | 14 | 16 | 18 | 20 ||---|---|---|---|---|---|---|---|---|---|---|| [H2O2]/moldm-3 | 6.23 | 4.84 | 3.76 | 3.20 | 2.60 | 2.16 | 1.85 | 1.49 | 1.27 | 1.01 |First, we will store these values as NumPy arrays.<jupyter_code>c = np.array([6.23, 4.84, 3.76, 3.20, 2.60, 2.16, 1.85, 1.49, 1.27, 1.01])
t = np.array([2, 4, 6, 8, 10, 12, 14, 16, 18, 20])
plt.plot(t, c, 'o')
plt.xlabel('Time/s')
plt.ylabel('[H$_2$O$_2$]/M')
plt.show()<jupyter_output><empty_output><jupyter_text>This reaction has a first-order rate equation, this will be our "model" that we will compare to the data, where the $[A]_0$ and $k$ are the varying parameters. This function is defined below.<jupyter_code>def first_order(t, a0, k):
    """
    The first-order rate equation. 
    
    Args:
        t (float): Time (s).
        a0 (float): Initial concentration (mol/dm3).
        k (float): Rate constant (s-1).
    
    Returns:
        (float): Concentration at time t (mol/dm3).
    """
    return a0 * np.exp(-k * t)<jupyter_output><empty_output><jupyter_text>We want to vary `a0` and `k` to minimise the $\chi^2$-value between this model and the data.<jupyter_code>def chi_squared(parameters, t, conc):
    """
    Determine the chi-squared value for a first-order rate equation.
    
    Args:
        parameters (list): The variable parameters of first_order (a0 and k).
        t (float): The experimental time (s) data.
        conc (float): The experimental concentration data which first_order aims to model.
    
    Returns:
        (float): chi^2 value.
    """
    a0 = parameters[0]
    k = parameters[1]
    return np.sum((conc - first_order(t, a0, k)) ** 2)<jupyter_output><empty_output><jupyter_text>Just like the minimisation of a potential model.We can use the `scipy.optimize.minimize` to minimize this chi-squared function.<jupyter_code>from scipy.optimize import minimize
result = minimize(chi_squared, [10, 1], args=(t, c))
result
x = np.linspace(2, 20, 1000)

plt.plot(t, c, 'o')
plt.plot(x, first_order(x, result.x[0], result.x[1]))
plt.xlabel('Time/s')
plt.ylabel('[H$_2$O$_2$]/M')
plt.show()
print(f"[A]_0 = {result.x[0]}; k = {result.x[1]}")<jupyter_output>[A]_0 = 7.421644144148406; k = 0.10308255579543868<jupyter_text>`scipy.optimize.curve_fit` offers a convenient wrapper around this methodology, and makes it easy to access the estimated parameter uncertainties.<jupyter_code>from scipy.optimize import curve_fit
popt, pcov = curve_fit(first_order, t, c)<jupyter_output>/tmp/ipykernel_30479/2319394663.py:13: RuntimeWarning: overflow encountered in exp
  return a0 * np.exp(-k * t)<jupyter_text>We can ignore the above error, this is just warning that during the optimization an unphysical value was found.However, the `curve_fit` function is robust enough to ignore this and continue.<jupyter_code>uncertainties = np.sqrt(pcov)

print(f"[A]_0 = {popt[0]} +/- {uncertainties[0][0]}; k = {popt[1]} +/- {uncertainties[1][1]}")<jupyter_output>[A]_0 = 7.42164476326434 +/- 0.15838876657483184; k = 0.10308257141927114 +/- 0.00308836766268339
{'fraction_non_alphanumeric': 0.13900967877682704, 'fraction_numerical': 0.05063823818207322, 'mean_word_length': 4.627466456195738, 'pattern_counts': {'":': 0, '<': 48, '<?xml version=': 0, '>': 48, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '390563', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 2703, 'n_tokens_neox': 2394, 'n_words': 762, 'path': 'CH40208/working_with_data/curve_fitting.ipynb', 'repo_name': 'pythoninchemistry/ch40208', 'size': '7129'}