# Variational Bayesian Inference

In this notebook, we will review the variational Bayes process, beginning with a technical introduction to the formalism and derivation, followed by a python implementation. 
The material covered here references Blei et al., 2018, Varitional Inference: A Review for Statisticians and Chappel et al., 2016. The FMRIB Variational Bayes Tutorial. 


## Technical Overview

### Problem Statement 

Similarly to problems addressed by sampling methods, the goal of variational inference (VI) is to approximate parameter distributions from data, specifically in cases where an analytical treatment is intractable. 

Consider the following example (from Blei et al., 2018): 

For latent variabels $\mathbf{z} = z_{1:m}$ and observations $\mathbf{x} = x_{1:n}$, the posterior conditional density is given by : 

$$ p(\mathbf{z} \mid \mathbf{x}) = \frac{p(\mathbf{x} \mid \mathbf{z}) * p(\mathbf{z})}{p(\mathbf{x})} $$

The denominator, ${p(\mathbf{x})}$, whose value is needed to compute the posterior, is calculated by: 

$$ {p(\mathbf{x})} = \int p(\mathbf{z} , \mathbf{x}) dz $$


This integral is often intractable or too computationally expensive to be feasible. For similar reasons (notably the number of latent variables)sampling methods are slow to converge.

VI aims to circumvent the large time complexity by approaching the problem through optimisation. 

The process begins by positing a contrived _approximate_ density, $\mathfrak{D}$ of latent variables $\mathbf{z}$. Then, using this density can find a set of valuers for $q(\mathbf{z}) \in \mathfrak{D}$ whose values maximise the Kullback-Liebler divergence between the approximate density and the true posterior. 

$$ q^{*}(\mathbf{z}) = \underset{q(\mathbf{z}) \in \mathfrak{D}}{argmin} \mathrm{KL}(q(\mathbf{z}) \mid \mid p(\mathbf{z} \mid \mathbf{x})) $$


### VI Computation 

To compute the optimsation problem stated above requires we do manipulation of the terms. This is because the problem is still dependent on the _evidence_ term ($p(\mathbf{x})$; recall, the intragral of this is often difficult/intractable). To see this clearly, we can rewrite out expression using Bayes rule: 

$$ q^{*}(\mathbf{z}) = \underset{q(\mathbf{z}) \in \mathfrak{D}}{argmin} \mathrm{KL}\bigg(q(\mathbf{z}) \mid \mid \frac{p(\mathbf{z}, \mathbf{x})}{p(\mathbf{x})}\bigg) $$

We can avoid evidence term using the following manipulations. First, let's break down the KL divergence, recalling that the KL divergence can otherwise be expressed as the expectated value of the first argument minus the expected value of the second argument: 

$$ KL(q(\mathbf{z}) \mid \mid p(\mathbf{z} \mid \mathbf{x})) = \mathbb{E}[\log q(\mathbf{z})] - \mathbb{E}[\log p(\mathbf{z} , \mathbf{x})] + \log p(\mathbf{x}) $$

We formulate the Evidence Lower Bound (ELBO) by multiplying both sides by -1: 


\begin{align}
 ELBO(q) &= - KL(q(\mathbf{z}) \mid \mid p(\mathbf{z} \mid \mathbf{x})) + \log p(\mathbf{x}) \\
 
         &= \mathbb{E}[\log p(\mathbf{z}, \mathbf{x})] - \mathbb{E}[\log q(\mathbf{z})]
\end{align}

Note that, maximsing the ELBO is the same as minimising the KL divergence, since the evidence is constant and additive. 

By exanding the joint probaibility, we can rewrite the ELBO in terms of the log likelihood and KL divergence between our prior $p(\mathbf{z})$ and our variational density $q(\mathbf{z})$: 

\begin{align}
ELBO(q) &= \mathbb{E}[\log p(\mathbf{x} \mid \mathbf{z})] + \mathbb{E}[\log p(\mathbf{z})] - \mathbb{E}[\log q(\mathbf{z})] \\

&= \mathbb{E}[\log p(\mathbf{x} \mid \mathbf{z})] - KL(q(\mathbf{z}) \mid \mid p(\mathbf{z}))
\end{align}

By maximising the ELBO or minimising the KL divergence, we can solve the optimisation problem. Now, we need a method to do this.

### Optimisation 

There are many method available for solving the optimsation problem. Here we appeal to variational calculus (hence, variational inference). 

Full derivation to follow... 

Optimisation equates to updating parameters w.r.t the expectation of all other parameters given the current posterior:

$$ \log q(z_i) \approx \int q(z_{-i}) \log p(x \mid z) p(w) dw_{-i}$$ 


## Example 

We will work with an exponential decay model to explore the variational Bayes framework and derive further equations needed. 

### Model Definition 

$$ y = ae^{-\lambda t} + \epsilon $$
$$ \epsilon = \mathcal{N}(0,\phi^{1}) $$

We can define this using a simple lambda function and some parameter values. 


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


```python
def f(theta, t): 
    """ Sipmle exponential decay model

    args: 
    theta : array, float
            parameters for function
        t : array, float 
            timesteps to evaluate at
    """
    a, l = theta
    f1 = a * np.exp(-l * t)

    return f1
```


```python
a = 42
l = 1.5
theta = [a, l]
t = np.linspace(0,20,100)

phi = 100 
```


```python
model = f(theta, t)
data = f(theta, t) + np.random.normal(0, np.sqrt(1/10), t.shape[0])
```


```python
plt.plot(t, model, linewidth = 3, alpha = 0.6)
plt.scatter(t, data, c = 'orange')
plt.show()
```

#### Forward Model Specification

We can set up a forard model in the following way:

$$ \mathbf{y} = \mathbf{g} + \epsilon $$ 

$$ \mathbf{g} = \{g(\theta_{1}), . . ., g(\theta_{k})\} $$

$$ \epsilon \approx \mathcal{N}(0, \Phi^{-1}) $$ 

$$ \Phi = \{\phi_1, . . ., \phi_n\} $$

Where $\mathbf{y}$ is a vector of N observation assumed to be generated by a nonlinear forward function $g(\theta) + \epsilon $ with $k$ parameters, \theta, and $k$ precision parameters, $\phi$. 



#### Inference Specification

Using this forward model, we can parameterise the log-likelihood as multivariate normal distribution, with the data conditioned on parameters:

$$ p(\mathbf{y} \mid \mathbf{\Theta}) = \frac{N}{2} \log(\Phi) - \frac{1}{2}(\mathbf{y} - \mathbf{g})^{T} \Phi (\mathbf{y} - \mathbf{g}) $$

Where $\Theta$ is the set of independent parameters. We factorise the approximate posterior accordingly: 

$$ q(\Theta \mid \mathbf{y}) = q(\theta_n)q(\phi_n) $$ 

and choose conjucate distributions for their form: 

$$ q(\theta \mid \mathbf{y}) \approx MVN(\theta; \mathbf{m}, \Lambda^{-1}) $$ 

$$ q(\phi \mid \mathbf{y}) \approx Ga(\phi; s, c) $$ 

Using these, we can construct the log posterior and thus update equations. 

\begin{align}
L &= \log p(\mathbf{y} \mid \mathbf{\Theta}) + \log p(\theta) + \log p(\phi) \\
&= \underbrace{\frac{N}{2} \log(\Phi) - \frac{1}{2}(\mathbf{y} - \mathbf{g})^{T} \Phi (\mathbf{y} - \mathbf{g})}_{\log p(\mathbf{y} \mid \mathbf{\Theta})} 
- \underbrace{\frac{1}{2}(\theta - \mathbf{m_0})^{T} \Lambda_0 (\theta - \mathbf{m_0})}_{\log p(\theta)}
+ \underbrace{(c_0 - 1)\log\phi - \frac{1}{s_0} \phi}_{\log p(\phi)}
+ \underbrace{const\{\theta, \phi\}}_{\text{constant terms}}
\end{align}


##### Compuational side note: 

We will need to make an approximation to the function $g(\theta)$ for mathematical and computational convenience. Simply, we will take a first order Taylor expansion around the mean of the posterior distribution. 

$$ g(\theta) \approx g(\mathbf{m}) + \mathbf{J}(\theta - \mathbf{m}) $$ 

Where $\mathbf{m}$ is the mean value and $\mathbf{J}$ is the Jacobian matrix. 

In practice, we will compute this approximation using a simply central difference scheme. Let's define some functions to do this for us. 


```python
def time_step(theta): 
    """Control time stepping interval proportional to the parameter magnitiude. 

    args: 
    theta:: float 
    """

    delta = theta * 1e-5

    if delta < 0:
            delta = -delta
    if delta < 1e-10:
            delta = 1e-10
    
    return delta

def central_difference(i, theta, delta, t):
    """ Calculate the derivate using a first order central difference approximation

    args: 
        i : int 
            parameter index
    theta : np.array float 
            array of paramters for function, f
    delta : float 
            time-step for central difference
        t : np.array float
            time steps to evaluate at
    returns: 
        df : np.array, float 
             first order approximation, df, for time steps, t
    """
    dtheta = np.array(theta), np.array(theta)

    dtheta[0][i] += delta
    dtheta[1][i] -= delta
    df = (f(dtheta[0], t) - f(dtheta[1], t)) / (2 * delta)
    return df

def Jacobian(theta, t):
    """Compute the Jacobian for globally defined function, f, with parameter set Theta 

    args:
    theta : np.array, float
            parameters of global function, f
        t : np.array, float 
            time steps to evaluate at

    returns: 
        J : np.array, float
            Jacobian vector/matrix evaluated at theta, t. 
    """
    J = None
    #f_n = len(f)
    p_n = len(theta)

    for i in range(p_n):
        delta = time_step(theta[i])
        if J is None:
                J = np.zeros([len(t), len(theta)], dtype=np.float32)
        J[:,i] = central_difference(i, theta, delta, t)
    
    return J
```

### Deriving updates

For our model, we will need to derive updates on our model parameters $\theta$ and noise parameters $phi$. 

Recall the the derivation for our updates w.r.t. the free energy functional: 

$$ p(\theta_{i}) = \int p(\Theta \mid y) p(\theta_{-i}) d\theta_{-i} $$ 

Where $\theta_{-i}$ are all parameters that are not $\theta$. In our case, we have two parameters $\theta$ and $\phi$. Which will be given by: 

$$ p(\theta) = \int p(\Theta \mid y) p(\phi) d\phi $$ 

and 

$$ p(\phi) = \int p(\Theta \mid y) p(\theta) d\theta $$


#### Update on theta 

We will first derive our update on model parameters $\theta$.

Recall the forms of the log posterior and parameters. With these, we can write down the full form of the update equation: 

$$ q(\theta) =  \int \bigg( \frac{N}{2} \log(\phi) - \frac{1}{2}(\mathbf{y} - \mathbf{g})^{T} \Phi (\mathbf{y} - \mathbf{g}) - \frac{1}{2}(\theta - \mathbf{m_0})^{T} \Lambda_0 (\theta - \mathbf{m_0}) + (c_0 - 1)\log\phi - \frac{1}{s_0} \phi + const\{\theta, \phi\} \bigg) Ga(\phi; s, c) d\phi $$

Where $\mathbf{g}$ is an abbreviation for the function $g(\theta)$. We can start simplifying this by first removing terms from that are constant w.r.t. $\theta$ (remember, we are trying to optimise $\theta$) and factoring out terms in the intregral constant w.r.t. $\phi$. 

In so doing, we arrive at the following: 
$$ q(\theta) =  - \frac{1}{2}(\theta - \mathbf{m_0})^{T} \Lambda_0 (\theta - \mathbf{m_0}) - \frac{1}{2}(\mathbf{y} - \mathbf{g})^{T}(\mathbf{y} - \mathbf{g}) \int \phi Ga(\phi; s, c) d\phi + const{\theta} $$ 

The integral term corresponds to the mode of the Gamme distribution, which, since we are using precision (the inverse of variance), is equal to $sc$. Thus, we have: 

$$ q(\theta) =  - \frac{1}{2}(\theta - \mathbf{m_0})^{T} \Lambda_0 (\theta - \mathbf{m_0}) - \frac{1}{2}sc(\mathbf{y} - \mathbf{g})^{T}(\mathbf{y} - \mathbf{g}) + const{\theta} $$

To complete the derivation, we substitute the function $\mathbf{g}$ for it's Taylor approximation, given above and restated here. 

$$ g(\theta) \approx g(\mathbf{m}) + \mathbf{J}(\theta - \mathbf{m}) $$ 
$$ \mathbf{y} - \mathbf{g} = \mathbf{k} + \mathbf{J}(\theta - \mathbf{m})$$ 

Substituting this in and droppping the constant term, we obtain: 
$$ q(\theta) = - \frac{1}{2}(\theta - \mathbf{m_0})^{T} \Lambda_0 (\theta - \mathbf{m_0}) - \frac{1}{2}sc(\mathbf{k} + \mathbf{J}(\theta - \mathbf{m}))^{T}(\mathbf{k} + \mathbf{J}(\theta - \mathbf{m})$$

Expanding this out: 
$$q(\theta) = -\frac{1}{2} \bigg( \theta^{T}\Lambda_0\theta + \theta^{T}\Lambda_0 \mathbf{m_0} + m_0^{T}\Lambda_0\theta + \mathbf{m_0}\Lambda_0 \mathbf{m_0} sc (\theta^{T}\mathbf{J}^{T}\mathbf{J}\theta + \theta^{T}\mathbf{J}^{T}\mathbf{k} - \theta^{T}\mathbf{J}^{T}\mathbf{Jm} + \mathbf{kJ}\theta + \mathbf{m}^{T}\mathbf{J}^{T}\mathbf{J}\theta) \bigg) $$

We can then exlude terms that do not depend on $\theta$ and group like terms together

$$ = -\frac{1}{2} \bigg(\theta^{T}\Lambda_0\theta + sc \theta^{T}\mathbf{J}^{T}\mathbf{J}\theta  + \theta^{T}\Lambda_0 \mathbf{m_0} + \theta^{T}\mathbf{J}^{T}k - \theta^{T}\mathbf{J}^{T}\mathbf{Jm} + \mathbf{m_0}^{T}\Lambda_0\theta + \mathbf{kJ}\theta + \mathbf{m}^{T}\mathbf{J}^{T}\mathbf{J}\theta \bigg) $$

Which can be factorised as: 

$$ = -\frac{1}{2} \bigg(\theta^T(\Lambda_0 + dc \mathbf{J}^T \mathbf{J})\theta - \theta^T(\Lambda_0 \mathbf{m_0} + sc\mathbf{J}(\mathbf{k + Jm})) -(\Lambda_0 \mathbf{m_0} + sc\mathbf{J}(\mathbf{k + Jm})^T \theta \bigg) $$

Remember we specify the log posterior as a multivariate normal distribution. Thus, analysing the corresponding coefficients given by the multivariate normal form, we arrive at the update equations: 

$$ 
\Lambda = sc \mathbf{J}^T \mathbf{J} + \Lambda_0
\\ 
\\
\text{and}
\\
\\ 
\mathbf{m} = \Lambda^{-1} ( sc \mathbf{J}^T (\mathbf{k + Jm}) + \Lambda_0 \mathbf{m_0}) 
$$

We can write a simple function to compute this update:


```python
def parameter_update(error, params, priors, J):
    """ Update forward model function parameters theta in accordance with the update equations above
    
    args: 
        error : array, float
                vector of the difference between observations and model prediction 
       params : tuple
                parameters values
       priors : tuple 
                priors values
            J : array, float 
                array of Jacbian values from the model given parameter updates (calculated above)
    returns: 
       params : tuple
                updated parameters values
    """
    m, p, s, c = params
    m0, p0, _, _ = priors

    p_new = s*c*np.dot(J.transpose(), J) + p0
    c_new = np.linalg.inv(p_new)
    m_new = np.dot(c_new, (s * c * np.dot(J.transpose(), (error +    np.dot(J, m))) + np.dot(p0, m0)))
    
    params[0][:], params[1][:] = m_new, p_new

    return params
```

#### Update on noise

We now move on to deriving update equations for noise. Expanding on the form above, the full form of the update equation will be: 


$$ q(\phi) =  \int \bigg( \frac{N}{2} \log(\phi) - \frac{1}{2}(\mathbf{y} - \mathbf{g})^{T} \Phi (\mathbf{y} - \mathbf{g}) - \frac{1}{2}(\theta - \mathbf{m_0})^{T} \Lambda_0 (\theta - \mathbf{m_0}) + (c_0 - 1)\log\phi - \frac{1}{s_0} \phi + const\{\theta, \phi\} \bigg) MVN(\theta; m, \Lambda) d\theta $$

As before, we can remove terms from the integral that are constant w.r.t. $\theta$ and drop terms that are constant w.r.t. $\phi$. Thus we have: 

$$ = \bigg( \frac{N}{2} + c_0 -1 \bigg)\log\phi - \frac{1}{s_0} \phi - \frac{1}{2} \phi \int (\mathbf{y} - \mathbf{g})^{T}(\mathbf{y} - \mathbf{g}) MVN(\theta; m, \Lambda) d\theta $$

Swapping in our linearisation for $g(\theta)$ and expanding out, we obtain: 

$$ = \bigg( \frac{N}{2} + c_0 -1 \bigg)\log\phi - \frac{1}{s_0} \phi - \frac{1}{2} \phi \int \bigg( \mathbf{k}^T\mathbf{k} - (\theta - \mathbf{m})^T\mathbf{J}^T\mathbf{k} - \mathbf{k}^T\mathbf{J}(\theta - \mathbf{m}) + (\theta - \mathbf{m})^T\mathbf{J}^T\mathbf{J}(\theta - \mathbf{m})\bigg)MVN(\theta; m, \Lambda) d\theta $$

When expanded out with the multivariate normal on $\theta$, the integrals of the secona and third terms cancel out, since it generates an expectation on the MVN (equalling the mean), m is constant with $d\theta$ and the integral of the MVN is equal 1; thus, we have $m - m = 0$. We also make use of the general rule result: 

$$ \int (\theta - \mathbf{m})^T\mathbf{U}(\theta - \mathbf{m})MVN(\theta; m, \Lambda) d\theta = Tr(\Lambda^{-1} U)  $$ 

Therefore, the result of intial integral is equal to: 
$$ q(\phi) = \bigg( \frac{N}{2} + c_0 -1 \bigg)\log\phi - \frac{1}{s_0} \phi - \frac{1}{2} \phi \bigg(\mathbf{k}^T\mathbf{k} + Tr(\Lambda^{-1} \mathbf{J}^T \mathbf{J})\bigg) $$ 


Comparing this result with the form of a gamma distribution, we can compare coefficients for its moments and see that the update equations will be: 

$$ c = \frac{N}{2} + c_0 $$ 
$$ \frac{1}{s} = \frac{1}{s_0} + \frac{1}{2}\mathbf{k}^T\mathbf{k} + Tr(\Lambda^{-1} \mathbf{J}^T \mathbf{J}) $$

Once again, we can compte a simple function to compute this. 


```python
def noise_update(error, params, priors, J):
    """ Update forward model function parameters phi in accordance with the update equations above
    
    args: 
        error : array, float
                vector of the difference between observations and model prediction 
       params : tuple
                parameters values
       priors : tuple 
                priors values
            J : array, float 
                array of Jacbian values from the model given parameter updates (calculated above)
    returns: 
       params : tuple
                updated parameters values
    """
    _, p, _, _ = params 
    _, _, s0, c0 = priors

    N = len(data)

    c = np.linalg.inv(p)
    c_new = N/2 + c0
    s_new = 1/(1/s0 + 1/2 * np.dot(error.transpose(), error) + 1/2 * np.trace(np.dot(c, np.dot(J.transpose(), J))))
    
    params[2][:], params[3][:] = c_new, s_new

    return params
```


```python
def error_update(y, theta, t):
    """Calculate difference between data and model with updated parameters

    args:
        y : array, float 
            vector of noisy data (observations)
    theta : array, float
            vector of parameter values for model, f
        t : array, float 
            vector of time steps at which to evaluate model 
    
    returns:
    error : array, float 
            vector of difference between noisy data and updated model
    """


    error = y - f(theta[:2], t)
    
    return error
```


```python
def fit(data, params, priors, t, n): 
    
    theta = np.zeros((n, len(params[0])))
    for i in range(n):
        theta[i,:] = params[0]
        error = error_update(data, theta[i,:], t)
        
        J = Jacobian(theta[i,:], t)
        params = parameter_update(error, params, priors, J)
        params = noise_update(error, params, priors, J)
        print(params)

    return params, theta

```


```python
m = np.array([1.0, 1.0])
C = np.array([[1.0, 0], [0.0, 1.0]])
p = np.linalg.inv(C)
c = np.array([1e-8])
s = np.array([50.0])

init_params = m, p, c, s
```


```python
m0 = np.array([1.0, 1.0])
a_var0, lam_var0 = np.array([1000, 1.0])
C0 = np.array([[a_var0, 0], [0, lam_var0]])
p0 = np.linalg.inv(C0)
 

beta_mean0 = 1
beta_var0 = 1000
# c=scale, s=shape parameters for Gamma distribution
c0 = beta_var0 / beta_mean0
s0 = beta_mean0**2 / beta_var0

priors = m0, p0, c0, s0 
```


```python
params, theta_n = fit(data, init_params, priors, t, n=100)
```

    (array([1.05144648, 0.99998455]), array([[ 1.00150429e-03, -6.10400558e-07],
           [-6.10400558e-07,  1.00000062e+00]]), array([50.001]), array([0.00029803]))
    (array([34.56800553,  1.15395915]), array([[ 0.04583344, -0.0191285 ],
           [-0.0191285 ,  1.02038578]]), array([50.001]), array([0.00053172]))
    (array([40.64869248,  1.44492058]), array([[ 0.07234559, -0.83878455],
           [-0.83878455, 26.58017469]]), array([50.001]), array([0.01089733]))
    (array([42.28314592,  1.52020275]), array([[  1.23311602, -12.76151673],
           [-12.76151673, 370.1479491 ]]), array([50.001]), array([0.1073437]))
    (array([42.36458576,  1.52407141]), array([[  11.69599848, -117.77508799],
           [-117.77508799, 3379.16226228]]), array([50.001]), array([0.18443964]))
    (array([42.36780555,  1.52424174]), array([[  20.05859863, -201.69264506],
           [-201.69264506, 5783.49832491]]), array([50.001]), array([0.1873373]))
    (array([42.36789104,  1.52424657]), array([[  20.37206743, -204.82975272],
           [-204.82975272, 5873.26545289]]), array([50.001]), array([0.1873973]))
    (array([42.36787601,  1.52424564]), array([[  20.37854732, -204.89440961],
           [-204.89440961, 5875.11537456]]), array([50.001]), array([0.18739853]))
    (array([42.36789893,  1.52424693]), array([[  20.37868699, -204.89597897],
           [-204.89597897, 5875.15942469]]), array([50.001]), array([0.18739855]))
    (array([42.36789078,  1.52424619]), array([[  20.3786758 , -204.89571521],
           [-204.89571521, 5875.1520541 ]]), array([50.001]), array([0.18739855]))
    (array([42.36788439,  1.52424595]), array([[  20.37868544, -204.89590102],
           [-204.89590102, 5875.15912025]]), array([50.001]), array([0.18739855]))
    (array([42.36788341,  1.52424597]), array([[  20.37868739, -204.89591606],
           [-204.89591606, 5875.1596111 ]]), array([50.001]), array([0.18739855]))
    (array([42.36788275,  1.52424588]), array([[  20.37868723, -204.89591442],
           [-204.89591442, 5875.15956415]]), array([50.001]), array([0.18739855]))
    (array([42.36788189,  1.52424574]), array([[  20.37868721, -204.89591426],
           [-204.89591426, 5875.1601314 ]]), array([50.001]), array([0.18739855]))
    (array([42.36789461,  1.52424654]), array([[  20.37868718, -204.89598545],
           [-204.89598545, 5875.16126659]]), array([50.001]), array([0.18739855]))
    (array([42.36788603,  1.52424637]), array([[  20.37868267, -204.89582414],
           [-204.89582414, 5875.15496248]]), array([50.001]), array([0.18739855]))
    (array([42.36788979,  1.52424638]), array([[  20.37868305, -204.895828  ],
           [-204.895828  , 5875.15507304]]), array([50.001]), array([0.18739855]))
    (array([42.36789127,  1.52424648]), array([[  20.37868288, -204.89584412],
           [-204.89584412, 5875.1555949 ]]), array([50.001]), array([0.18739855]))
    (array([42.36789314,  1.52424666]), array([[  20.37868294, -204.89584475],
           [-204.89584475, 5875.15504102]]), array([50.001]), array([0.18739855]))
    (array([42.36788   ,  1.52424599]), array([[  20.37868299, -204.89577375],
           [-204.89577375, 5875.15333914]]), array([50.001]), array([0.18739855]))
    (array([42.36788298,  1.52424591]), array([[  20.37868753, -204.89589957],
           [-204.89589957, 5875.15850675]]), array([50.001]), array([0.18739855]))
    (array([42.36788158,  1.52424569]), array([[  20.37868716, -204.89591378],
           [-204.89591378, 5875.16011772]]), array([50.001]), array([0.18739855]))
    (array([42.36789744,  1.52424691]), array([[  20.37868942, -204.89602129],
           [-204.89602129, 5875.16126927]]), array([50.001]), array([0.18739855]))
    (array([42.3678873 ,  1.52424594]), array([[  20.37867595, -204.89569888],
           [-204.89569888, 5875.15209841]]), array([50.001]), array([0.18739855]))
    (array([42.36788342,  1.52424588]), array([[  20.37868764, -204.89593648],
           [-204.89593648, 5875.16082799]]), array([50.001]), array([0.18739855]))
    (array([42.3678822 ,  1.52424578]), array([[  20.37868728, -204.89591494],
           [-204.89591494, 5875.16015096]]), array([50.001]), array([0.18739855]))
    (array([42.36789267,  1.52424644]), array([[  20.37868719, -204.8959677 ],
           [-204.8959677 , 5875.16069824]]), array([50.001]), array([0.18739855]))
    (array([42.36789481,  1.52424675]), array([[  20.37868268, -204.89586002],
           [-204.89586002, 5875.15553835]]), array([50.001]), array([0.18739855]))
    (array([42.36787806,  1.52424563]), array([[  20.37868076, -204.8957381 ],
           [-204.8957381 , 5875.15334198]]), array([50.001]), array([0.18739855]))
    (array([42.36790125,  1.5242471 ]), array([[  20.37868981, -204.89602515],
           [-204.89602515, 5875.16080803]]), array([50.001]), array([0.18739855]))
    (array([42.3678952 ,  1.52424669]), array([[  20.37867365, -204.89568037],
           [-204.89568037, 5875.14979245]]), array([50.001]), array([0.18739855]))
    (array([42.36787708,  1.5242458 ]), array([[  20.37868328, -204.89577668],
           [-204.89577668, 5875.15399516]]), array([50.001]), array([0.18739855]))
    (array([42.3678885 ,  1.52424611]), array([[  20.37868757, -204.89593577],
           [-204.89593577, 5875.15966363]]), array([50.001]), array([0.18739855]))
    (array([42.36787837,  1.52424565]), array([[  20.37868708, -204.89589503],
           [-204.89589503, 5875.15952039]]), array([50.001]), array([0.18739855]))
    (array([42.36790124,  1.5242471 ]), array([[  20.37868956, -204.89602263],
           [-204.89602263, 5875.16073582]]), array([50.001]), array([0.18739855]))
    (array([42.3678952 ,  1.52424669]), array([[  20.37867366, -204.8956804 ],
           [-204.8956804 , 5875.14979346]]), array([50.001]), array([0.18739855]))
    (array([42.36787714,  1.52424581]), array([[  20.37868328, -204.89577668],
           [-204.89577668, 5875.15399505]]), array([50.001]), array([0.18739855]))
    (array([42.36788829,  1.52424609]), array([[  20.37868757, -204.89593577],
           [-204.89593577, 5875.15966367]]), array([50.001]), array([0.18739855]))
    (array([42.36787801,  1.52424563]), array([[  20.37868708, -204.89589505],
           [-204.89589505, 5875.15952098]]), array([50.001]), array([0.18739855]))
    (array([42.36790126,  1.5242471 ]), array([[  20.37868955, -204.89602254],
           [-204.89602254, 5875.16073317]]), array([50.001]), array([0.18739855]))
    (array([42.3678952 ,  1.52424669]), array([[  20.37867364, -204.8956803 ],
           [-204.8956803 , 5875.14979046]]), array([50.001]), array([0.18739855]))
    (array([42.36787708,  1.5242458 ]), array([[  20.37868328, -204.89577668],
           [-204.89577668, 5875.15399522]]), array([50.001]), array([0.18739855]))
    (array([42.3678885 ,  1.52424611]), array([[  20.37868757, -204.89593577],
           [-204.89593577, 5875.15966363]]), array([50.001]), array([0.18739855]))
    (array([42.36787837,  1.52424565]), array([[  20.37868708, -204.89589503],
           [-204.89589503, 5875.15952039]]), array([50.001]), array([0.18739855]))
    (array([42.36790124,  1.5242471 ]), array([[  20.37868956, -204.89602263],
           [-204.89602263, 5875.16073582]]), array([50.001]), array([0.18739855]))
    (array([42.3678952 ,  1.52424669]), array([[  20.37867366, -204.8956804 ],
           [-204.8956804 , 5875.14979346]]), array([50.001]), array([0.18739855]))
    (array([42.36787714,  1.52424581]), array([[  20.37868328, -204.89577668],
           [-204.89577668, 5875.15399505]]), array([50.001]), array([0.18739855]))
    (array([42.36788829,  1.52424609]), array([[  20.37868757, -204.89593577],
           [-204.89593577, 5875.15966367]]), array([50.001]), array([0.18739855]))
    (array([42.36787801,  1.52424563]), array([[  20.37868708, -204.89589505],
           [-204.89589505, 5875.15952098]]), array([50.001]), array([0.18739855]))
    (array([42.36790126,  1.5242471 ]), array([[  20.37868955, -204.89602254],
           [-204.89602254, 5875.16073317]]), array([50.001]), array([0.18739855]))
    (array([42.3678952 ,  1.52424669]), array([[  20.37867364, -204.8956803 ],
           [-204.8956803 , 5875.14979046]]), array([50.001]), array([0.18739855]))
    (array([42.36787708,  1.5242458 ]), array([[  20.37868328, -204.89577668],
           [-204.89577668, 5875.15399522]]), array([50.001]), array([0.18739855]))
    (array([42.3678885 ,  1.52424611]), array([[  20.37868757, -204.89593577],
           [-204.89593577, 5875.15966363]]), array([50.001]), array([0.18739855]))
    (array([42.36787837,  1.52424565]), array([[  20.37868708, -204.89589503],
           [-204.89589503, 5875.15952039]]), array([50.001]), array([0.18739855]))
    (array([42.36790124,  1.5242471 ]), array([[  20.37868956, -204.89602263],
           [-204.89602263, 5875.16073582]]), array([50.001]), array([0.18739855]))
    (array([42.3678952 ,  1.52424669]), array([[  20.37867366, -204.8956804 ],
           [-204.8956804 , 5875.14979346]]), array([50.001]), array([0.18739855]))
    (array([42.36787714,  1.52424581]), array([[  20.37868328, -204.89577668],
           [-204.89577668, 5875.15399505]]), array([50.001]), array([0.18739855]))
    (array([42.36788829,  1.52424609]), array([[  20.37868757, -204.89593577],
           [-204.89593577, 5875.15966367]]), array([50.001]), array([0.18739855]))
    (array([42.36787801,  1.52424563]), array([[  20.37868708, -204.89589505],
           [-204.89589505, 5875.15952098]]), array([50.001]), array([0.18739855]))
    (array([42.36790126,  1.5242471 ]), array([[  20.37868955, -204.89602254],
           [-204.89602254, 5875.16073317]]), array([50.001]), array([0.18739855]))
    (array([42.3678952 ,  1.52424669]), array([[  20.37867364, -204.8956803 ],
           [-204.8956803 , 5875.14979046]]), array([50.001]), array([0.18739855]))
    (array([42.36787708,  1.5242458 ]), array([[  20.37868328, -204.89577668],
           [-204.89577668, 5875.15399522]]), array([50.001]), array([0.18739855]))
    (array([42.3678885 ,  1.52424611]), array([[  20.37868757, -204.89593577],
           [-204.89593577, 5875.15966363]]), array([50.001]), array([0.18739855]))
    (array([42.36787837,  1.52424565]), array([[  20.37868708, -204.89589503],
           [-204.89589503, 5875.15952039]]), array([50.001]), array([0.18739855]))
    (array([42.36790124,  1.5242471 ]), array([[  20.37868956, -204.89602263],
           [-204.89602263, 5875.16073582]]), array([50.001]), array([0.18739855]))
    (array([42.3678952 ,  1.52424669]), array([[  20.37867366, -204.8956804 ],
           [-204.8956804 , 5875.14979346]]), array([50.001]), array([0.18739855]))
    (array([42.36787714,  1.52424581]), array([[  20.37868328, -204.89577668],
           [-204.89577668, 5875.15399505]]), array([50.001]), array([0.18739855]))
    (array([42.36788829,  1.52424609]), array([[  20.37868757, -204.89593577],
           [-204.89593577, 5875.15966367]]), array([50.001]), array([0.18739855]))
    (array([42.36787801,  1.52424563]), array([[  20.37868708, -204.89589505],
           [-204.89589505, 5875.15952098]]), array([50.001]), array([0.18739855]))
    (array([42.36790126,  1.5242471 ]), array([[  20.37868955, -204.89602254],
           [-204.89602254, 5875.16073317]]), array([50.001]), array([0.18739855]))
    (array([42.3678952 ,  1.52424669]), array([[  20.37867364, -204.8956803 ],
           [-204.8956803 , 5875.14979046]]), array([50.001]), array([0.18739855]))
    (array([42.36787708,  1.5242458 ]), array([[  20.37868328, -204.89577668],
           [-204.89577668, 5875.15399522]]), array([50.001]), array([0.18739855]))
    (array([42.3678885 ,  1.52424611]), array([[  20.37868757, -204.89593577],
           [-204.89593577, 5875.15966363]]), array([50.001]), array([0.18739855]))
    (array([42.36787837,  1.52424565]), array([[  20.37868708, -204.89589503],
           [-204.89589503, 5875.15952039]]), array([50.001]), array([0.18739855]))
    (array([42.36790124,  1.5242471 ]), array([[  20.37868956, -204.89602263],
           [-204.89602263, 5875.16073582]]), array([50.001]), array([0.18739855]))
    (array([42.3678952 ,  1.52424669]), array([[  20.37867366, -204.8956804 ],
           [-204.8956804 , 5875.14979346]]), array([50.001]), array([0.18739855]))
    (array([42.36787714,  1.52424581]), array([[  20.37868328, -204.89577668],
           [-204.89577668, 5875.15399505]]), array([50.001]), array([0.18739855]))
    (array([42.36788829,  1.52424609]), array([[  20.37868757, -204.89593577],
           [-204.89593577, 5875.15966367]]), array([50.001]), array([0.18739855]))
    (array([42.36787801,  1.52424563]), array([[  20.37868708, -204.89589505],
           [-204.89589505, 5875.15952098]]), array([50.001]), array([0.18739855]))
    (array([42.36790126,  1.5242471 ]), array([[  20.37868955, -204.89602254],
           [-204.89602254, 5875.16073317]]), array([50.001]), array([0.18739855]))
    (array([42.3678952 ,  1.52424669]), array([[  20.37867364, -204.8956803 ],
           [-204.8956803 , 5875.14979046]]), array([50.001]), array([0.18739855]))
    (array([42.36787708,  1.5242458 ]), array([[  20.37868328, -204.89577668],
           [-204.89577668, 5875.15399522]]), array([50.001]), array([0.18739855]))
    (array([42.3678885 ,  1.52424611]), array([[  20.37868757, -204.89593577],
           [-204.89593577, 5875.15966363]]), array([50.001]), array([0.18739855]))
    (array([42.36787837,  1.52424565]), array([[  20.37868708, -204.89589503],
           [-204.89589503, 5875.15952039]]), array([50.001]), array([0.18739855]))
    (array([42.36790124,  1.5242471 ]), array([[  20.37868956, -204.89602263],
           [-204.89602263, 5875.16073582]]), array([50.001]), array([0.18739855]))
    (array([42.3678952 ,  1.52424669]), array([[  20.37867366, -204.8956804 ],
           [-204.8956804 , 5875.14979346]]), array([50.001]), array([0.18739855]))
    (array([42.36787714,  1.52424581]), array([[  20.37868328, -204.89577668],
           [-204.89577668, 5875.15399505]]), array([50.001]), array([0.18739855]))
    (array([42.36788829,  1.52424609]), array([[  20.37868757, -204.89593577],
           [-204.89593577, 5875.15966367]]), array([50.001]), array([0.18739855]))
    (array([42.36787801,  1.52424563]), array([[  20.37868708, -204.89589505],
           [-204.89589505, 5875.15952098]]), array([50.001]), array([0.18739855]))
    (array([42.36790126,  1.5242471 ]), array([[  20.37868955, -204.89602254],
           [-204.89602254, 5875.16073317]]), array([50.001]), array([0.18739855]))
    (array([42.3678952 ,  1.52424669]), array([[  20.37867364, -204.8956803 ],
           [-204.8956803 , 5875.14979046]]), array([50.001]), array([0.18739855]))
    (array([42.36787708,  1.5242458 ]), array([[  20.37868328, -204.89577668],
           [-204.89577668, 5875.15399522]]), array([50.001]), array([0.18739855]))
    (array([42.3678885 ,  1.52424611]), array([[  20.37868757, -204.89593577],
           [-204.89593577, 5875.15966363]]), array([50.001]), array([0.18739855]))
    (array([42.36787837,  1.52424565]), array([[  20.37868708, -204.89589503],
           [-204.89589503, 5875.15952039]]), array([50.001]), array([0.18739855]))
    (array([42.36790124,  1.5242471 ]), array([[  20.37868956, -204.89602263],
           [-204.89602263, 5875.16073582]]), array([50.001]), array([0.18739855]))
    (array([42.3678952 ,  1.52424669]), array([[  20.37867366, -204.8956804 ],
           [-204.8956804 , 5875.14979346]]), array([50.001]), array([0.18739855]))
    (array([42.36787714,  1.52424581]), array([[  20.37868328, -204.89577668],
           [-204.89577668, 5875.15399505]]), array([50.001]), array([0.18739855]))
    (array([42.36788829,  1.52424609]), array([[  20.37868757, -204.89593577],
           [-204.89593577, 5875.15966367]]), array([50.001]), array([0.18739855]))
    (array([42.36787801,  1.52424563]), array([[  20.37868708, -204.89589505],
           [-204.89589505, 5875.15952098]]), array([50.001]), array([0.18739855]))
    (array([42.36790126,  1.5242471 ]), array([[  20.37868955, -204.89602254],
           [-204.89602254, 5875.16073317]]), array([50.001]), array([0.18739855]))



```python
estimate = np.zeros((100,5))
for i in range(5): 
    estimate[:,i] = f(theta_n[i*20], t)
```


```python
plt.plot(t, model, linewidth = 3, alpha = 0.6)
plt.scatter(t, data, c = 'orange')
plt.show()

plt.plot(t, estimate, linewidth = 1, alpha = 0.6)
plt.scatter(t, data, c = 'orange')
plt.show()
```

7902de7c05e5485fade5c160fa2453f9b9a83306
ipynb
Jupyter Notebook
python_tutorials/Variational_Inference.ipynb
PavanChaggar/Bayesian_inference
2a00d8f9fc391c6529c001be617a776a6b60c87b
Apache-2.0
2020-07-22T18:47:33.000Z
2021-07-26T12:06:50.000Z
python_tutorials/Variational_Inference.ipynb
PavanChaggar/Bayesian_inference
2a00d8f9fc391c6529c001be617a776a6b60c87b
Apache-2.0
python_tutorials/Variational_Inference.ipynb
PavanChaggar/Bayesian_inference
2a00d8f9fc391c6529c001be617a776a6b60c87b
Apache-2.0
Qwen/Qwen-72B
1. YES
2. YES