<jupyter_start><jupyter_text>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.<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>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)
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)
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]))
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()<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>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)
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)
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.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()<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>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)

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)
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]))
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()<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>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

import scipy.spatial.distance
dist = scipy.spatial.distance.squareform(scipy.spatial.distance.pdist(np.expand_dims(xdata, 1)))
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)
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]))
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()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14533506985555292, 'fraction_numerical': 0.03528297418896519, 'mean_word_length': 3.666574585635359, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '204253', 'license': 'MIT', 'n_tokens_mistral': 6387, 'n_tokens_neox': 5975, 'n_words': 1656, 'path': 'assets/notebooks/Untitled.ipynb', 'repo_name': 'mfkasim1/mfkasim1.github.io', 'size': '16892'}