```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
import emcee
```

# Fitting with model uncertainty

In fitting spectrum routine, given an observed spectrum and a computational model to simulate the spectrum, we are interested in determining the parameters (with its error) that gives the most similar spectrum to the observed one.

If we know the "true" underlying model of the spectrum (i.e. we can simulate the spectrum **exactly**, given a set of parameters), obtaining the error of the parameters are just like a linear fit. However, if we only have an approximate model, the error can be severely underdetermined.

In this notebook, I will show the problem of fitting with an approximate model and how to solve it (hopefully).

## The true model and the true data

In this notebook, we consider a simplified model of generating a spectrum. Let's say we have 4 parameters in which determining the values at \\(x = (0, 0.3, 0.7, 1.0)\\). All the other values are obtained using the cubic spline interpolation. We call this the "**cubic**" model which is the true model. Suppose in an experiment, we can observe 1500 data points linearly spaced from \\([0, 1.0]\\), with some noise.


```python
np.random.seed(123)
ndata = 1500
xparameters = [0, 0.3, 0.7, 1.0]
parameters = [0.1, 0.7, 0.2, 0.4]
sigma_noise = 0.05

finterp1 = interp1d(xparameters, parameters, kind="cubic")
xdata = np.linspace(0, 1, ndata)
ydata = finterp1(xdata) + np.random.randn(*xdata.shape) * sigma_noise

plt.plot(xdata, ydata, '.')
plt.show()
```

This is the data we observe (as an example).

## Using the true model and known noise level

Assume that we are in an ideal world where we know the true model and we also know the noise level. The Bayesian inference of the parameters, $\mathbf{p}$, knowing the data, $\mathcal{D}$, can be written as

$$\begin{equation}
P(\mathbf{p} | \mathcal{D}) \propto P(\mathcal{D} | \mathbf{p}) P(\mathbf{p}).
\end{equation}$$

Take a uniform prior for the parameters and the likelihood is just a Gaussian distribution. Therefore,

$$\begin{equation}
P(\mathbf{p} | \mathcal{D}) \propto \exp\left[-\sum_i\frac{\left(\mathbf{\hat{y_i}} - \mathbf{y_i(\mathbf{p})} \right)^2}{2\sigma^2}\right].
\end{equation}$$

We can use the `emcee` module to draw samples.



```python
def get_ymodel1(xparams, params, xdata):
    finterp1 = interp1d(xparams, params, kind="cubic")
    return finterp1(xdata)

def lnprob1(params, xparams, xdata, ydata):
    ymodel = get_ymodel1(xparams, params, xdata)
    sigma = 0.05 # known noise level
    return -np.sum((ydata - ymodel)**2) / (2 * sigma**2)
```


```python
ndim, nwalkers = 4, 100
p0 = np.random.random((nwalkers, ndim))
sampler1 = emcee.EnsembleSampler(nwalkers, ndim, lnprob1, args=[xparameters, xdata, ydata])
_ = sampler1.run_mcmc(p0, 1000)
```


```python
nburn = 10000
params1 = sampler1.flatchain[nburn:,:]
plt.figure(figsize=(12,4))
for i in range(ndim):
    plt.subplot(1, ndim, i+1)
    ax = plt.gca()
    plt.hist(params1[:,i], 50, normed=True, range=parameters[i] + np.array([-1, 1])*0.05)
    ylim = ax.get_ylim()
    plt.plot(np.zeros(2) + parameters[i], ylim, 'k--')
    ax.set_ylim(ylim)
    plt.title("parameter #%d" % (i+1))

plt.show()

for i in range(ndim):
    print("parameter #%d: %f +- %f (true: %f)" % \
         (i+1, np.median(params1[:,i]),
          np.percentile(params1[:,i], 84) - np.percentile(params1[:,i], 16),
          parameters[i]))
```


```python
plt.figure(figsize=(8,4))
plt.plot(xdata, ydata, '.')
plt.plot(xdata, get_ymodel1(xparameters, np.median(params1, 0), xdata), '-', linewidth=2)

# unseen on the graph because it is too narrow
plt.fill_between(xdata, 
                 get_ymodel1(xparameters, np.percentile(params1, 16, axis=0), xdata),
                 get_ymodel1(xparameters, np.percentile(params1, 84, axis=0), xdata), facecolor="orange")
plt.show()
```

As we can see, the retrieved parameters (i.e. the median) are very close to the true values.

## Using an approximate model and known noise level

Now we are going to a slightly less ideal world where we used an approximate model with known noise level. The approximate model here uses the quadratic and linear spline interpolation (the "**quadratic**" and "**linear**" model) instead of the cubic spline interpolation (the "**cubic**" model).
The inference is still the same as before, we just change how we calculate the model data in the `lnprob1` function.


```python
def get_ymodel2(xparams, params, xdata):
    finterp1 = interp1d(xparams, params, kind="linear")
    return finterp1(xdata)

def get_ymodel3(xparams, params, xdata):
    finterp1 = interp1d(xparams, params, kind="quadratic")
    return finterp1(xdata)

def lnprob2(params, xparams, xdata, ydata):
    ymodel = get_ymodel2(xparams, params, xdata)
    sigma = 0.05 # known noise level
    return -np.sum((ydata - ymodel)**2) / (2 * sigma**2)

def lnprob3(params, xparams, xdata, ydata):
    ymodel = get_ymodel3(xparams, params, xdata)
    sigma = 0.05 # known noise level
    return -np.sum((ydata - ymodel)**2) / (2 * sigma**2)
```


```python
ndim, nwalkers = 4, 100
p0 = np.random.random((nwalkers, ndim))
sampler2 = emcee.EnsembleSampler(nwalkers, ndim, lnprob2, args=[xparameters, xdata, ydata])
_ = sampler2.run_mcmc(p0, 1000)
sampler3 = emcee.EnsembleSampler(nwalkers, ndim, lnprob3, args=[xparameters, xdata, ydata])
_ = sampler3.run_mcmc(p0, 1000)
```


```python
nburn = 10000
params2 = sampler2.flatchain[nburn:,:]
params3 = sampler3.flatchain[nburn:,:]
plt.figure(figsize=(16,6))
for i in range(ndim):
    for j in range(2):
        paramss = params2 if j == 0 else params3
        
        plt.subplot(2, ndim, j*ndim+i+1)
        ax = plt.gca()
        plt.hist(paramss[:,i], 50, normed=True, range=parameters[i] + np.array([-1, 1])*0.2)
        ylim = ax.get_ylim()
        plt.plot(np.zeros(2) + parameters[i], ylim, 'k--')
        ax.set_ylim(ylim)
        if j == 0: 
            plt.title("parameter #%d" % (i+1))
        if i == 0:
            plt.ylabel("%s model" % ("linear" if j == 0 else "quadratic"))

plt.show()

for j in range(2):
    print("%s model" % ("linear" if j == 0 else "quadratic"))
    paramss = params2 if j == 0 else params3
    for i in range(ndim):
        print("parameter #%d: %f +- %f (true: %f)" % \
             (i+1, np.median(paramss[:,i]),
              np.percentile(paramss[:,i], 84) - np.percentile(paramss[:,i], 16),
              parameters[i]))
```


```python
plt.figure(figsize=(16,4))
plt.subplot(1,2,1)
plt.plot(xdata, ydata, '.')
plt.plot(xdata, get_ymodel2(xparameters, np.median(params2, 0), xdata), '-', linewidth=2)
# unseen on the graph because it is too narrow
plt.fill_between(xdata, 
                 get_ymodel2(xparameters, np.percentile(params2, 16, axis=0), xdata),
                 get_ymodel2(xparameters, np.percentile(params2, 84, axis=0), xdata), facecolor="orange")
plt.subplot(1,2,2)
plt.plot(xdata, ydata, '.')
plt.plot(xdata, get_ymodel3(xparameters, np.median(params3, 0), xdata), '-', linewidth=2)
# unseen on the graph because it is too narrow
plt.fill_between(xdata, 
                 get_ymodel3(xparameters, np.percentile(params3, 16, axis=0), xdata),
                 get_ymodel3(xparameters, np.percentile(params3, 84, axis=0), xdata), facecolor="orange")
plt.show()
```

As we can see here, if we are using approximate models, the retrieved values of the parameteres are off while the produced error is still small. In other words, it is confident, but it is wrong. If we are using less accurate model, the effect is even worse: the predicted values are off by much and the error remains small.

The problem with the inference we used in this case is that we only assume that the data differ only by the Gaussian noise. This is not true since the data also differ due to the model inadequacy.

#### Max deviation inadequacy

In my latest paper, I assume the model inadequacy as stochastic variables with distribution

$$\begin{equation}
P(\eta(x)) \propto \exp\left[-\frac{\max\left(\mathbf{\hat{y_i}} - \mathbf{y_i}\right)^2}{2d^2}\right]
\end{equation}$$

with value of $d$ approximately the values of the noise.

Unfortunately, I didn't take into account the effect of the Gaussian noise. So the distribution above becomes the likelihood of observing the data.


```python
def lnprob2_model1(params, xparams, xdata, ydata):
    ymodel = get_ymodel2(xparams, params, xdata)
    d = 0.05 # d approximately the noise level
    return -np.max((ydata - ymodel)**2) / (2 * d**2)

def lnprob3_model1(params, xparams, xdata, ydata):
    ymodel = get_ymodel3(xparams, params, xdata)
    d = 0.05 # d approximately the noise level
    return -np.max((ydata - ymodel)**2) / (2 * d**2)

```


```python
ndim, nwalkers = 4, 100
p0 = np.random.random((nwalkers, ndim))
sampler2_model1 = emcee.EnsembleSampler(nwalkers, ndim, lnprob2_model1, args=[xparameters, xdata, ydata])
_ = sampler2_model1.run_mcmc(p0, 1000)
sampler3_model1 = emcee.EnsembleSampler(nwalkers, ndim, lnprob3_model1, args=[xparameters, xdata, ydata])
_ = sampler3_model1.run_mcmc(p0, 1000)
```


```python
nburn = 10000
params2_model1 = sampler2_model1.flatchain[nburn:,:]
params3_model1 = sampler3_model1.flatchain[nburn:,:]
plt.figure(figsize=(16,6))
for i in range(ndim):
    for j in range(2):
        paramss = params2_model1 if j == 0 else params3_model1
        
        plt.subplot(2, ndim, j*ndim+i+1)
        ax = plt.gca()
        plt.hist(paramss[:,i], 50, normed=True, range=parameters[i] + np.array([-1, 1])*0.2)
        ylim = ax.get_ylim()
        plt.plot(np.zeros(2) + parameters[i], ylim, 'k--')
        ax.set_ylim(ylim)
        if j == 0: 
            plt.title("parameter #%d" % (i+1))
        if i == 0:
            plt.ylabel("%s model" % ("linear" if j == 0 else "quadratic"))

plt.show()

for j in range(2):
    print("%s model" % ("linear" if j == 0 else "quadratic"))
    paramss = params2_model1 if j == 0 else params3_model1
    for i in range(ndim):
        print("parameter #%d: %f +- %f (true: %f)" % \
             (i+1, np.median(paramss[:,i]),
              np.percentile(paramss[:,i], 84) - np.percentile(paramss[:,i], 16),
              parameters[i]))
```


```python
plt.figure(figsize=(16,4))
plt.subplot(1,2,1)
plt.plot(xdata, ydata, '.', alpha=0.3)
plt.plot(xdata, get_ymodel2(xparameters, np.median(params2_model1, 0), xdata), '-', linewidth=2)
# unseen on the graph because it is too narrow
plt.fill_between(xdata, 
                 get_ymodel2(xparameters, np.percentile(params2_model1, 16, axis=0), xdata),
                 get_ymodel2(xparameters, np.percentile(params2_model1, 84, axis=0), xdata),
                 facecolor="orange", alpha=0.5)
plt.subplot(1,2,2)
plt.plot(xdata, ydata, '.', alpha=0.3)
plt.plot(xdata, get_ymodel3(xparameters, np.median(params3_model1, 0), xdata), '-', linewidth=2)
# unseen on the graph because it is too narrow
plt.fill_between(xdata, 
                 get_ymodel3(xparameters, np.percentile(params3_model1, 16, axis=0), xdata),
                 get_ymodel3(xparameters, np.percentile(params3_model1, 84, axis=0), xdata),
                 facecolor="orange", alpha=0.5)
plt.plot(xdata, get_ymodel1(xparameters, parameters, xdata), '-', linewidth=2)

nburn = 10000
params2 = sampler2.flatchain[nburn:,:]
params3 = sampler3.flatchain[nburn:,:]
plt.figure(figsize=(16,6))
for i in range(ndim):
    for j in range(2):
        paramss = params2 if j == 0 else params3
        
        plt.subplot(2, ndim, j*ndim+i+1)
        ax = plt.gca()
        plt.hist(paramss[:,i], 50, normed=True, range=parameters[i] + np.array([-1, 1])*0.2)
        ylim = ax.get_ylim()
        plt.plot(np.zeros(2) + parameters[i], ylim, 'k--')
        ax.set_ylim(ylim)
        if j == 0: 
            plt.title("parameter #%d" % (i+1))
        if i == 0:
            plt.ylabel("%s model" % ("linear" if j == 0 else "quadratic"))

plt.show()

for j in range(2):
    print("%s model" % ("linear" if j == 0 else "quadratic"))
    paramss = params2 if j == 0 else params3
    for i in range(ndim):
        print("parameter #%d: %f +- %f (true: %f)" % \
             (i+1, np.median(paramss[:,i]),
              np.percentile(paramss[:,i], 84) - np.percentile(paramss[:,i], 16),
              parameters[i]))
plt.show()
```

From the results above, this model of the model inadequacy works well if the model is approximately accurate. If the model is quite far from the true model, then the result is quite far off (compared to its standard deviation). Moreover, this model requires us to make an assumption on how *right* (or *wrong*) the model is by setting the value of $d$ manually.

What we want is minimal input from users and the model is automatically adjust the error on the retrieved parameters. If the model is less accurate, then the error on the retrieved parameters should be larger. If the model is more accurate, then the error should be smaller.

#### Gaussian process model inadequacy

Another way to model the model inadequacy is using Gaussian Process.


```python
from scipy.stats import multivariate_normal

def lnprob2_model2(params, xparams, xdata, ydata, dist):
    d = params[-1]
    m = params[-2]
    sigma = 0.05
    if d > 1 or d < 0 or m < 0 or sigma < 0: return -np.inf
    
    ncut = 10
    cov = m*m * np.exp(-dist[::ncut,::ncut]**2/(2*d*d)) + np.eye(len(xdata[::ncut])) * (sigma*sigma + 1e-7)
    
    ymodel = get_ymodel2(xparams, params[:-3], xdata)
    obs = ymodel - ydata
    obs = obs[::ncut]
    try:
        return multivariate_normal.logpdf(obs, mean=np.zeros(obs.shape[0]), cov=cov) - np.log(sigma * m * d)
    except:
        return -np.inf

def lnprob3_model2(params, xparams, xdata, ydata, dist):
    d = params[-1]
    m = params[-2]
    sigma = 0.05
    if d > 1 or d < 0 or m < 0 or sigma < 0: return -np.inf
    
    ncut = 10
    cov = m*m * np.exp(-dist[::ncut,::ncut]**2/(2*d*d)) + np.eye(len(xdata[::ncut])) * (sigma*sigma + 1e-7)
    
    ymodel = get_ymodel3(xparams, params[:-3], xdata)
    obs = ymodel - ydata
    obs = obs[::ncut]
    try:
        return multivariate_normal.logpdf(obs, mean=np.zeros(obs.shape[0]), cov=cov) - np.log(sigma * m * d)
    except:
        return -np.inf

```


```python
import scipy.spatial.distance
dist = scipy.spatial.distance.squareform(scipy.spatial.distance.pdist(np.expand_dims(xdata, 1)))
```


```python
ndim, nwalkers = 4, 100
nhparams = 3
p0 = np.random.random((nwalkers, ndim+nhparams))
sampler2_model2 = emcee.EnsembleSampler(nwalkers, ndim+nhparams, lnprob2_model2, args=[xparameters, xdata, ydata, dist])
%time _ = sampler2_model2.run_mcmc(p0, 1000)
sampler3_model2 = emcee.EnsembleSampler(nwalkers, ndim+nhparams, lnprob3_model2, args=[xparameters, xdata, ydata, dist])
%time _ = sampler3_model2.run_mcmc(p0, 1000)
```

    CPU times: user 20min 17s, sys: 17min 33s, total: 37min 51s
    Wall time: 5min 14s
    CPU times: user 21min, sys: 18min 37s, total: 39min 37s
    Wall time: 5min 34s



```python
nburn = 10000
params2_model2 = sampler2_model2.flatchain[nburn:,:]
params3_model2 = sampler3_model2.flatchain[nburn:,:]
plt.figure(figsize=(16,6))
for i in range(ndim):
    for j in range(2):
        paramss = params2_model2 if j == 0 else params3_model2
        
        plt.subplot(2, ndim, j*ndim+i+1)
        ax = plt.gca()
        plt.hist(paramss[:,i], 50, normed=True, range=parameters[i] + np.array([-1, 1])*2)
        ylim = ax.get_ylim()
        plt.plot(np.zeros(2) + parameters[i], ylim, 'k--')
        ax.set_ylim(ylim)
        if j == 0: 
            plt.title("parameter #%d" % (i+1))
        if i == 0:
            plt.ylabel("%s model" % ("linear" if j == 0 else "quadratic"))

plt.show()

for j in range(2):
    print("%s model" % ("linear" if j == 0 else "quadratic"))
    paramss = params2_model2 if j == 0 else params3_model2
    for i in range(ndim):
        print("parameter #%d: %f +- %f (true: %f)" % \
             (i+1, np.median(paramss[:,i]),
              np.percentile(paramss[:,i], 84) - np.percentile(paramss[:,i], 16),
              parameters[i]))
```


```python
plt.figure(figsize=(16,4))
plt.subplot(1,2,1)
plt.plot(xdata, ydata, '.', alpha=0.3)
plt.plot(xdata, get_ymodel2(xparameters, np.median(params2_model2[:,:-nhparams], 0), xdata), '-', linewidth=2)
# unseen on the graph because it is too narrow
plt.fill_between(xdata, 
                 get_ymodel2(xparameters, np.percentile(params2_model2[:,:-nhparams], 16, axis=0), xdata),
                 get_ymodel2(xparameters, np.percentile(params2_model2[:,:-nhparams], 84, axis=0), xdata),
                 facecolor="orange", alpha=0.5)
plt.subplot(1,2,2)
plt.plot(xdata, ydata, '.', alpha=0.3)
plt.plot(xdata, get_ymodel3(xparameters, np.median(params3_model2[:,:-nhparams], 0), xdata), '-', linewidth=2)
# unseen on the graph because it is too narrow
plt.fill_between(xdata, 
                 get_ymodel3(xparameters, np.percentile(params3_model2[:,:-nhparams], 16, axis=0), xdata),
                 get_ymodel3(xparameters, np.percentile(params3_model2[:,:-nhparams], 84, axis=0), xdata),
                 facecolor="orange", alpha=0.5)
plt.plot(xdata, get_ymodel1(xparameters, parameters, xdata), '-', linewidth=2)
plt.show()
```


```python

```

ed8d8f6b03452a226c240d55cee42942efdba888
ipynb
Jupyter Notebook
assets/notebooks/Untitled.ipynb
mfkasim1/mfkasim1.github.io
2755ad0fa2a27e6e790fbe67d3635454ef790e49
MIT
2020-07-22T07:28:39.000Z
2020-07-22T07:28:39.000Z
assets/notebooks/Untitled.ipynb
mfkasim1/mfkasim1.github.io
2755ad0fa2a27e6e790fbe67d3635454ef790e49
MIT
assets/notebooks/Untitled.ipynb
mfkasim1/mfkasim1.github.io
2755ad0fa2a27e6e790fbe67d3635454ef790e49
MIT
2020-07-22T07:28:40.000Z
2020-07-22T07:28:40.000Z
Qwen/Qwen-72B
1. YES
2. YES