```python
from Model_GPR import Model,GPR,gamma,normal,RBF_2
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import animation, cm
import pandas as pd
```

## Generate Data


```python
import numpy as np
np.random.seed(404) 
noise = 0.25
x, y = np.arange(-5, 5, 0.3), np.arange(-5, 5, 0.3)
xx, yy = np.meshgrid(x, y)
X_2D = np.column_stack((xx.reshape(-1),yy.reshape(-1)))
X_obs = np.random.uniform(-4, 4, (100, 2))
Y_obs = np.sin(np.linalg.norm(X_obs, axis=1))-np.tanh(0.5*np.linalg.norm(X_obs,axis=1)) + \
             noise * np.random.randn(len(X_obs))
```


```python
fig_1 ,(ax1,ax2) = plt.subplots(1,2,figsize=(10,4))
ax1.scatter(X_obs[:,0],Y_obs, c=Y_obs, cmap=cm.jet)
ax1.set_xlabel('X1')
ax1.set_ylabel('Y')
ax2.scatter(X_obs[:,1],Y_obs,c=Y_obs,cmap=cm.jet)
ax2.set_xlabel('X2');
```


```python
from MLE import MLE,plot_MLE_posterior_dist
from scipy.optimize import minimize
from numpy.linalg import cholesky, det
from scipy.linalg import solve_triangular
```

## Without parameter optimization


```python
fig = plt.figure(figsize=(7,7))
plot_MLE_posterior_dist(xx,yy,X_obs,Y_obs,X_2D,[1.0,1.0,0.3],"Before parameter optimization: l={} sigma={} noise={}".format(1.0,1.0,0.3))
```

## Let's start with an MLE of the parameters

We ne to define an objective function, which is the negative loglikelihood of the data

\begin{equation}
 -l(y|\theta)= \frac{1}{2}(y-\mu)^{T}K_{y}(\theta)^{-1}(y-\mu)+ \frac{m}{2}\log{2\pi}+ \frac{1}{2}\log{|K_{y}(\theta)|}.
\end{equation}

a more stable version following http://www.gaussianprocess.org/gpml/chapters/RW.pdf (algorithm 2.1) is

\begin{align}
    L &= \text{cholesky}(K_y(\theta))\\
    \alpha &= L^{-T}(L^{-1}y)\\
    -l(y|\theta) &=\frac{1}{2}y^{T}\alpha + \sum_{i}^{m}\log{L_{ii}} + \frac{m}{2}\log{2\pi}
\end{align}

which holds for a zero mean Gaussian process.


```python
res_MLE = minimize(MLE(X_obs, Y_obs),[1.0,1.0,0.3], 
               bounds=((1e-5,None),(1e-5,None),(1e-5,None)),
               method='L-BFGS-B')
res_MLE
```




          fun: 40.60040642011458
     hess_inv: <3x3 LbfgsInvHessProduct with dtype=float64>
          jac: array([ 0.00406146, -0.00175433, -0.0138165 ])
      message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
         nfev: 72
          nit: 14
       status: 0
      success: True
            x: array([2.18534535, 1.35116544, 0.23744747])




```python
fig = plt.figure(figsize=(7,7))
plot_MLE_posterior_dist(xx,yy,X_obs,Y_obs,X_2D,res_MLE.x,"After parameter optimization (MLE): l={} sigma={} noise={}".format(np.round(res_MLE.x[0],2),np.round(res_MLE.x[1],2),np.round(res_MLE.x[2],2)))
```

## Define Model

We use three gamma priors for the unkown hyperparameters $\ell$, $\sigma$ and noise level $\sigma_{\text{noise}}$

\begin{align}
    \ell \sim \Gamma(2.0,3.0)\\
    \sigma  \sim \Gamma(2.0,3.0)\\
    \sigma\_{\text{noise}}\sim \Gamma(2.0,3.0)
\end{align}

and define the model chronologically. This means we start by defining the prior assumptions for the hyperparameters and noise level. After doing this we choose a covariance function - in this example we use the RBF Kernel to model the covariance between the input points. We can print a summary of the defined model (optional).
The last step is to use the MAP function of the model to find the optimized parameters. The initial guess was arbitrairly set to 1,1,0.3.


```python
test_model = GPR("Test_1")
with test_model as model_2:
    model_2.add_prior("ℓ",gamma(2,3))
    model_2.add_prior("σ",gamma(2,3))
    model_2.add_prior("noise", gamma(2,3))
    model_2.add_cov(RBF_2)
    model_2.summarize()
    res = model_2.MAP(X_obs,Y_obs,theta=[1, 1,0.3])
print(res)
```

    model:Test_1
     self.prior :{'ℓ': <Model_GPR.gamma object at 0x0000027074A4C5F8>, 'σ': <Model_GPR.gamma object at 0x0000027074A4CD30>, 'noise': <Model_GPR.gamma object at 0x0000027074A4CBE0>}
          fun: 45.25451798676848
     hess_inv: <3x3 LbfgsInvHessProduct with dtype=float64>
          jac: array([ 1.69109171e-04,  2.77111667e-05, -4.49063009e-04])
      message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
         nfev: 60
          nit: 11
       status: 0
      success: True
            x: array([2.03257791, 1.17100866, 0.23610388])


## The posterior distribution with optimized parameters


```python
model_2.posterior_dist(X_2D,X_obs,Y_obs)
fig_2 = plt.figure(figsize=(7,7))
model_2.plot_post_dist(xx,yy,X_obs,Y_obs,"After parameter optimization (MAP): l={} sigma={} noise={}".format(np.round(res.x[0],2),np.round(res.x[1],2),np.round(res.x[2],2)))
```

## The variation of the parameters while optimization


```python
model_2.theta_df.plot(subplots=True,layout=(3,1),figsize=(10,10), color='dodgerblue');
```

## Average trainings loss for 10 function evaluations


```python
fig_3,ax = plt.subplots(figsize=(10,7))
pd.DataFrame(model_2.training_loss, columns=['logp']).rolling(10).mean()[10:].plot(color='dodgerblue', ax=ax)
ax.set_yscale('log')
```

## Second approach with pymc3


```python
import pymc3 as pm
import theano
import theano.tensor as tt
```

We define the same model as before. Again we start by defining the prior distributions for the kernel hyperparameter. After doing this we select a covariance function. In contrast to above, we first need to define a Gaussian process with the selected covariance function. Note that the covariance function itself stores the prior assumptions as random variables. Then we define a prior distribution for the noise level and finally construct the loglikelihood of the modell. Finally we can use the find_MAP method to optimize the parameters.


```python
with pm.Model() as model:
    ℓ = pm.Gamma("ℓ", alpha=2.0, beta=3.0)
    σ = pm.Gamma("σ", alpha=2.0,beta=3.0)
    
    cov = σ ** 2 * pm.gp.cov.ExpQuad(2, ℓ)
    gp = pm.gp.Marginal(cov_func=cov)

    σ_noise = pm.Gamma("σ_noise", alpha=2.0,beta=3.0)
    y_ = gp.marginal_likelihood("y", X=X_obs, y=Y_obs, noise=σ_noise)

    mp = pm.find_MAP(method='BFGS')
```



<div>
    <style>
        /* Turns off some styling */
        progress {
            /* gets rid of default border in Firefox and Opera. */
            border: none;
            /* Needs to be in here for Safari polyfill so background images work as expected. */
            background-size: auto;
        }
        .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {
            background: #F44336;
        }
    </style>
  <progress value='15' class='' max='15' style='width:300px; height:20px; vertical-align: middle;'></progress>
  100.00% [15/15 00:00<00:00 logp = -45.256, ||grad|| = 0.72597]
</div>



    


We can now check if the pymc3 optimized hyperparameter lead to any significant differences to the optimized solution with the help of our model.


```python
np.array([mp['ℓ'], mp['σ'], mp['σ_noise']])-model_2.theta
```




    array([-6.03719906e-06, -4.55456266e-06,  7.65385917e-08])



## MCMC parameter sampling


```python
with model:
    trace = pm.sample(1000, chains=2)
```

    C:\ProgramData\Anaconda3\lib\site-packages\pymc3\sampling.py:468: FutureWarning: In an upcoming release, pm.sample will return an `arviz.InferenceData` object instead of a `MultiTrace` by default. You can pass return_inferencedata=True or return_inferencedata=False to be safe and silence this warning.
      FutureWarning,
    Auto-assigning NUTS sampler...
    Initializing NUTS using jitter+adapt_diag...
    Multiprocess sampling (2 chains in 2 jobs)
    NUTS: [σ_noise, σ, ℓ]




<div>
    <style>
        /* Turns off some styling */
        progress {
            /* gets rid of default border in Firefox and Opera. */
            border: none;
            /* Needs to be in here for Safari polyfill so background images work as expected. */
            background-size: auto;
        }
        .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {
            background: #F44336;
        }
    </style>
  <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>
  100.00% [4000/4000 00:26<00:00 Sampling 2 chains, 0 divergences]
</div>



    Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 41 seconds.



```python
import arviz as az
```


```python
with model:
    az.plot_trace(trace, var_names=["ℓ", "σ", "σ_noise"]);
    display(az.summary(trace, round_to=2))
```


```python

```

2d1277f52626c2f01e854ac4680724b86477056d
ipynb
Jupyter Notebook
Example_1.ipynb
jsve96/MA_GPR
adffdb6326d0cbde638995bde5ea031e4c75d408
MIT
Example_1.ipynb
jsve96/MA_GPR
adffdb6326d0cbde638995bde5ea031e4c75d408
MIT
Example_1.ipynb
jsve96/MA_GPR
adffdb6326d0cbde638995bde5ea031e4c75d408
MIT
Qwen/Qwen-72B
1. YES
2. YES