09-LinearRegression.Rmd

# Linear Regression



In linear regression, the goal is to find a function \( f \) that maps inputs \( \mathbf{x} \in \mathbb{R}^D \) to outputs \( f(\mathbf{x}) \in \mathbb{R} \).  We are given training data \( \{(\mathbf{\mathbf{x}}_n, y_n)\}_{n=1}^N \), where each observation is modeled as:
\[
y_n = f(\mathbf{x}_n) + \epsilon_n.
\]
Here, \( \epsilon_n \) represents independent and identically distributed (i.i.d.) Gaussian noise with zero mean.

With linear regression, we have two basic goals:

- **Modeling the data:** Estimate the function \( f \) that explains the observed data.  
- **Generalization:** Ensure \( f \) predicts well for unseen inputs, not just the training data.


When using linear regression, we have to make some decisions apriori.

1. **Model Choice:**  
   - Decide on the *type* and *parametrization* of the regression function (e.g., polynomial degree).  
   - Model selection (see Section 8.6) helps identify the simplest model that explains the data effectively.

2. **Parameter Estimation:**  
   - Determine the optimal model parameters using appropriate *loss functions* and *optimization algorithms*.

3. **Overfitting and Model Selection:**  
   - Overfitting occurs when the model fits training data too closely, failing to generalize.  
   - This typically arises from overly flexible or complex models.

4. **Connection Between Loss Functions and Priors:**  
   - Many optimization objectives can be derived from probabilistic assumptions about the data.  
   - Understanding this relationship clarifies how prior beliefs influence model behavior.

5. **Uncertainty Modeling:**  
   - Since training data are finite, predictions carry uncertainty.  
   - Modeling uncertainty provides *confidence bounds* for predictions, which are especially important with limited data.

---

## Problem Formulation

In regression, we model the relationship between inputs \( \mathbf{x} \in \mathbb{R}^D \) and outputs \( y \in \mathbb{R} \) in the presence of observation noise.  We assume that each observation follows a probabilistic model:
\[
p(y | \mathbf{x}) = \mathcal{N}(y \mid f(\mathbf{\mathbf{x}}), \sigma^2).
\]
This means the data are generated by:
\[
y = f(\mathbf{x}) + \epsilon, \quad \epsilon \sim \mathcal{N}(0, \sigma^2),
\]
where \( \epsilon \) is i.i.d. Gaussian noise with zero mean and variance \( \sigma^2 \).



The goal is to find a function \( f(\mathbf{x}) \) that:

- Appro\mathbf{x}imates the true (unknown) data-generating function,  
- Generalizes well to unseen data.

---

### Parametric Models

We restrict our attention to parametric models, where the function depends on a set of parameters \( \theta \). 
In linear regression the model is *linear in the parameters*:
\[
p(y | \mathbf{x}, \theta) = \mathcal{N}(y \mid \mathbf{x}^T \theta, \sigma^2),
\]
or equivalently,
\[
y = \mathbf{x}^T \theta + \epsilon, \quad \epsilon \sim \mathcal{N}(0, \sigma^2)
\]
Here:

- \( \theta \in \mathbb{R}^D \) represents the model parameters,  
- The likelihood e\mathbf{x}presses how probable a given \( y \) is for known \( \mathbf{x} \) and \( \theta \),  
- Without noise (\( \sigma^2 \to 0 \)), the model becomes deterministic (a Dirac delta).  




<div class="example">
Linear Regression Model

For \( \mathbf{x}, \theta \in \mathbb{R} \):

- The model \( y = \theta \mathbf{x} \) describes a **straight line** through the origin.  
- The slope of the line is given by \( \theta \).  
- Different values of \( \theta \) yield different linear functions.  
</div>

Although this model is linear in both \( \mathbf{x} \) and \( \theta \), we can generalize it by introducing **nonlinear feature transformations**:
\[
y = \mathbf{\phi}(\mathbf{x})^T \theta.
\]
Here, \( \mathbf{\phi}(\mathbf{x}) \) represents a **feature mapping** of the input \( \mathbf{x} \).  
The model remains linear **in the parameters** \( \theta \), even if \( \mathbf{\phi}(\mathbf{x}) \) is nonlinear in \( \mathbf{x} \).


For now, we assume that the noise variance \( \sigma^2 \) is known and focus on learning the optimal parameters \( \theta \).

---



### Exercises {.unnumbered .unlisted}



<div class="exercise">
Finding a regression function requires solving a variety of problems.  List and discuss them. 

        
<div style="text-align: right;">
[Solution]( )
</div> 
</div>


<div class="exercise">
Consider the following data: \[(1,3),  (2,4),  (3,8),  (4,9).\]  Find the estimated regression line \[y = \beta_0 + \beta_1 x,\] based on the data.   For each $x_i$, compute both the estimated value of $y$ and the residuals. 

<div style="text-align: right;">
[Solution]( )
</div> 
</div>



<div class="exercise">
A simple linear regression model is fit, relating plant growth over 1 year (y) to amount of fertilizer provided (x). Twenty five plants are selected, 5 each assigned to each of the fertilizer levels (12, 15, 18, 21, 24).  The results of the model fit are given below:

**Regression Coefficients**

| Model    |   B    | Std. Error |   t    | Sig |
|----------|--------|------------|--------|-----|
| Constant | 8.624  | 1.81       | 4.764  | 0   |
| $x$      | 0.527  | 0.098      | 5.386  | 0   |

Can we conclude that there is an association between fertilizer and plant growth at the 0.05 significance level?
        
<div style="text-align: right;">
[Solution]( )
</div> 
</div>




<div class="exercise">

A multiple regression model is fit, relating salary (Y) to the following predictor variables: experience ($X_1$, in years), accounts in charge of ($X_2$) and gender ($X_3$ is 1 if female, 0 if male). The following ANOVA table and output gives the results for fitting the model. Conduct all tests at the 0.05 significance level:

**ANOVA Table**

|            | df |   SS   |   MS   |   F   | P-value |
|------------|----|--------|--------|-------|---------|
| Regression |  3 | 2470.4 | 823.5  | 76.9  | 0       |
| Residual   | 21 | 224.7  | 10.7   |       |         |
| Total      | 24 | 2695.1 |        |       |         |

**Regression Coefficients**

| Model       |   B    | Std. Error |   t    |  Sig   |
|-------------|--------|------------|--------|--------|
| Constant    | 39.58  | 1.89       | 21.00  | 0      |
| Experience  | 3.61   | 0.36       | 10.04  | 0      |
| Accounts    | -0.28  | 0.36       | -0.79  | 0.4389 |
| Gender      | -3.92  | 1.48       | -2.65  | 0.0149 |

Test whether salary is associated with any of the predictor variables.
        
<div style="text-align: right;">
[Solution]( )
</div> 
</div>







##  Parameter Estimation

In linear regression, we are given a training dataset  
\[
\mathcal{D} = \{(\mathbf{x}_1, y_1), \ldots, (\mathbf{x}_N, y_N)\},
\]
where \( \mathbf{x}_n \in \mathbb{R}^D \) are inputs and \( y_n \in \mathbb{R} \) are corresponding observations.  Let $\mathcal{X}$ and $\mathcal{Y}$ be the sets of training inputs and traning targets respectively. We assume each target is generated according to a probabilistic model:
\[
p(y_n | \mathbf{x}_n, \theta) = \mathcal{N}(y_n| \mathbf{x}_n^\top \theta, \sigma^2),
\]
where \( \theta \) are model parameters and \( \sigma^2 \) is the noise variance.  Because each data point is conditionally independent given the parameters, the likelihood factorizes as:
\[
p(\mathcal{Y} | \mathcal{X}, \theta) = \prod_{n=1}^N p(y_n | \mathbf{x}_n, \theta) = \prod_{n=1}^N \mathcal{N}(y_n| \mathbf{x}_n^\top \theta, \sigma^2).
\]

---

### Maximum Likelihood Estimation (MLE)

The **Maximum Likelihood Estimation** approach finds parameters that maximize the likelihood:
\[
\theta_{\text{ML}} = \arg \max_\theta p(\mathcal{Y} | \mathcal{X}, \theta).
\]

Since taking the logarithm does not change the location of the optimum, we minimize the **negative log-likelihood** instead:
\[
L(\theta) = -\log p(\mathcal{Y} | \mathcal{X}, \theta) = \frac{1}{2\sigma^2} \| y - X\theta \|^2 + const.
\]
This expression is quadratic in \( \theta \), so setting its gradient to zero yields the closed-form solution:
\[
\theta_{\text{ML}} = (X^\top X)^{-1} X^\top y.
\]

Here:

- \( X \in \mathbb{R}^{N \times D} \) is the **design matrix**, where each row is \( x_n^\top \).  
- The matrix \( X^\top X \) must be invertible (i.e., full column rank).

This solution minimizes the sum of squared errors between the predictions \( X\theta \) and the observed targets \( y \).

---

### MLE with Nonlinear Features

Although the term "linear regression" refers to models that are linear in the parameters, the inputs can undergo nonlinear transformations through a feature mapping:
\[
f(x) = \mathbf{\phi}(x)^\top \theta,
\]
where \( \mathbf{\phi}(x) \) is a feature vector (e.g., polynomial terms, basis functions).

We define the **feature matrix** as:
\[
\mathbf{\Phi} =
\begin{bmatrix}
\mathbf{\phi}(x_1)^\top \\
\mathbf{\phi}(x_2)^\top \\
\vdots \\
\mathbf{\phi}(x_N)^\top
\end{bmatrix} \in \mathbb{R}^{N \times K}.
\]

For example, with second-order polynomial features:
\[
\mathbf{\phi}(x) = [1, x, x^2]^\top,
\quad
\mathbf{\Phi} =
\begin{bmatrix}
1 & x_1 & x_1^2 \\
1 & x_2 & x_2^2 \\
\vdots & \vdots & \vdots \\
1 & x_N & x_N^2
\end{bmatrix}.
\]
The corresponding MLE solution becomes:
\[
\theta_{\text{ML}} = (\mathbf{\Phi}^\top \mathbf{\Phi})^{-1} \mathbf{\Phi}^\top y.
\]
The estimated noise variance can also be obtained as:
\[
\sigma^2_{\text{ML}} = \frac{1}{N} \sum_{n=1}^N (y_n - \mathbf{\phi}(x_n)^\top \theta_{\text{ML}})^2.
\]

---

### Overfitting in Linear Regression

Overfitting occurs when a model fits the training data too closely, capturing noise rather than the underlying trend.  To assess performance, we often compute the **Root Mean Square Error (RMSE)**:
\[
\text{RMSE} = \sqrt{\frac{1}{N} \| \mathbf{y} - \mathbf{\mathbf{\Phi}}\mathbf{\theta} \|^2}.
\]

- **Training error** tends to decrease as model complexity (e.g., polynomial degree \( M \)) increases.  
- **Test error** initially decreases but then increases beyond an optimal model complexity, indicating overfitting.

For polynomial regression:

- Low-degree polynomials underfit (poor fit to data).  
- Moderate-degree polynomials (e.g., \( M = 4 \)) generalize well.  
- High-degree polynomials (\( M \geq N - 1 \)) fit all training points but fail to generalize.

---

###  Maximum A Posteriori (MAP) Estimation

Maximum likelihood estimation (MLE) can lead to overfitting, often resulting in very large parameter values.  To mitigate this, we introduce a prior distribution \( p(\theta) \) on the parameters, which encodes our beliefs about plausible parameter values before observing any data.

For example:
\[
p(\theta) = \mathcal{N}(0, 1)
\]
suggests that \( \theta \) is likely to lie within approximately \([-2, 2]\).


Once data \( (\mathcal{X}, \mathcal{Y}) \) are available, we seek parameters that maximize the posterior distribution:
\[
p(\mathbf{\theta} | \mathcal{X}, \mathcal{Y}) = \frac{p(\mathcal{Y} | \mathcal{X}, \mathbf{\theta}) \, p(\mathbf{\theta})}{p(\mathcal{Y} | \mathcal{X})}
\]

The MAP estimate is the maximizer of this posterior:
\[
\mathbf{\theta}_{\text{MAP}} = \arg \max_{\mathbf{\theta}} \, p(\mathbf{\theta} | \mathcal{X}, \mathcal{Y}).
\]
Taking the logarithm, we get:
\[
\log p(\mathbf{\theta} | \mathcal{X}, \mathcal{Y}) = \log p(\mathcal{Y} | \mathcal{X}, \mathbf{\theta}) + \log p(\mathbf{\theta}) + \text{const.}
\]

This shows that the MAP estimate balances:

- The **log-likelihood** term (fit to data), and  
- The **log-prior** term (penalizing unlikely parameters).

Hence, MAP can be viewed as a compromise between fitting the data and respecting prior beliefs.

---

### Optimization
MAP estimation is equivalent to minimizing the **negative log-posterior**:
\[
\mathbf{\theta}_{\text{MAP}} = \arg \min_{\mathbf{\theta}} \{-\log p(\mathcal{Y} | \mathcal{X}, \mathbf{\theta}) - \log p(\mathbf{\theta})\}.
\]


If we assume a Gaussian prior \( p(\mathbf{\theta}) = \mathcal{N}(\mathbf{0}, b^2 \mathbf{I}) \), the negative log-posterior becomes:
\[
- \log p(\mathbf{\theta} | \mathcal{X}, \mathcal{Y}) =
\frac{1}{2\sigma^2}(y - \mathbf{\Phi} \mathbf{\theta})^T (y - \mathbf{\Phi} \mathbf{\theta})
+ \frac{1}{2b^2}\mathbf{\theta}^T \mathbf{\theta} + \text{const.}
\]
Solving for \( \mathbf{\theta} \) yields the MAP estimate:
\[
\mathbf{\theta}_{\text{MAP}} =
\left(\mathbf{\Phi}^T \mathbf{\Phi} + \frac{\sigma^2}{b^2} \mathbf{I} \right)^{-1} \mathbf{\Phi}^T y.
\]

---

### Comparison to MLE
Compared to the MLE solution:
\[
\mathbf{\theta}_{\text{MLE}} = (\mathbf{\Phi}^T \mathbf{\Phi})^{-1} \mathbf{\Phi}^T y,
\]
MAP adds the regularization term \( \frac{\sigma^2}{b^2} \mathbf{I} \), which:

- Ensures invertibility of the matrix,  
- Reduces overfitting,  
- Shrinks parameter values toward zero.  




<div class="example">
In polynomial regression:

- Using a Gaussian prior \( p(\mathbf{\theta}) = \mathcal{N}(0, \mathbf{I}) \),  
- The MAP estimate smooths high-degree polynomials (reducing overfitting),  
- The effect is minimal for low-degree models.  
</div>

However, while MAP helps, it is not a complete solution to overfitting.

---

### MAP Estimation as Regularization


Instead of using a prior, we can add a **regularization term** directly to the loss function:
\[
\min_{\mathbf{\theta}} \|y - \mathbf{\Phi} \mathbf{\theta}\|^2 + \lambda \|\mathbf{\theta}\|_2^2.
\]
Here:

- The first term is the **data-fit (misfit)** term,  
- The second term is the **regularizer**, which penalizes large parameter values,  
- \( \lambda \ge 0 \) controls the strength of regularization.  


The regularization term corresponds to the negative log of a Gaussian prior:
\[
-\log p(\mathbf{\theta}) = \frac{1}{2b^2}\|\mathbf{\theta}\|_2^2 + \text{const.}
\]
If we set \( \lambda = \frac{\sigma^2}{b^2} \), the regularized least squares and MAP estimation solutions are equivalent:
\[
\mathbf{\theta}_{\text{RLS}} = (\mathbf{\Phi}^T \mathbf{\Phi} + \lambda I)^{-1} \mathbf{\Phi}^T y.
\]



<div class="example">
Linear Regression Example (4–6 data points)

We consider a small dataset with 5 observations.  
Let the input be a scalar \( x \) and the response be \( y \):

| \(x\) | 0 | 1 | 2 | 3 | 4 |
|------|---|---|---|---|---|
| \(y\) | 1 | 2 | 2 | 3 | 5 |

We assume **Gaussian noise** throughout.

**Model with Nonlinear Features**

Instead of a purely linear model, we use nonlinear features:
\[
\phi(x) =
\begin{bmatrix}
1 \\
x \\
x^2
\end{bmatrix}
\]

The regression model is:
\[
y = \mathbf{w}^\top \phi(x) + \varepsilon,
\quad \varepsilon \sim \mathcal{N}(0, \sigma^2)
\]

Let:
\[
\mathbf{w} =
\begin{bmatrix}
w_0 \\ w_1 \\ w_2
\end{bmatrix}
\]

**Design Matrix**

The design matrix is:
\[
\mathbf{X} =
\begin{bmatrix}
1 & 0 & 0 \\
1 & 1 & 1 \\
1 & 2 & 4 \\
1 & 3 & 9 \\
1 & 4 & 16
\end{bmatrix},
\quad
\mathbf{y} =
\begin{bmatrix}
1 \\ 2 \\ 2 \\ 3 \\ 5
\end{bmatrix}
\]

**Maximum Likelihood Estimation (MLE)**

The likelihood is:
\[
p(\mathbf{y} \mid \mathbf{X}, \mathbf{w})
=
\mathcal{N}(\mathbf{X}\mathbf{w}, \sigma^2 \mathbf{I})
\]
Maximizing the likelihood is equivalent to minimizing:
\[
\mathcal{L}(\mathbf{w})
=
\|\mathbf{y} - \mathbf{X}\mathbf{w}\|^2
\]

The MLE solution is:
\[
\hat{\mathbf{w}}_{\text{MLE}}
=
(\mathbf{X}^\top \mathbf{X})^{-1}
\mathbf{X}^\top \mathbf{y}
\]

This shows:

- MLE = least squares  
- Nonlinear features handled via \( \mathbf{X} \)  

**Regularization (Ridge Regression)**

To prevent overfitting, add an \( \ell_2 \) penalty:
\[
\mathcal{L}_{\text{reg}}(\mathbf{w})
=
\|\mathbf{y} - \mathbf{X}\mathbf{w}\|^2
+
\lambda \|\mathbf{w}\|^2
\]

The solution becomes:
\[
\hat{\mathbf{w}}_{\text{ridge}}
=
(\mathbf{X}^\top \mathbf{X} + \lambda \mathbf{I})^{-1}
\mathbf{X}^\top \mathbf{y}
\]
This is regularized MLE.

**MAP Estimation**

Assume a Gaussian prior on the parameters:
\[
p(\mathbf{w}) = \mathcal{N}(\mathbf{0}, \tau^2 \mathbf{I})
\]
The posterior is:
\[
p(\mathbf{w} \mid \mathbf{y})
\propto
p(\mathbf{y} \mid \mathbf{w}) p(\mathbf{w})
\]
Maximizing the posterior is equivalent to minimizing:
\[
\|\mathbf{y} - \mathbf{X}\mathbf{w}\|^2
+
\frac{\sigma^2}{\tau^2}\|\mathbf{w}\|^2
\]
Thus the MAP estimator is:
\[
\hat{\mathbf{w}}_{\text{MAP}}
=
(\mathbf{X}^\top \mathbf{X} + \lambda \mathbf{I})^{-1}
\mathbf{X}^\top \mathbf{y}
\quad
\text{where }
\lambda = \frac{\sigma^2}{\tau^2}
\]

</div>



---



### Exercises {.unnumbered .unlisted}




<div class="exercise">
 Why should we maximize the log-likelihood rather than the likelihood?
 
<div style="text-align: right;">
[Solution]( )
</div> 
</div>



<div class="exercise">
 Show $p(Y|X, \mathbf{\theta}) = \prod_{n=1}^N N(y_n|\mathbf{x}_n^T\mathbf{\theta}, \sigma^2)$.
 
<div style="text-align: right;">
[Solution]( )
</div> 
</div>



<div class="exercise">
We know that \[-\log p(Y|X, \mathbf{\theta}) = -\sum_{n=1}^N \log p(y_n|\mathbf{x}_n, \mathbf{\theta}).\] 

<div style="text-align: right;">
[Solution]( )
</div> 
</div>



<div class="exercise">
State the formula for $p(y_n|\mathbf{x}_n, \mathbf{\theta})$.  Assume that $\sigma^2$ is known.  Use it to show \[\log p(y_n|\mathbf{x}_n, \mathbf{\theta}) = -\dfrac{1}{2\sigma^2}(y_n - \mathbf{x}_n^T \mathbf{\theta})^2 + const.\]  What is the constant?  What does it matter that it has no $\theta$ in it? 

<div style="text-align: right;">
[Solution]( )
</div> 
</div>



<div class="exercise">
Suppose $y = 1 + 2x_1 + 3x_2$, $\mathbf{X} = \begin{bmatrix}1 & 4 & -3 \\ 1 & 1 & 1 \end{bmatrix}$, and $\mathbf{y} = \begin{bmatrix}0\\5 \end{bmatrix}$.  Show that
    \[\dfrac{1}{2\sigma^2} \sum_{n=1}^N (y_n - \mathbf{x}_n^T \mathbf{\theta})^2 =  \dfrac{1}{2\sigma^2} ||\mathbf{y} - \mathbf{X}\mathbf{\theta}||^2. \] 
    
<div style="text-align: right;">
[Solution]( )
</div> 
</div>



<div class="exercise">
Define $L(\theta)$ as the negative log-likelihood function.  Show that 
    \begin{align*} 
        L(\mathbf{\theta}) &= \dfrac{1}{2\sigma^2} \sum_{n=1}^N (y_n - \mathbf{x}_n^T \mathbf{\theta})^2 \\ 
        &=  \dfrac{1}{2\sigma^2} ||\mathbf{y} - \mathbf{X}\mathbf{\theta}||^2. \end{align*}
    
<div style="text-align: right;">
[Solution]( )
</div> 
</div>



<div class="exercise">
$\dfrac{dL}{d\mathbf{\theta}} = \dfrac{1}{\sigma^2} \left(-\mathbf{y}\mathbf{X} + \mathbf{\theta}^T \mathbf{X}^T \mathbf{X} \right) \in \bbR^{1 \times x} $.  Justify each step. 
    
<div style="text-align: right;">
[Solution]( )
</div> 
</div>



<div class="exercise">
Solve $\dfrac{dL}{d\mathbf{\theta}} = 0$ to find $\mathbf{\theta}_{ML}$. 
    
<div style="text-align: right;">
[Solution]( )
</div> 
</div>



<div class="exercise">
To find the MAP estimate, begin with \[p(\mathbf{\theta}|X,Y) = \dfrac{p(Y|X, \mathbf{\theta})p(\mathbf{\theta})}{p(Y|X)}.\]  Find the maximum log-likelihood.
    
<div style="text-align: right;">
[Solution]( )
</div> 
</div>









## Bayesian Linear Regression


Bayesian linear regression extends the concept of linear regression by treating parameters as random variables instead of estimating a single best set of parameters (as in MLE or MAP).  Rather than finding a point estimate for parameters \( \mathbf{\theta} \), it computes the full posterior distribution over them and uses this distribution to make predictions.

This approach:

- Incorporates prior beliefs about parameters,  
- Naturally accounts for uncertainty,  
- Prevents overfitting, especially with limited data.  

---

### The Model

We define the model as:
\[
p(\mathbf{\theta}) = \mathcal{N}(\mathbf{m}_0, \mathbf{S}_0) \quad \text{(prior on parameters)}
\]
\[
p(y | \mathbf{x}, \mathbf{\theta}) = \mathcal{N}(y | \phi(\mathbf{x})^\top \mathbf{\theta}, \sigma^2) \quad \text{(likelihood)}.
\]
The joint distribution is:
\[
p(y, \mathbf{\theta} | \mathbf{x}) = p(y | \mathbf{x}, \mathbf{\theta}) \, p(\mathbf{\theta}).
\]
Here:

- \( \mathbf{\theta} \) is now a random variable,  
- \( \mathbf{m}_0 \) and \( \mathbf{S}_0 \) are the prior mean and covariance.


Predictions are obtained by integrating out the parameter uncertainty:
\[
p(y_* | \mathbf{x}_*) = \int p(y_* | \mathbf{x}_*, \mathbf{\theta}) p(\mathbf{\theta}) d\mathbf{\theta}.
\]
Since both likelihood and prior are Gaussian, the predictive distribution is also Gaussian:
\[
p(y_* | \mathbf{x}_*) = \mathcal{N}(\phi(\mathbf{x}_*)^\top \mathbf{m}_0, \phi(\mathbf{x}_*)^\top \mathbf{S}_0 \phi(\mathbf{x}_*) + \sigma^2)
\]

The term \( \phi(\mathbf{x}_*)^\top \mathbf{S}_0 \phi(\mathbf{x}_*) \) reflects uncertainty due to parameter variability. The term \( \sigma^2 \) reflects observation noise.

For noise-free function values \( f(\mathbf{x}_*) = \phi(\mathbf{x}_*)^\top \mathbf{\theta} \):
\[
p(f(\mathbf{x}_*)) = \mathcal{N}(\phi(\mathbf{x}_*)^\top \mathbf{m}_0, \phi(\mathbf{x}_*)^\top \mathbf{S}_0 \phi(\mathbf{x}_*))
\]


The parameter prior \( p(\mathbf{\theta}) \) induces a distribution over functions:

- Each sampled parameter vector \( \mathbf{\theta}_i \sim p(\mathbf{\theta}) \) defines a function \( f_i(\mathbf{x}) = \phi(\mathbf{x})^\top \mathbf{\theta}_i \).  
- The collection of these defines \( p(f(\cdot)) \), a distribution over possible functions.







<div class="example">
Bayesian Linear Regression with a Single Feature


We model a simple linear relationship
\[
y = \theta_0 + \theta_1 x + \varepsilon,
\quad \varepsilon \sim \mathcal{N}(0, \sigma^2).
\]
Define the feature map
\[
\phi(x) =
\begin{bmatrix}
1 \\
x
\end{bmatrix},
\quad
\mathbf{\theta} =
\begin{bmatrix}
\theta_0 \\
\theta_1
\end{bmatrix}.
\]
Then
\[
p(y \mid x, \mathbf{\theta})
= \mathcal{N}(y \mid \phi(x)^\top \mathbf{\theta}, \sigma^2).
\]


Assume a **Gaussian prior**:
\[
p(\mathbf{\theta})
= \mathcal{N}(\mathbf{m}_0, \mathbf{S}_0),
\quad
\mathbf{m}_0 =
\begin{bmatrix}
0 \\ 0
\end{bmatrix},
\quad
\mathbf{S}_0 =
\begin{bmatrix}
1 & 0 \\
0 & 1
\end{bmatrix}.
\]

**Interpretation**:

- Before seeing data, we believe slopes and intercepts are near 0.  
- Large uncertainty allows many possible linear functions.  

Each draw \( \mathbf{\theta}_i \sim p(\mathbf{\theta}) \) defines a function
\[
f_i(x) = \theta_{0,i} + \theta_{1,i} x.
\]
Thus the prior defines a **distribution over lines**.



For a new input \( x_* \),
\[
p(y_* \mid x_*)
= \mathcal{N}
\left(
\phi(x_*)^\top \mathbf{m}_0,\;
\phi(x_*)^\top \mathbf{S}_0 \phi(x_*) + \sigma^2
\right).
\]
Because \( \mathbf{m}_0 = \mathbf{0} \):
\[
\mathbb{E}[y_*] = 0.
\]
The variance
\[
\phi(x_*)^\top \mathbf{S}_0 \phi(x_*)
= 1 + x_*^2
\]
grows with distance from the origin — **uncertainty increases away from data**.

</div>


---

### Posterior Distribution

After observing training data \( (\mathcal{X}, \mathcal{Y}) \), we compute the **posterior** over parameters:

\[
p(\mathbf{\theta} | \mathcal{X}, \mathcal{Y}) = \frac{p(\mathcal{Y} | \mathcal{X}, \mathbf{\theta}) p(\mathbf{\theta})}{p(\mathcal{Y} | \mathcal{X})}
\]

Where:

- \( p(\mathcal{Y} | \mathcal{X}, \mathbf{\theta}) = \mathcal{N}(\mathcal{Y} | \mathbf{\Phi} \mathbf{\theta}, \sigma^2 I) \)  
- \( p(\mathcal{Y} | \mathcal{X}) = \int p(\mathcal{Y} | \mathcal{X}, \mathbf{\theta}) p(\mathbf{\theta}) d\mathbf{\theta} \) is the marginal likelihood (normalizing constant).


Because both prior and likelihood are Gaussian, the posterior is also Gaussian:
\[
p(\mathbf{\theta} | X, Y) = \mathcal{N}(\mathbf{\theta} | \mathbf{m}_N, \mathbf{S}_N),
\]
with
\[
\mathbf{S}_N = (\mathbf{S}_0^{-1} + \sigma^{-2} \mathbf{\Phi}^\top \mathbf{\Phi})^{-1}, \quad
\mathbf{m}_N = \mathbf{S}_N (\mathbf{S}_0^{-1} \mathbf{m}_0 + \sigma^{-2} \mathbf{\Phi}^\top Y).
\]
This is derived by completing the square in the exponent of the unnormalized posterior.




<div class="example">
Posterior After Observing Data


Suppose we observe:
\[
\mathcal{X} =
\begin{bmatrix}
0 \\
1 \\
2
\end{bmatrix},
\quad
\mathcal{Y} =
\begin{bmatrix}
1 \\
2 \\
2
\end{bmatrix}.
\]
The design matrix is
\[
\mathbf{\Phi} =
\begin{bmatrix}
1 & 0 \\
1 & 1 \\
1 & 2
\end{bmatrix}.
\]



Because the prior and likelihood are Gaussian, the posterior is:
\[
p(\mathbf{\theta} \mid \mathcal{X}, \mathcal{Y})
= \mathcal{N}(\mathbf{m}_N, \mathbf{S}_N),
\]
with
\[
\mathbf{S}_N
= (\mathbf{S}_0^{-1} + \sigma^{-2} \mathbf{\Phi}^\top \mathbf{\Phi})^{-1},
\quad
\mathbf{m}_N
= \mathbf{S}_N
(\mathbf{S}_0^{-1} \mathbf{m}_0 + \sigma^{-2} \mathbf{\Phi}^\top \mathcal{Y}).
\]

**Interpretation**:

- Posterior mean \( \mathbf{m}_N \): best estimate of parameters (MAP).  
- Posterior covariance \( \mathbf{S}_N \): remaining uncertainty after seeing data.  

Uncertainty shrinks in directions well-supported by data.

</div>




---

### Posterior Predictions

To predict at a new point \( \mathbf{x}_* \), we again integrate over \( \mathbf{\theta} \), but now using the posterior instead of the prior:
\[
p(y_* | \mathcal{X}, \mathcal{Y}, \mathbf{x}_*) = \int p(y_* | \mathbf{x}_*, \mathbf{\theta}) p(\mathbf{\theta} | \mathcal{X}, \mathcal{Y}) d\mathbf{\theta}.
\]
This yields:
\[
p(y_* | \mathcal{X}, \mathcal{Y}, \mathbf{x}_*) = \mathcal{N}(\phi(\mathbf{x}_*)^\top \mathbf{m}_N, \phi(\mathbf{x}_*)^\top \mathbf{S}_N \phi(\mathbf{x}_*) + \sigma^2).
\]

- The predictive mean \( \phi(\mathbf{x}_*)^\top \mathbf{m}_N \) equals the MAP estimate prediction.  
- The predictive variance accounts for both model and observation uncertainty.


For \( f(\mathbf{x}_*) = \phi(\mathbf{x}_*)^\top \mathbf{\theta} \):
\[
\mathbb{E}[f(\mathbf{x}_*) | \mathcal{X}, \mathcal{Y}] = \phi(\mathbf{x}_*)^\top \mathbf{m}_N, \quad
\text{Var}[f(\mathbf{x}_*) | \mathcal{X}, \mathcal{Y}] = \phi(\mathbf{x}_*)^\top \mathbf{S}_N \phi(\mathbf{x}_*)
\]


Sampling \( \mathbf{\theta}_i \sim p(\mathbf{\theta} | \mathcal{X}, \mathcal{Y}) \) gives functions \( f_i(\mathbf{x}) = \phi(\mathbf{x})^\top \mathbf{\theta}_i \).  These represent the posterior distribution over functions, with:

- Mean function: \( \mathbf{m}_N^\top \phi(\mathbf{x}) \),  
- Variance: \( \phi(\mathbf{x})^\top \mathbf{S}_N \phi(\mathbf{x}) \).

**Interpretation:**  
Higher uncertainty (larger variance) reflects regions with sparse or no data.







<div class="example">
Posterior Predictive Distribution

For a new input \( x_* = 1.5 \):

\[
p(y_* \mid \mathcal{X}, \mathcal{Y}, x_*)
= \mathcal{N}
\left(
\phi(x_*)^\top \mathbf{m}_N,\;
\phi(x_*)^\top \mathbf{S}_N \phi(x_*) + \sigma^2
\right).
\]

This variance has two components:

1. \( \phi(x_*)^\top \mathbf{S}_N \phi(x_*) \): parameter uncertainty  
2. \( \sigma^2 \): observation noise  

Even with infinite data, noise remains.

</div>



<div class="example">
Noise-Free Function Values

For the latent function
\[
f(x_*) = \phi(x_*)^\top \mathbf{\theta},
\]
we have
\[
p(f(x_*) \mid \mathcal{X}, \mathcal{Y})
= \mathcal{N}
\left(
\phi(x_*)^\top \mathbf{m}_N,\;
\phi(x_*)^\top \mathbf{S}_N \phi(x_*)
\right).
\]
This is a **distribution over functions**, not just numbers.

</div>





---

### Marginal Likelihood

The **marginal likelihood** (or **evidence**) measures how well the model explains the observed data:
\[
p(\mathcal{Y} | \mathcal{X}) = \int p(\mathcal{Y} |\mathcal{X}, \mathbf{\theta}) p(\mathbf{\theta}) d\mathbf{\theta}.
\]
Using Gaussian conjugacy, the result is:
\[
p(\mathcal{Y} | \mathcal{X}) = \mathcal{N}(\mathbf{y} | \mathbf{X} \mathbf{m}_0, \mathbf{X} \mathbf{S}_0 \mathbf{X}^\top + \sigma^2 \mathbf{I}),
\]
with:
\[
\mathbb{E}[\mathcal{Y} | \mathcal{X}] = \mathbf{X} \mathbf{m}_0, \quad
\text{Cov}[\mathcal{Y} | \mathcal{X}] = \mathbf{X} \mathbf{S}_0 X^\top + \sigma^2 I
\]
The marginal likelihood is crucial for model comparison and selection, as it measures the model’s fit after integrating over parameter uncertainty.





<div class="example">
Marginal Likelihood (Evidence)

The marginal likelihood integrates out parameters:
\[
p(\mathcal{Y} \mid \mathcal{X})
= \int p(\mathcal{Y} \mid \mathcal{X}, \mathbf{\theta}) p(\mathbf{\theta}) d\mathbf{\theta}.
\]

Result:
\[
p(\mathcal{Y} \mid \mathcal{X})
= \mathcal{N}
(\mathcal{Y} \mid \mathbf{\Phi} \mathbf{m}_0,
\mathbf{\Phi} \mathbf{S}_0 \mathbf{\Phi}^\top + \sigma^2 \mathbf{I}).
\]

**Interpretation**:

- Measures how well the model explains data *on average* over parameters  
- Automatically penalizes overly flexible models  
- Central to Bayesian model selection  

</div>



---

#### Key Insights

| Concept | Description |
|----------|--------------|
| **Parameter Treatment** | Bayesian regression treats parameters \( \mathbf{\theta} \) as random variables. |
| **Prior → Posterior** | Prior \( p(\mathbf{\theta}) \) is updated using data to yield posterior \( p(\mathbf{\theta} | \mathcal{X}, \mathcal{Y}) \). |
| **Predictions** | Are averages over all plausible parameter values (not point estimates). |
| **Uncertainty** | Captured via predictive variance combining model and observation noise. |
| **Marginal Likelihood** | Used for model evidence and Bayesian model selection. |

---

#### Summary of Distributions

| Distribution | E\mathbf{x}pression | Description |
|---------------|-------------|-------------|
| **Prior** | \( p(\mathbf{\theta}) = \mathcal{N}(\mathbf{m}_0, \mathbf{S}_0) \) | Beliefs before seeing data |
| **Likelihood** | \( p(\mathcal{Y} | \mathcal{X}, \mathbf{\theta}) = \mathcal{N}(Y | \mathbf{\Phi} \mathbf{\theta}, \sigma^2 I) \) | Data model |
| **Posterior** | \( p(\mathbf{\theta} | \mathcal{X}, \mathcal{Y}) = \mathcal{N}(\mathbf{\theta} | \mathbf{m}_N, \mathbf{S}_N) \) | Updated parameter beliefs |
| **Predictive** | \( p(y_* | \mathcal{X}, \mathcal{Y}, \mathbf{x}_*) = \mathcal{N}(\phi(\mathbf{x}_*)^\top \mathbf{m}_N, \phi(\mathbf{x}_*)^\top \mathbf{S}_N \phi(\mathbf{x}_*) + \sigma^2) \) | Predictions at new points |
| **Marginal Likelihood** | \( p(\mathcal{Y} | \mathcal{X}) = \mathcal{N}(\mathcal{Y} | \mathbf{X} \mathbf{m}_0, \mathbf{X}  \mathbf{S}_0 \mathbf{X} ^\top + \sigma^2 \mathbf{I} ) \) | Model evidence |



---



### Exercises {.unnumbered .unlisted}


Put some exercises here.








##  Maximum Likelihood as Orthogonal Projection


Maximum Likelihood Estimation (MLE) in linear regression has a clear geometric interpretation.  It corresponds to the orthogonal projection of the target vector \( \mathbf{y} \) onto the subspace spanned by the input data.

---

### Simple Linear Regression Case

Consider the model:
\[
y = \mathbf{x}\mathbf{\theta} + \epsilon, \quad \epsilon \sim \mathcal{N}(0, \sigma^2).
\]
Given training data \( \{(x_1, y_1), \ldots, (x_N, y_N)\} \), the MLE for \( \theta \) is:
\[
\theta_{ML} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y} = \frac{\mathbf{X}^\top \mathbf{y}}{\mathbf{X}^\top \mathbf{X}}
\]
where:

- \( \mathbf{X} = [x_1, \ldots, x_N]^\top \in \mathbb{R}^N \),  
- \( \mathbf{y} = [y_1, \ldots, y_N]^\top \in \mathbb{R}^N \).  

This gives the fitted (reconstructed) outputs:
\[
\mathbf{X} \theta_{ML} = \frac{\mathbf{X}\mathbf{X}^\top}{\mathbf{X}^\top \mathbf{X}} \mathbf{y}
\]


Geometrically, we can view the matrix
  \[
  P = \frac{\mathbf{X}\mathbf{X}^\top}{\mathbf{X}^\top \mathbf{X}}.
  \]
as the projection matrix onto the one-dimensional subspace spanned by \( \mathbf{X} \).  The MLE solution \( \mathbf{X}\theta_{ML} \) is thus the orthogonal projection of \( \mathbf{y} \) onto this subspace.  This projection minimizes the squared distance between \( \mathbf{y} \) and the model predictions \( \mathbf{X}\theta \):
\[
\min_{\theta} \|\mathbf{y} - \mathbf{X}\theta\|^2.
\]
Hence, MLE not only gives the best statistical fit (in the least-squares sense) but also the geometrically optimal fit — the closest point to \( \mathbf{y} \) within the space of linear combinations of \( \mathbf{X} \).

---

### General Linear Regression Case

For the more general model:
\[
y = \phi(\mathbf{x})^\top \theta + \epsilon, \quad \epsilon \sim \mathcal{N}(0, \sigma^2),
\]
where \( \phi(\mathbf{x}) \in \mathbb{R}^K \) is a vector of feature functions, the results extend naturally.  The MLE estimate is:
\[
\theta_{ML} = (\Phi^\top \Phi)^{-1} \Phi^\top y,
\]
and the fitted outputs are:
\[
\mathbf{y} \approx \Phi \theta_{ML}.
\]

Here:

- \( \Phi \in \mathbb{R}^{N \times K} \) is the feature matrix,  
- The column space of \( \Phi \) defines a K-dimensional subspace of \( \mathbb{R}^N \),  
- The projection matrix is:
  \[
  P = \Phi (\Phi^\top \Phi)^{-1} \Phi^\top.
  \]
Thus, MLE corresponds to orthogonal projection of \( \mathbf{y} \) onto the subspace spanned by the columns of \( \Phi \).

---

### Special Case: Orthonormal Basis

If the feature vectors \( \phi_k \) form an orthonormal basis, then:
\[
\Phi^\top \Phi = I
\]
and the projection simplifies to:
\[
P = \Phi \Phi^\top = \sum_{k=1}^K \phi_k \phi_k^\top.
\]
In this case:

- The projection of \( \mathbf{y} \) is simply the sum of individual projections onto each basis vector.  
- The coupling between features disappears because of orthogonality.  



<div class="example"> 
Fourier bases and wavelets are examples of orthogonal bases commonly used in signal processing.
</div>







---



### Exercises {.unnumbered .unlisted}


Put some exercises here.