{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Guide to Data Fitting in Python\n",
    "\n",
    "<div class=\"overview admonition\"> \n",
    "<p class=\"admonition-title\">Overview</p>\n",
    "\n",
    "Questions:\n",
    "\n",
    "* How can I use Python to fit a linear model?\n",
    "\n",
    "* How can I use Python to fit arbitrary functional forms to data?\n",
    "\n",
    "* What libraries can be used for data fitting, and what are the differences between them?\n",
    "\n",
    "Objectives:\n",
    "\n",
    "* Introduce different libraries that can be used for data fitting and the differences between them.\n",
    "\n",
    "* Fit a model using NumPy polyfit.\n",
    "\n",
    "* Fit a model using SciPy linregress.\n",
    "\n",
    "* Fit a model using statsmodels.\n",
    "\n",
    "* Fit a model using scikit-learn.\n",
    "\n",
    "* Fit an arbitrary functional form to data using SciPy curve_fit.\n",
    "\n",
    "</div>\n",
    "\n",
    "## Purpose of this lesson\n",
    "This lesson will demonstrate how to fit models to data using Python.\n",
    "We will explore different libraries that can be used for linear fitting and show the different ways to fit a linear model using these libraries.\n",
    "This lesson is aimed at educators who want to teach students how to fit models to data using Python, and its purpose is to provide an overview of the different libraries available for data fitting, the differences between them, and what kind of fitting can be done with each option..\n",
    "\n",
    "## Choosing the Right Python Library for Fitting\n",
    "\n",
    "In this section, we will explore different Python libraries that can be used for fitting models.\n",
    "As with most tasks when programming, there are many ways to achieve the same goal.\n",
    "When it comes to fitting a model, you have several options.\n",
    "The approach you choose will depend on the complexity of the code you'd like to write and the statistical information you need. \n",
    "\n",
    "Below are some common Python libraries used for linear fitting.\n",
    "They are generally ordered from simplest to most complex in terms of use and the statistical information they provide.\n",
    "\n",
    "### NumPy (`polyfit`)\n",
    "\n",
    "- **What is it?**: NumPy is a fundamental package for numerical computing in Python. The [`np.polyfit`](https://numpy.org/doc/stable/reference/generated/numpy.polyfit.html) function can be used to fit polynomial functions of specified degree. \n",
    "- **When to use it?**: Use the `np.polyfit` function when you need a simple linear fit or if you are doing a polynomial fit. It can be especially useful if you are already working with NumPy and do not want to add extra dependencies to your lesson. \n",
    "\n",
    "### SciPy (`linregress`)\n",
    "\n",
    "- **What is it?**: SciPy is a Python library used for scientific and technical computing. The `linregress` function from SciPy provides a quick way to perform linear regression. It returns important statistics such as the slope, intercept, correlation coefficient, p-value, and standard error of the regression, which can be helpful for straightforward analyses.\n",
    "- **When to use it?**: Use `linregress` when you need a fast way to perform linear regression with some basic statistical insights. \n",
    "\n",
    "### Statsmodels\n",
    "\n",
    "- **What is it?**: Statsmodels is a Python module that provides classes and functions for the estimation of many different statistical models, as well as for conducting statistical tests, and statistical data exploration. The formula API allows for specifying models using R-style formulas.\n",
    "- **When to use it?**: Choose statsmodels when you need more detailed statistical information about the model, such as p-values, confidence intervals, and diagnostic tests.\n",
    "### SciKit Learn\n",
    "\n",
    "- **What is it?**: Scikit-learn is a free machine learning library for Python. It features various classification, regression, and clustering algorithms. It is designed to interoperate with the Python numerical and scientific libraries NumPy and SciPy. \n",
    "- **When to use it?**: Use scikit-learn when you might be interested in exploring more advanced machine learning models beyond linear regression. \n",
    "\n",
    "## Fitting Linear Models to Data\n",
    "\n",
    "In the first section, we consider \"linear models\" to be models where the relationship between the independent and dependent variables is linear. This means that the dependent variable can be expressed as a linear combination of the independent variables.\n",
    "This might encompass linear regression, but can also apply to polynomial regression, where the relationship between the variables is still linear, even if the model itself is polynomial. \n",
    "Our example data happens to be one that can be fit with linear regression. However, the many concepts in this lesson can be applied to other functional forms.\n",
    "\n",
    "A common exercise in introductory chemistry is fitting a linear model to absorbance measurements at a particular wavelength vs. concentration of an analyte using spectrophotometry. \n",
    "In these cases, the absorbance of a sample is described using [Beer-Lambert Law](https://en.wikipedia.org/wiki/Beer%E2%80%93Lambert_law)\n",
    "\n",
    "$$\n",
    "A = \\varepsilon \\cdot l \\cdot c\n",
    "$$\n",
    "\n",
    "where $A$ is the absorbance, $\\varepsilon$ is the absorptivity constant, $l$ is path length, and $c$ is the concentration of the sample.\n",
    "\n",
    "In the first section of this notebook, we will read in absorbance data, fit a linear model, and calculate the concentrations of unknowns using our model using many different Python libraries.\n",
    "\n",
    "The cell below imports the libraries we will use for our analysis. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd # for reading data from a file and putting in a table\n",
    "import matplotlib.pyplot as plt # for plotting\n",
    "\n",
    "import seaborn as sns # for making plots look nice\n",
    "# apply the seaborn theme to plots\n",
    "sns.set_theme()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this lesson, the data we are working with represents absorbance data recorded in a Bradford Assay for determining protein concentration.\n",
    "\n",
    "After we have imported our libraries, we will next use `pandas` to read in our data.\n",
    "Our data is stored in a comma separated value (CSV) file, though pandas can also read from Excel files.\n",
    "\n",
    "In the example in this notebook, we are reading the CSV directly from a web URL instead of from a local file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in file - reading a CSV from a web url.\n",
    "df = pd.read_csv(\"https://raw.githubusercontent.com/MolSSI-Education/python-experimental-chem/main/data/protein_assay.csv\")\n",
    "\n",
    "# View the first five rows.\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing Data\n",
    "\n",
    "After reading in our data, we might wish to inspect it visually. \n",
    "We can do this using `matplotlib`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the data\n",
    "plt.scatter(df['concentration'], df['absorbance'])\n",
    "plt.xlabel('Concentration')\n",
    "plt.ylabel('Absorbance')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting using NumPy Polynomials\n",
    "\n",
    "The [polyfit](https://numpy.org/doc/stable/reference/generated/numpy.polyfit.html) function in NumPy can be used to fit polynomial models. \n",
    "When fitting a polynomial using NumPy, you specify the degree of the polynomial you would like to fit as well as your x and y data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_data = df[\"concentration\"]\n",
    "y_data = df[\"absorbance\"]\n",
    "\n",
    "my_fit = np.polyfit(x_data, y_data, deg=1, full=True)\n",
    "coeff = my_fit[0]\n",
    "fit_data = my_fit[1:]\n",
    "\n",
    "print(f\"The fit slope is {coeff[0]} and the y-intercept is {coeff[1]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Measures like the $R^2$ value, which is a measure of how well the model fits the data, are not directly available when using `polyfit`.\n",
    "However, you can calculate the $R^2$ value by hand using the residuals, though we will not show it here.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NumPy allows you to create a polynomial function from your fit that you can use to predict new values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a polynomial function from fit values\n",
    "poly_fit = np.poly1d(coeff)\n",
    "\n",
    "predicted = poly_fit(x_data)\n",
    "print(predicted)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cell below shows a plot of the data and the polynomial fit using NumPy's `polyfit` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_label = f\"{coeff[0]:.3f}x + {coeff[1]:.3f}\"\n",
    "\n",
    "plt.scatter(x_data, y_data)\n",
    "plt.plot(x_data, predicted, color=\"black\", linestyle=\"--\", label=fit_label)\n",
    "plt.xlabel(\"Concentration\")\n",
    "plt.ylabel(\"Absorbance\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that an advantage of using `polyfit` is that we could also easily try fitting a higher-order polynomial to the data if we wanted to.\n",
    "This can be an advantage if you are not sure what the best model is for your data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting a Linear Model using `scipy.linregress`\n",
    "\n",
    "The [linregress](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.linregress.html) function from SciPy can be used to perform linear regression.\n",
    "This method is best if you are sure that your data is linear and you need simple statistics such as the slope, intercept, and correlation coefficient.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "\n",
    "lin_fit = stats.linregress(x_data, y_data)\n",
    "\n",
    "fit_label = f\"{lin_fit.slope:.3f}x + {lin_fit.intercept:.3f}\\n\"\n",
    "\n",
    "print(f\"Fit equation: {fit_label}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SciPy `linregress` does not provide a built-in way to calculate predictions from the model, but you can easily calculate these by hand using the slope and intercept values returned by the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get predicted y values\n",
    "predicted = lin_fit.slope * x_data + lin_fit.intercept\n",
    "print(predicted)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_label += \"\\n R-squared: \" + f\"{lin_fit.rvalue**2:.3f}\"\n",
    "plt.scatter(x_data, y_data)\n",
    "plt.plot(x_data, predicted, color=\"black\", linestyle=\"--\", label=fit_label)\n",
    "plt.xlabel(\"Concentration\")\n",
    "plt.ylabel(\"Absorbance\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting a Linear Model using `statsmodels`\n",
    "\n",
    "To fit a linear equation to our data, we can also use a library called [statsmodels](https://www.statsmodels.org/stable/index.html). \n",
    "`statsmodels` is a Python module that provides classes and functions for the estimation of many different statistical models, as well as for conducting statistical tests, and statistical data exploration.\n",
    "\n",
    "In particular, we are using a part of `statsmodels` called the formula API.\n",
    "The formula API lets you define a formula as a string for fitting and is specifically designed to work with dataframes.\n",
    "When defining a formula, you use the column names in a string to define the relationship.\n",
    "Note that the column names should not have spaces, or entering the relationship is a bit more complicated.\n",
    "\n",
    "For example, if you have a dataframe `df` with columns `x` and `y`, you can define a formula for a linear model as `y ~ x`. \n",
    "As a slightly more complicatd example one could also fit something like a polynomial using `\"y ~ np.power(x, 2) + x\"` if you had imported NumPy (`import numpy as np`). \n",
    "\n",
    "To use the formula API, we will use `smf` (imported in first cell). \n",
    "We will use ordinary least squares (`ols`) for our fit, though [a number of other options are offered](https://www.statsmodels.org/dev/api.html#statsmodels-formula-api).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.formula.api as smf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "regression = smf.ols(\"absorbance ~ concentration\", data=df).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see a summary of the fit including the `R-squared` by using the `.summary()` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "regression.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "regression.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see what our model predicts for our input concentration values by using the `predict` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "statsmodel_predicted = regression.predict()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_label = f\"{regression.params.Intercept:.3f} + {regression.params.concentration:.3f}x\\n\"\n",
    "fit_label += f\"R-squared: {regression.rsquared:.3f}\"\n",
    "\n",
    "plt.scatter(x_data, y_data)\n",
    "plt.plot(x_data, predicted, color=\"black\", linestyle=\"--\", label=fit_label)\n",
    "plt.xlabel(\"Concentration\")\n",
    "plt.ylabel(\"Absorbance\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting a Linear Model using `scikit-learn`\n",
    "\n",
    "The final library we will use to fit a linear model is `scikit-learn`.\n",
    "`scikit-learn` is a free machine learning library for Python that features various classification, regression, and clustering algorithms. Conveniently, SciKit-Learn also uses a particular format for its models that allows you to switch out models.\n",
    "\n",
    "In this notebook, we will see how to use scikitlearn to fit chemical data to linear models. However, the concepts outlined here will be relevant to fitting data to other models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression # import the linear regression model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "linear_model = LinearRegression()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now need to get our data in a format that is friendly for Scikit-learn. \n",
    "So far, we have been working with data that are single pandas columns (also called a Series), but scikit-learn requires a NumPy array or a pandas DataFrame.\n",
    "One way we can get our data into the right format is to use double brackets when selecting our columns from the DataFrame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_data = df[[\"concentration\"]] # use double brackets to make it a 2D array\n",
    "y_data = df[['absorbance']] # use double brackets to make it a 2D array\n",
    "\n",
    "linear_model.fit(x_data, y_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After `fit` is called, the model is trained on the data. Usually, when using SciKit Learn models, you wouldn't really \"print\" out the model coefficients, but we are doing so here for demonstration purposes.\n",
    "Instead, you would use the trained model to predict new values based on the model fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print model fit\n",
    "print(f\"Scikit Learn Fit: {linear_model.coef_[0][0]:.3f}x + {linear_model.intercept_[0]:.3f}\")\n",
    "\n",
    "# print the model score\n",
    "print(f\"R-squared: {linear_model.score(x_data, y_data):.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Usually a trained SciKit Learn model is used for prediction (as shown below), instead of for printing out the coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_predictions = linear_model.predict(x_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_label = f\"{linear_model.coef_[0][0]:.3f}x + {linear_model.intercept_[0]:.3f}\\n\"\n",
    "fit_label += f\"R-squared: {linear_model.score(x_data, y_data):.3f}\"\n",
    "\n",
    "plt.scatter(x_data, y_data)\n",
    "plt.plot(x_data, model_predictions, color=\"black\", linestyle=\"--\", label=fit_label)\n",
    "plt.xlabel(\"Concentration\")\n",
    "plt.ylabel(\"Absorbance\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting Nonlinear Models with SciKit Learn\n",
    "\n",
    "All scikit learn models use the same API, or interface. This means to switch from a linear model to a more sophisticated model like a random forest model, one need only change the model creation.\n",
    "\n",
    "For example, recall our code to fit a linear model and use it for prediction:\n",
    "\n",
    "```python\n",
    "from sklearn.linear_model import LinearRegression \n",
    "\n",
    "model = LinearRegression()\n",
    "model.fit(X,Y)\n",
    "predictions = model.predict(X)\n",
    "To do the same thing with a random forest regresso, the code would be:\n",
    "\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "\n",
    "model = RandomForestRegressor()\n",
    "model.fit(X, Y)\n",
    "predictions = model.predict()\n",
    "```\n",
    "\n",
    "This makes SciKit Learn especially appropriate if you are interested in exploring more advanced machine learning models beyond linear regression."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting Arbitrary Functional Forms to Data\n",
    "\n",
    "In the previous sections, we fit linear models to data including linear and polynomial fits.\n",
    "However, there are times when it will be beneficial to fit data to arbitary functional forms that may not represent linear relationships between the independent and dependent variables.\n",
    "\n",
    "One example of this type of function is Gaussian distribution. \n",
    "To see how to use SciPy `curve_fit` to fit arbitrary functional forms (example Gaussian) to data, see this [Data Fitting lesson](https://education.molssi.org/python-data-analysis/03-data-fitting/index.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}