```python
import numpy as np
import matplotlib.pyplot as plt
import math
from numpy.linalg import inv
from numpy import matmul as mm
import scipy
from scipy.linalg import lu_factor,lu_solve
from scipy.linalg import cho_factor,cho_solve
from scipy.stats import multivariate_normal
```


```python
x = np.linspace(-4,2,20)
y = (x+1)**3+10*np.random.normal(0,1,20)
```

### Posterior distribution

$$p(\mathbf{w|x,y},\mathcal{M})=\frac{p(\mathbf{w}|\mathcal{M})p(\mathbf{y|x,w},\mathcal{M})}{p(\mathbf{y|x},\mathcal{M})}=\mathcal{N}(\mathbf{w|\mu,\Sigma})$$
where 

\begin{align}
\mathbf{\Sigma} &= (\sigma^{-2}_{noise}\mathbf{\Phi^\top\Phi}+\sigma^{-2}_{\mathbf{w}}\mathbf{I})^{-1}\\
\mathbf{\mu} &= (\mathbf{\Phi^\top\Phi}+\frac{\sigma^2_{noise}}{\sigma^2_\mathbf{w}}\mathbf{I})^{-1}\mathbf{\Phi^\top y}
\end{align}

### Predictive distribution

\begin{align}
p(y_*|x_*,\mathbf{x},\mathbf{y},\mathcal{M}) &= \int p(y_*|x_*,\mathbf{w},\mathcal{M})p(\mathbf{y}|\mathbf{x},\mathbf{w},\mathcal{M})d\mathbf{w}\\
&= \mathcal{N}(y_*|\mathbf{\phi}(x_*)^\top\mathbf{\mu},\mathbf{\phi}(x_*)^\top\mathbf{\Sigma}\mathbf{\phi}(x_*)+\sigma^2_{noise}\mathbf{I})
\end{align}

### Marginal likelihood (Evidence)

\begin{align}
p(\mathbf{y}|\mathbf{x},\mathcal{M})&= \int p(\mathbf{w}|\mathbf{x},\mathcal{M})p(\mathbf{y}|\mathbf{x},\mathbf{w},\mathcal{M})d\mathbf{w}\\
&= \mathcal{N}(\mathbf{y};\mathbf{0},\sigma^2_\mathbf{w}\mathbf{\Phi\Phi^\top}+\sigma^2_{noise}\mathbf{I})
\end{align}


```python
def BayesFitLU(x,y,order=3):
    for i in range(order+1):
        if i == 0:
            Phi = x**i
        else:
            Phi = np.vstack((Phi,x**i))
    Phi = Phi.T
    if order == 0:
        Phi = Phi.reshape(-1,1)
    margcov = mm(Phi,Phi.T)+100*np.eye(Phi.shape[0])
    lu,piv = lu_factor(mm(Phi.T,Phi)+100*np.eye(Phi.shape[1]))
    tmp = lu_solve((lu,piv),Phi.T)
    w = mm(tmp,y)
    f = mm(Phi,w)
    logmarglik = multivariate_normal.logpdf(y,None,margcov)
    
    covlu,covpiv = lu_factor(mm(Phi.T,Phi)/100+np.eye(Phi.shape[1]))
    return f,covlu,covpiv,w,logmarglik
```


```python
mllistLU = []
plt.figure(figsize=(20,20))
for i in range(10):
    plt.subplot(5,2,i+1)
    f,lu,piv,w,lml = BayesFitLU(x,y,i)
    if i == 3:
        print(w)
    mllistLU.append(lml)
    plt.plot(x,y,'b.')
    plt.plot(x,f,'r-')
    plt.title('Order '+str(i)+': '+str(lml))
```


```python
plt.plot(np.arange(10),mllistLU)
plt.title('Log-Marginal Evidence')
```


```python
def BayesPredictLU(x,mu,lu,piv,order=3):
    for i in range(order+1):
        if i == 0:
            Phi = x**i
        else:
            Phi = np.vstack((Phi,x**i))
    Phi = Phi.T
    if order == 0:
        Phi = Phi.reshape(-1,1)
    tmu = mm(Phi,mu)
    tcov = mm(Phi,lu_solve((lu,piv),Phi.T))+100*np.eye(Phi.shape[0])
    return tmu,tcov
```


```python
mu,lu,piv,w,lml = BayesFitLU(x,y,3)
xt = np.linspace(-6,4,100)
ytmu,ytcov = BayesPredictLU(xt,w,lu,piv,3)
```


```python
unc = np.sqrt(np.diag(ytcov))
plt.fill_between(xt,ytmu-unc,ytmu+unc,color='lime')
plt.plot(xt,ytmu,'r-')
plt.plot(x,y,'r.')
```

### Summary


```python
mllistLU = []
plt.figure(figsize=(10,15))
plt.subplot(3,1,1)
for i in range(10):
    f,lu,piv,w,lml = BayesFitLU(x,y,i)
    mllistLU.append(lml)
    plt.plot(x,y,'b.')
    plt.plot(x,f,'r-')

plt.subplot(3,1,2)
plt.plot(np.arange(10),mllistLU)
plt.title('Log-Marginal Evidence')

mu,lu,piv,w,lml = BayesFitLU(x,y,3)
xt = np.linspace(-6,4,100)
ytmu,ytcov = BayesPredictLU(xt,w,lu,piv,3)
plt.subplot(3,1,3)
unc = np.sqrt(np.diag(ytcov))
plt.fill_between(xt,ytmu-unc,ytmu+unc,color='lime')
plt.plot(xt,ytmu,'r-')
plt.plot(x,y,'r.')
plt.title('Order 3 log-evidence: '+str(lml))
```


```python

```

cd9af9a3d959681601114ea7b4ea93388a834355
ipynb
Jupyter Notebook
Probabilistic Machine Learning/Introduction to Probabilistic Machine Learning/Bayesian inference and prediction with finite regression models.ipynb
zcemycl/ProbabilisticPerspectiveMachineLearning
8291bc6cb935c5b5f9a88f7b436e6e42716c21ae
MIT
2019-11-20T10:20:29.000Z
2021-11-09T11:15:23.000Z
Probabilistic Machine Learning/Introduction to Probabilistic Machine Learning/.ipynb_checkpoints/Bayesian inference and prediction with finite regression models-checkpoint.ipynb
zcemycl/ProbabilisticPerspectiveMachineLearning
8291bc6cb935c5b5f9a88f7b436e6e42716c21ae
MIT
Probabilistic Machine Learning/Introduction to Probabilistic Machine Learning/.ipynb_checkpoints/Bayesian inference and prediction with finite regression models-checkpoint.ipynb
zcemycl/ProbabilisticPerspectiveMachineLearning
8291bc6cb935c5b5f9a88f7b436e6e42716c21ae
MIT
2020-05-27T03:56:38.000Z
2021-05-02T13:15:42.000Z
Qwen/Qwen-72B
1. YES
2. YES