<jupyter_start><jupyter_text>Nonlinear Regression in [`SciPy`](https://docs.scipy.org/doc/scipy/reference/) and [R](https://www.r-project.org/about.html)We often need to find a function $y=f(x,\beta)$ of variable $x$ and $p$ unknown parameters $\beta$ which fits a given set of $n$ predictor, {$x_1,...,x_n$}, and response, {$y_1,...,y_n$}, values.The best fit can be obtained by minimizing the residual sums of squares (RSS) with respect to $\beta$,\begin{equation}\mathrm{RSS}(\beta) = \sum_{n=1}^{\infty} \left[ y_i - f(x_i,\beta) \right ] ^2.\end{equation}For nonlinear functions, numerical optimization methods are needed to minimize RSS and find the best-fit parameters, $\hat{\beta}$.The most common numerical optimization method for nonlinear regression is the [Gauss-Newton method](https://en.wikipedia.org/wiki/Gauss%E2%80%93Newton_algorithm).The assumption in nonlinear regression is that the measurement errors are standard normally distributed.A summary measure of the model fit is the residual variance which is defined as the ratio of the minimum of RSS devided by the number of degrees of freedom:\begin{equation}\chi_R^2 = \frac{\mathrm{RSS}(\hat{\beta})}{n-p}.\end{equation}Sometimes, we have several curves that we want to fit with the same function $f$. These curves have some parameters in common that we want to be determined on all the available data.In this Notebook, we show how to perform least-sqares fits in `scipy` and `R`.<jupyter_code>from scipy.optimize import least_squares
from scipy.optimize import curve_fit
import pandas as pd
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
plt.rcParams.update({'figure.figsize': [3.5, 3], 'font.size':8, 'legend.fontsize':8, 'figure.dpi': 150,
                     'legend.handlelength':1, 'legend.frameon':False,
                     'xtick.major.size':6,'ytick.major.size':6,
                     'xtick.minor.size':3,'ytick.minor.size':3,
                     'xtick.direction':'out','ytick.direction':'out',
                     'xtick.labelsize':8,'ytick.labelsize':8,
                     'xtick.major.pad':2,'ytick.major.pad':2,
                     'axes.linewidth':1,'axes.labelsize':10,})
from cryptography.fernet import Fernet
import base64
def encrypt(string, key):
    keygen = lambda x : base64.urlsafe_b64encode(x.encode() + b' '*(32 - len(x)))
    cipher = Fernet(keygen(key))
    return cipher.encrypt(string.encode())
def decrypt(string, key):
    keygen = lambda x : base64.urlsafe_b64encode(x.encode() + b' '*(32 - len(x)))
    cipher = Fernet(keygen(key))
    return print(cipher.decrypt(string.encode()).decode())<jupyter_output><empty_output><jupyter_text>Case Study: Time-Resolved Fluorescence Spectroscopy SDS Micelles in the Presence of a Hydrophobic Fluorophore and a Fluorescence Quencher- **The micelle:**    Sodium dodecyl sulfate (SDS) is an anionic surfactant with a hydrophobic tail of 12 carbon atoms and a hydrophilic sulphate head group. Due to its amphipatic nature, SDS self-assembles into micelles in aqueous solution. Micelle formation is a cooperative process which occurs in a narrow range of surfactant concentrations around the critical micellar concentration (CMC).At concentrations below the CMC., surfactant molecules are mostly in the monomer state.At the  (CMC), surfactant molecules form micelles with a well-defined aggregation number. The size distribution of the aggregates is rather narrow and depends on the relative dimensions of the head and tail groups.Above the CMC, the number of micelles increases while the size of the micelles stays the same.       - **The fluorophore:**Pyrene is a hydrophobic fluorescent molecule.In an aqueous SDS solution at [SDS]>CMC, pyrene tends to partition into the hydrophobic interior of the micelles.- **The quencher:**Benzophenone is a hydrophobic compound which also tend to be in the micelles where it quenches the fluorescence of pyrene (it makes it decay faster). - **The experiment:**We measure the number of photons emitted by pyrene molecules in SDS micelles for increasing time intervals after the excitation pulse. The fluorescence intensity decay of pyrene following an excitation pulse is given by\begin{equation}F(t)=F(0)\exp(-t/\tau_0)\end{equation}where $\tau_0$ is the lifetime of the fluorophore.In the presence of $j$ benzophenone molecules per micelle, the lifetime of pyrene decreases to\begin{equation}\tau_j = \frac{\tau_0}{1 + j k_q \tau_0}.\end{equation}where $k_q$ is the quenching constant.The distribution of quenchers in the micelles follows Poisson's statistics, therefore the fraction of micelles containing $j$ benzophenone molecules is given by\begin{equation}P(j) = \frac{m^j}{j!}\exp{(-m)}\end{equation}where $m$ is the average number of benzophenone molecules per micelle.The fluorescence intensity decay of pyrene in the presence of benzophenone is given by the weighted average over the difstribution of micelles with different numbers of quenchers\begin{equation}F(t)=F(0)\sum_{j=0}^{\infty} P(j) \exp{(-t/\tau_j)}=F(0)\exp\left\{-t/\tau_0-m\left [1-\exp{(-k_qt)}\right ]\right\},\end{equation}where we have used the converging series\begin{equation}\sum_{j=0}^{\infty} \frac{x^j}{j!} = \exp{(x)}.\end{equation}We used a solution 70 mM SDS and 10 $\mu$M pyrene in 5mM phosphate buffer at pH 7.4. The emission wavelength of the spectrometer is set to 391 nm, corresponding to the maximum fluorescence intensity for pyrene.- First, we measure the intensity decay of the solution in the absence of quencher. A least-squares fit on this curve yields $F(0)$ and the pyrene lifetime, $\tau_0$. - Then, we measure intensity decays for increasing concentrations of benzophenone, _i.e._, $[$quencher$]$ $=$ 0.1, 0.2, 0.5, 1 e 1.5 mM.We perform a global fit on these five curves. Each corve has a different average number of quencher molecules per micelle, $m$, and common parameters $F(0)$, $\tau_0$, and $k_q$. Since the excitation pulse occurs in a shorter time scale compared to the decay of the fluorescence signal, we disregard the width of the excitation pulse in our analysis. In the least-squares fitting procedures, we weight each photon count, $F(t)$, by its error which is simply $\sqrt{F(t)}$ because fluorescence counts follow Poisson's statistics.  The DataThe data are tab-separated values (TSV) files (collected at [Prof. L. Stella lab](http://www.stc.uniroma2.it/physchem/stella.html))<jupyter_code>%cat 'data/fluorescence/quencher_0.0.txt'
%ls 'data/fluorescence'<jupyter_output><empty_output><jupyter_text>Let's plot the intensity decays for increasing concentrations of quencher.The intensity signal decays faster with increasing concentration of benzophenone.<jupyter_code>concs = [0.0,0.2,0.5,1.0,1.5]
for conc in concs:
    t,F = np.loadtxt('data/fluorescence/quencher_{:1.1f}.txt'.format(conc), unpack=True, skiprows=9)
    plt.errorbar(t, F, np.sqrt(F), lw=0, marker='o', label='{:1.1f}'.format(conc)+' mM',
                 ms=2, elinewidth=1., capsize=2, capthick=1., mfc='w')
plt.legend(title='[quencher]',labelspacing=.7,loc=1)
plt.ylabel('$F(t)$  /  counts'); plt.xlabel('$t$  /  ns');<jupyter_output><empty_output><jupyter_text>Here, we use boolean masking to select the data points after the excitation pulse and we shift the x values so that they measure the time intervals from the excitation pulse.<jupyter_code>for conc in concs:
    d = np.loadtxt('data/fluorescence/quencher_{:1.1f}.txt'.format(conc), unpack=False, skiprows=9)
    d = d[ d[:,0] > d[d[:,1]==d[:,1].max(), 0] ] # points before the intesity max are discarded
    plt.errorbar(d[:,0]-d[0,0], d[:,1], np.sqrt(d[:,1]), lw=0, marker='o', label='{:1.1f}'.format(conc)+' mM',
                 ms=2, elinewidth=1., capsize=2, capthick=1., mfc='w')
plt.legend(title='[quencher]',labelspacing=.7)
plt.ylabel('$F(t)$ / counts'); plt.xlabel('$t$ / ns'); plt.show()<jupyter_output><empty_output><jupyter_text>A convenient way to save data sets to file is to use [pickling](https://docs.python.org/3/library/pickle.html). We create a  `pandas` DataFrame and we save it as a binary file and also as a csv file.<jupyter_code>data = pd.DataFrame()
for conc in concs:
    t,F = np.loadtxt('data/fluorescence/quencher_{:1.1f}.txt'.format(conc), unpack=True, skiprows=9)
    # points before the intesity max are discarded
    t = t[F.argmax():]
    F = F[F.argmax():]
    t = t - t[0] # t at the intesity max is t zero
    d = pd.DataFrame(data=np.c_[t,F,np.sqrt(F)],columns=['t','F','errF'])
    d['conc'] = conc
    d.set_index('conc',inplace=True)
    data = data.append(d)
data.to_pickle('data/fluorescence.p')
data.to_csv('data/fluorescence.csv')
data.head()<jupyter_output><empty_output><jupyter_text>The curve for $[$quencher$]$=0 mM is fitted to a single exponential decay function using [`optimize.least_square`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html) which requires the following arguments:- the function which computes the vector of residuals- the initial parameters- the data- the lower and the upper bounds of the fitting parameters`least_squares(function_to_minimize, initial_params,     [x, y, yerr], bounds=[(lower_bounds),(upper_bounds)])``least_squares` returns an OptimizeResult object containig the best-fit parameters (`x`), the residuals calculated for the best-fit parameters (`fun`), RSS/2 (`cost`), and the modified Jacobian matrix, $J$, at the solution (`jac`), $J^T J$ is a [Gauss-Newton approximation of the Hessian](https://reference.wolfram.com/language/tutorial/UnconstrainedOptimizationGaussNewtonMethods.html) of the residual function.Since we know that the error on the counts follows Poisson's distribution, we use the following [nonlinear regression model](https://link.springer.com/book/10.1007%2F978-0-387-09616-2)\begin{equation}y_i = f(x_i,\beta) + \frac{\epsilon_i}{\sqrt{w_i}},\end{equation}where the weight is the inverse of the error on the counts, $1 / \sqrt{y_i}$.We minimize, with respect to $\beta$, the residual sums of squares (RSS) weighted by $y_i$:\begin{equation}\mathrm{RSS}(\beta) = \sum_{n=1}^{\infty} \frac{\left[ y_i - f(x_i,\beta) \right ] ^2}{y_i}.\end{equation}For further details read chapter 6 of [Nonlinear Regression with R](https://link.springer.com/book/10.1007%2F978-0-387-09616-2) by <NAME> and <NAME>.<jupyter_code># we define the exponential decay
def exp_decay(x, param):
    return param[0] * np.exp( - x / param[1] ) 

# we define the function that calculates the residuals
def residuals(param, *args):
    return ( args[1] - exp_decay(args[0], param) ) / args[2] # residuals weighted by the error on y
# we create a list with the data
args = [data.loc[0.0][col] for col in ['t','F','errF']]

# we fit the data
res = least_squares(residuals, x0=[3000,150], args=args, bounds=[(6,150),(3500,200)])

# we plot the data
plt.errorbar(args[0],args[1],args[2],lw=0,marker='o',ms=1,elinewidth=1.,capsize=3,capthick=1.,
             color='k',alpha=0.3)

# we plot the fitting curve
plt.plot(args[0],exp_decay(args[0],res['x']),lw=3,color='r')

# we set the log scale for the y-axis
plt.yscale('log')
plt.ylabel('$\ln [F(t)]$  /  counts'); plt.xlabel('$t$  /  ns'); 
print('Is one of the convergence criteria satisfied?', res.success) 
print('The best-fit parameters are',res.x[0],'and',res.x[1])
from IPython import display

table = r"""| Parameter     | Value         | Error |
| :---------------- |:-------------:| :----:|
| $F(0)$ (counts)   | %.2f          | ?     |
| $\tau_0$ (ns)     | %.2f          | ?     |"""

display.Markdown(table % tuple(res.x) )<jupyter_output><empty_output><jupyter_text>\begin{equation}\chi_R^2 = \frac{\mathrm{RSS}(\hat{\beta})}{n-p}\end{equation}<jupyter_code>print('The residual variance is', 2 * res.cost / (res.fun.size - res.x.size))<jupyter_output><empty_output><jupyter_text>We can check whether the model is appropriate and we have obtained a good fit by comparing the Gaussian with $\sigma=1$ and $\mu=0$ with the histogram of the residuals at solution.We can see that the distribution of the residuals is compatible with a normal distribution. Would the comparison be this favorable if we didn't weight the error by the square root of the counts?<jupyter_code>bins = np.linspace(res.fun.min(),res.fun.max()+1,100)
plt.plot(bins, np.exp(-bins*bins/2)/np.sqrt(2*np.pi), label='Gaussian', color='k')
plt.hist(res.fun, bins, lw=2, histtype='step', density=True, label='residuals', color='r')
plt.xlabel('Residuals'); plt.ylabel('Probability Distribution'); plt.legend();<jupyter_output><empty_output><jupyter_text>The uncertainties on the best-fit parameters are estimated from the variance-covariance matrix.Here is the implementation in the function [`optimize.curve_fit`)](https://github.com/scipy/scipy/blob/2526df72e5d4ca8bad6e2f4b3cbdfbc33e805865/scipy/optimize/minpack.pyL739) which uses `least_squares` under the hood.Further information [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.htmlscipy.optimize.curve_fit) and [here](http://scipy-cookbook.readthedocs.io/items/FittingData.html).<jupyter_code>from scipy.linalg import svd
# copied from curve_fit 
def error_analysis(res):
    rss = 2 * res.cost  # res.cost is half sum of squares!
    s_sq = rss / (res.fun.size - res.x.size) # residual variance or reduced chi squared
    # Do Moore-Penrose inverse discarding zero singular values.
    _, s, VT = svd(res.jac, full_matrices=False)
    threshold = np.finfo(float).eps * max(res.jac.shape) * s[0]
    s = s[s > threshold]
    VT = VT[:s.size]
    pcov = np.dot(VT.T / s**2, VT) * s_sq # variance-covariance matrix
    return np.sqrt(np.diag(pcov)), s_sq # errors on the params are the diagonal elements of pcov
table = r"""| Parameter     | Value         | Error |
| :---------------- |:-------------:| :----:|
| $F(0)$ (counts)   | %.0f          | %.0f  |
| $\tau_0$ (ns)     | %.1f          | %.1f  |"""
p_and_err = np.repeat(res.x, 2)
p_and_err[1::2] = error_analysis(res)[0]
print('Reduced chi squared:',error_analysis(res)[1])
display.Markdown(table % tuple(p_and_err) )<jupyter_output><empty_output><jupyter_text>Task 1- Fit the data for $[$quencher$]$ $=$ 0 using [`optimize.curve_fit`)](https://github.com/scipy/scipy/blob/2526df72e5d4ca8bad6e2f4b3cbdfbc33e805865/scipy/optimize/minpack.pyL739).- Plot the data and the function with the best-fit parameters.- View a markdown table with the best-fit parameters and their uncertainties.<jupyter_code>answer = '<KEY>

decrypt(answer,)<jupyter_output><empty_output><jupyter_text>Interfacing Python and [R](https://www.r-project.org/about.html)This part of the Notebook requires the packages r-essentials and rpy2. To install them, issue the following commands in the terminal```bashsource activate LUcomputeconda install -c conda-forge r-essentialsconda install -c conda-forge rpy2```It is also possible to install rpy2 using pip<jupyter_code>!pip install rpy2
We load the extension [rpy2](https://rpy2.readthedocs.io/en/version_2.8.x/)
%load_ext rpy2.ipython
X = args[0]
Y = args[1]
errY = args[2]<jupyter_output><empty_output><jupyter_text>We send the data from Python to R<jupyter_code>%Rpush X Y errY<jupyter_output><empty_output><jupyter_text>We plot the data in R<jupyter_code>%%R
plot(X, Y, ylim = range(c(Y-errY, Y+errY)), xlab='t  /  ns', ylab='ln [F(t)]  /  counts', log='y')
arrows(X, Y-errY, X, Y+errY, length=0.05, angle=90, code=3)<jupyter_output><empty_output><jupyter_text>We create a data frame in R<jupyter_code>%%R
df <- data.frame('X' = X, 'Y' = Y)<jupyter_output><empty_output><jupyter_text>We fit the data without weighting by the error<jupyter_code>%%R
fit1 <- nls(Y ~ F0 * exp( - X / tau0 ), data=df, start = list(F0 = 2000, tau0 = 100))
summary(fit1)<jupyter_output><empty_output><jupyter_text>We weight by the error and we compare the results of the two approaches<jupyter_code>%%R
plot(Y ~ X, data=df, xlab='t  /  ns', ylab='ln [F(t)]  /  counts', log='y' )

fit2 <- nls(Y ~ F0 * exp( - X / tau0 ), data=df, start = list(F0 = 2000, tau0 = 100), weights=1/Y)

curve(predict(fit2, newdata = data.frame(X=x)), add=TRUE, lwd=4, col='red')

curve(predict(fit1, newdata = data.frame(X=x)), add=TRUE, lwd=4, col='blue')

legend(400, 1500, legend=c('weighting','w/o weighting'),
       col=c("red", "blue"), lwd=4, cex=1.5)

summary(fit2)<jupyter_output><empty_output><jupyter_text>We calculate the residual variance<jupyter_code>%%R
sum(residuals(fit2))/952<jupyter_output><empty_output><jupyter_text>Global Least-Squares FitNow we performa a global fit of the 5 curves with increasing [quencher]: $F(0)$, $\tau_0$, and $k_q$ are global parameters whereas $m$ depends on [quencher].* We create 3 lists (for x, y, and yerr) each containing 5 `NumPy` arrays (for the 5 curves).* We create a list called args containing the 3 lists (allx, ally, allyerr)<jupyter_code># args is a list of three lists containing 5 numpy arrays
args = [ [data.loc[conc][col] for conc in concs] for col in ['t','F','errF'] ] 
param0 = [3000, 160, 0.04, 0] + [0.5]*4 # initial parameters
lbound = [0, 0, 0] + [0]*5              # lower bound
ubound = [4000, 200, 1, 1e-5] + [3]*4   # upper bound<jupyter_output><empty_output><jupyter_text>In the residuals function, we calculate the vector of residuals for each curve which appended to a single 1D `NumPy` array.<jupyter_code># we define the function for the intensity decay
def decay(x, param):
    return param[0]*np.exp( - x/param[1] + param[3]*np.expm1(-param[2]*x) ) 
# we define the function that calculates the least squares
def residuals(p, *args):
    residuals = np.empty(0)
    for i,m in enumerate(p[3:]):
        param = p[0], p[1], p[2], m
        r = ( args[1][i] - decay(args[0][i], param) ) / args[2][i]
        residuals = np.append(residuals, r)
    return residuals
res = least_squares(residuals, x0=param0, args=args, bounds=[lbound,ubound])
c = plt.rcParams['axes.prop_cycle'].by_key()['color']
p = res['x']
for i,m in enumerate(p[3:]):
        param = p[0], p[1], p[2], m
        plt.fill_between(args[0][i],args[1][i]-args[2][i],args[1][i]+args[2][i],alpha=.7)
        plt.plot(args[0][i],decay(args[0][i],param),lw=1,color=c[i],label='{:1.1f}'.format(concs[i])+' mM')
plt.yscale('log'); plt.legend(title='[quencher]',labelspacing=.7)
plt.ylabel('$\ln F(t)$ / counts'); plt.xlabel('$t$ / ns'); 
print('Reduced chi squared:',error_analysis(res)[1])<jupyter_output><empty_output><jupyter_text>Instead of using markdown to view the results in a table, we create a `pandas` DataFrame.<jupyter_code>table = pd.DataFrame({'Value':res.x,'Error':error_analysis(res)[0]},
             index=['$F(0)$ (counts)','$\tau_0$ (ns)','$k_q$ (ns$^{-1}$)','$m_0$','$m_1$','$m_2$','$m_3$','$m_4$'],
            columns=['Value','Error'])
table
print(table.to_latex())<jupyter_output><empty_output><jupyter_text>Task 2 Take the logarithm of each curve and perform a global fit to $\ln F(t) = \ln F(0) - t/\tau_0 + m [\exp(-k_q t) -1]$ using `least_square`.1. Modify the `NumPy` arrays in args[1] by taking the `np.log`.```.pydata = pd.read_pickle('data/fluorescence.p')allx = []; ally = []; allyerr = []  three empty listsfor conc in data.columns:    allx.append(data[conc]['t']); ally.append(data[conc]['F']); allyerr.append(data[conc]['Ferr'])args = allx, ally, allyerr  args is a list of three lists containing 5 numpy arraysfor i in range( data.columns.size ):    args[1][i] = ?    args[2][i] = ? ```2. Modify the `NumPy` arrays in args[2] according to the error propagation rule: error$(\ln F)$ = error($F$) / $F$.3. Define a new fitting function for the linearized equation. 4. Fitting and plotting can be carried out exactly as in the previous cell.<jupyter_code>answer = '<KEY>'

decrypt(answer,)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1491494532199271, 'fraction_numerical': 0.02303564196030782, 'mean_word_length': 4.730490281404119, 'pattern_counts': {'":': 0, '<': 92, '<?xml version=': 0, '>': 92, 'https://': 13, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '58885', 'license': 'MIT', 'n_tokens_mistral': 7132, 'n_tokens_neox': 6599, 'n_words': 1948, 'path': 'lectures/least_squares.ipynb', 'repo_name': 'urania277/jupyter-course', 'size': '24055'}