> For the complete documentation index, see [llms.txt](https://zavalab.gitbook.io/statschemebook/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://zavalab.gitbook.io/statschemebook/chapter-6-statistical-data-analysis-and-learning/01-classification-models.md).

# Part III: Machine Learning Models

{% hint style="success" %}
**Slides:** [stats\_ch6\_partIII.pdf](https://github.com/vzavala/StatsBook/blob/master/Slides/stats_ch6_partIII.pdf) · [4-up handout](https://github.com/vzavala/StatsBook/blob/master/Slides/handouts/stats_ch6_partIII_4up.pdf) · **Code:** [Python notebooks](https://github.com/vzavala/StatsBook/tree/master/Python) · [MATLAB scripts](https://github.com/vzavala/StatsBook/tree/master/Matlab)
{% endhint %}

* Statistical learning provides powerful tools to construct predictive models from data.
* Models explored here are generalizations of models that we have analyzed.
* We begin discussion with *classification* models.
* Have inputs $$X=(X\_1,...,X\_n)$$ with obs $$x\_\omega$$ that lie in continuous domain $$\mathcal{D}\_X \subseteq \mathbb{R}^n$$.
* Have output $$Y$$ with scalar obs $$y\_{\omega}$$ that lie in (binary) domain $$\mathcal{D}\_Y={0,1}$$.
* We postulate relationship between $$X$$ and $$Y$$:

  $$
  Y=m(X,\theta)+\epsilon
  $$

  where $$\theta \in \mathbb{R}^{n+1}$$ are params, $$m(X,\theta)$$ is model, and $$\epsilon$$ is a noise RV.
* Goal is to find parameters that best explain our observations (data):

  $$
  y\_\omega=m(x\_\omega,\theta)+\epsilon\_\omega,; \omega \in \mathcal{S}
  $$
* Estimation problem is known as a *classification* problem.
* Inputs $$X$$ are known as features or descriptors.
* Output $$Y$$ is known as label or class.
* What is different about class models is binary nature of output $$Y$$.
* Class models have many apps and binary nature makes them intuitive to analyze.
* Think about, for example, a chemical classification problem: given features of chemical (e.g., its molecular structure), can we predict if this is toxic or not?
* Similarly, imagine you want to detect illness from blood pressure and glucose levels.
* To solve class problem, we use technique called logistic regression.
* Postulate model of form:

  $$
  m(x,\theta)=\frac{1}{1+e^{-\theta^Tx}}
  $$
* Mixture $$\theta^Tx=\theta\_0+\sum\_{i=1}^n\theta\_ix\_i$$ is known as *evidence*.
* Params attribute importance (weights) to different features.
* Postulated model is logistic function and captures binary 0-1 logic.
* Logistic function has this behavior:
  * $$m(x,\theta)\to 1;$$ as $$;\theta^Tx\to +\infty;$$ (large evidence)
  * $$m(x,\theta)\to 0;$$ as $$;\theta^Tx\to -\infty;$$ (small evidence)
  * $$m(x,\theta)=\frac{1}{2};$$ if $$;\theta^Tx=0;$$ (no evidence)
* Reversing sign of param has effect of reversing logic.
* Param $$\theta\_0$$ does not weight any input; known as bias (ignores evidence).

<figure><img src="https://1753452229-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FakKE106umUCSWFtTnl31%2Fuploads%2Fgit-blob-204c466bebeb73f44019b7b4a71decd1412b18fd%2Fch6_logistic_function.png?alt=media" alt="" width="630"><figcaption><p>Behavior of logistic function for θ = −1, θ = 0, and θ = +1.</p></figcaption></figure>

#### Python implementation

{% hint style="success" %}
Notebook [`ch6_logistic_function.ipynb`](https://colab.research.google.com/github/vzavala/StatsBook/blob/master/Python/ch6_logistic_function.ipynb) — open it in Colab and run the cells alongside this section.
{% endhint %}

The logistic function $$m(x,\theta)=1/(1+e^{-\theta^Tx})$$ for a scalar input and $$\theta=-1,0,+1$$ (written for the book to accompany the slide figure). The evidence is $$\theta x$$.

**Note:** This notebook reproduces the figure of the logistic function in Section 6.4 of *Statistics for Chemical Engineers*.

```python
import numpy as np
import matplotlib.pyplot as plt

def m(x, theta):
    return 1 / (1 + np.exp(-theta * x))

x = np.linspace(-10, 10, 400)
fig, axs = plt.subplots(1, 3, figsize=(13, 3.5))
for ax, theta in zip(axs, (-1, 0, 1)):
    ax.plot(x, m(x, theta), 'k', linewidth=1.5)
    ax.set_title(rf'$\theta={theta}$')
    ax.set_xlabel('$x$'); ax.set_ylabel(r'$m(x,\theta)$')
    ax.set_ylim(-0.05, 1.05)
    ax.grid(True)
plt.tight_layout()
plt.show()
```

<figure><img src="https://1753452229-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FakKE106umUCSWFtTnl31%2Fuploads%2Fgit-blob-c055fe5785d0c3762dec7851b721c933129dfec5%2Fch6_logistic_function-03-0.png?alt=media" alt="" width="520"><figcaption></figcaption></figure>

**Note:** The logistic function maps any input to a probability between 0 and 1; its parameters shift and sharpen the transition.

* Function $$m(x,\theta)$$ is used to model conditional probs: $$\mathbb{P}(Y=1,|,x,\theta)$$ and $$\mathbb{P}(Y=0,|,x,\theta)=1-\mathbb{P}(Y=1,|,x,\theta)$$.
* Specifically, want to capture logic:
  * If evidence $$\theta^Tx$$ is large, there is a high prob that $$Y=1$$
  * If evidence $$\theta^Tx$$ is small, there is a high prob that $$Y=0$$
  * If no evidence ($$\theta^Tx=0$$), there is ambiguity (equally probable that $$Y=0$$ or $$Y=1$$)
* This is captured by defining conditional pdf:

  $$
  f(y,,|x,\theta)=\mathbb{P}(Y=y|x,\theta)=m(x,\theta)^y(1-m(x,\theta))^{1-y}.
  $$
* We apply MLE to estimate the parameters $$\theta$$.
* Goal is to find $$\hat{\theta}$$ that maximizes the joint probability $$\prod\_{\omega \in \mathcal{S}}f(y\_\omega,|,x\_\omega,\theta)$$:

  $$
  \max\_{\theta}; \sum\_{\omega \in \mathcal{S}}y\_\omega\cdot \log(m(x\_\omega,\theta))+(1-y\_\omega)\cdot \log(1-m(x\_\omega,\theta))
  $$

## Example: Classification for Fault Detection

* Imagine we have data $$(P\_\omega,C\_\omega,T\_\omega)$$ for Gibbs reactor.
* Data is separated into two classes; one class represents normal mode (when $$T\_\omega$$ is low) and we define corresponding labels as $$y\_\omega=1$$.
* The other class represents failure mode (when $$T\_\omega$$ is high) and we define $$y\_\omega=0$$.
* Our goal is to develop a classification model that predicts mode of operation by using data on pressure and conversion $$x\_\omega=(P\_\omega,C\_\omega)$$.
* In other words, model should be able to detect if reactor is failing by looking at a data pair of pressure and concentration. Is that possible?

<figure><img src="https://1753452229-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FakKE106umUCSWFtTnl31%2Fuploads%2Fgit-blob-65c0c9307401be152eea5b11e7a95013d7048d9b%2Fch6_gibbs_logistic_data.png?alt=media" alt="" width="420"><figcaption><p>Gibbs reactor data showing separation between normal and failure operating modes.</p></figcaption></figure>

* Model correctly predicts the class 85% of time (214 out of 250 obs).
* Logistic class also gives estimate of prob that predictions are in a given class.

  <figure><img src="https://1753452229-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FakKE106umUCSWFtTnl31%2Fuploads%2Fgit-blob-fba40423ebaddcd4e2f0fdf68e875bb90621ce5a%2Fch6_gibbs_logistic_confusion.png?alt=media" alt="" width="489"><figcaption><p>Confusion matrix (left) and classification probabilities (right) for Gibbs reactor data.</p></figcaption></figure>

#### Python implementation

{% hint style="success" %}
Notebook [`ch6_gibbs_logistic.ipynb`](https://colab.research.google.com/github/vzavala/StatsBook/blob/master/Python/ch6_gibbs_logistic.ipynb) — open it in Colab and run the cells alongside this section.
{% endhint %}

Logistic classification for fault detection in the Gibbs reactor (Victor Z and Ugochukwu Ikegwu, UW-Madison, 2024).

| Slides                                       | Python                 | Meaning                                                 |
| -------------------------------------------- | ---------------------- | ------------------------------------------------------- |
| $$x\_\omega=(P\_\omega,C\_\omega)$$          | `X` (columns `P`, `C`) | features: pressure and conversion                       |
| $$y\_\omega\in{0,1}$$                        | `y`                    | label: 1 = normal (low $$T$$), 0 = failure (high $$T$$) |
| $$\theta$$, $$\theta^Tx$$                    | `theta`, `evidence`    | parameters and evidence                                 |
| $$m(x,\theta)=\mathbb{P}(Y=1\mid x,\theta)$$ | `m_hat`                | classification probability                              |
| $$\hat{y}\_\omega$$                          | `y_hat`                | predicted label                                         |

**Note:** This notebook reproduces Example 6.22 (Section 6.4) of *Statistics for Chemical Engineers*: logistic regression to detect faults of the Gibbs reactor.

```python
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import confusion_matrix
from sklearn.model_selection import train_test_split

data = np.vstack((np.loadtxt(path+'Data/gibbs_hightemp_class.dat'),
                  np.loadtxt(path+'Data/gibbs_lowtemp_class.dat')))
S = len(data)
y = (data[:, 0] == 583.15).astype(int)          # 1 = normal (low T), 0 = failure
np.random.seed(1)
P = data[:, 1] + np.random.randn(S) * 0.1
C = data[:, 2] + np.random.randn(S) * 0.1
X = np.column_stack((P, C))

plt.scatter(P[y == 0], C[y == 0], c='k', s=20, label='$y=0$ (failure)')
plt.scatter(P[y == 1], C[y == 1], c='w', edgecolors='0.5', s=20, label='$y=1$ (normal)')
plt.xlabel('$P$'); plt.ylabel('$C$'); plt.legend(loc='lower right'); plt.grid(True)
plt.show()
```

<figure><img src="https://1753452229-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FakKE106umUCSWFtTnl31%2Fuploads%2Fgit-blob-77ff0f8283a8da72f9929ee1da83fc8ee30e3510%2Fch6_gibbs_logistic-03-0.png?alt=media" alt="" width="520"><figcaption></figcaption></figure>

MLE: $$\max\_\theta \sum\_\omega y\_\omega\log m(x\_\omega,\theta)+(1-y\_\omega)\log(1-m(x\_\omega,\theta))$$. (`penalty=None` gives plain MLE; scikit-learn regularizes by default.)

```python
model = LogisticRegression(penalty=None, max_iter=1000).fit(X, y)
theta = np.concatenate([model.intercept_, model.coef_.ravel()])   # (theta_0, theta_P, theta_C)
evidence = theta[0] + X @ theta[1:]
m_hat = 1 / (1 + np.exp(-evidence))                               # = model.predict_proba(X)[:, 1]
y_hat = (m_hat >= 0.5).astype(int)
print("theta_hat =", np.round(theta, 3))
print(f"accuracy = {np.mean(y_hat == y):.2%}  ({np.sum(y_hat == y)} of {S})")

def plot_confusion(ax, cm, title):
    ax.imshow(cm, cmap='Greys', alpha=0.3)
    for (i, j), v in np.ndenumerate(cm):
        ax.text(j, i, v, ha='center', va='center', fontsize=14)
    ax.set_xticks([0, 1]); ax.set_yticks([0, 1])
    ax.set_xlabel('predicted'); ax.set_ylabel('true'); ax.set_title(title)

fig, axs = plt.subplots(1, 2, figsize=(11, 4))
plot_confusion(axs[0], confusion_matrix(y, y_hat), 'confusion matrix')
axs[1].plot(m_hat, 'ko', markerfacecolor='w', markersize=3)
axs[1].set_xlabel(r'$\omega$'); axs[1].set_ylabel(r'$\mathbb{P}(y_\omega=1\mid x_\omega,\hat{\theta})$')
axs[1].grid(True)
plt.tight_layout()
plt.show()
```

```
theta_hat = [-2.229 -0.071 21.873]
accuracy = 87.80%  (439 of 500)
```

<figure><img src="https://1753452229-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FakKE106umUCSWFtTnl31%2Fuploads%2Fgit-blob-c58a726ff416fc71375d61ec76e4c52d79800061%2Fch6_gibbs_logistic-05-1.png?alt=media" alt="" width="520"><figcaption></figcaption></figure>

Train/test split (60/40) to check that the classifier generalizes.

```python
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.4, random_state=42)
model_tr = LogisticRegression(penalty=None, max_iter=1000).fit(X_tr, y_tr)
fig, axs = plt.subplots(1, 2, figsize=(9, 4))
plot_confusion(axs[0], confusion_matrix(y_tr, model_tr.predict(X_tr)), 'training')
plot_confusion(axs[1], confusion_matrix(y_te, model_tr.predict(X_te)), 'test')
plt.tight_layout()
plt.show()
```

<figure><img src="https://1753452229-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FakKE106umUCSWFtTnl31%2Fuploads%2Fgit-blob-ca306ec7818bbcd844aaed98a6d1911a03b2af74%2Fch6_gibbs_logistic-07-0.png?alt=media" alt="" width="520"><figcaption></figcaption></figure>

**Note:** As in the book, the model separates the classes well; cross-validation should be used to check that it generalizes to new data.