<jupyter_start><jupyter_text>Observation Error DemonstrationThis notebook will show a simple example to show the incorporation ofobservation error, and the effect of the associated parameter thatscales the observation error.This notebook continues with the example in the linear_response_example.ipynbnotebook; a simple linear model makes validating the result simple.Read in observations and modeling functions. The main function herefacilitates setting and fixing various model parameters. In this case,fixing the parameters that correspond to a linear emulator response.<jupyter_code>import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import seaborn as sns

from sepia import SepiaData, SepiaModel, SepiaPlot, \
    SepiaEmulatorPrediction, SepiaFullPrediction

def run_model(m=4, fix_lamOs=None, use_Sigy=None):
    t_sim=np.linspace(1,4,m).reshape((-1,1))
    y_sim=np.linspace(4,1,m).reshape((-1,1))
    y_obs=np.ones((1,1))*2.5

    # Use the observation error if supplied
    if not use_Sigy:
        dat = SepiaData(t_sim=t_sim, y_sim=y_sim, y_obs=y_obs)
    else:
        dat = SepiaData(t_sim=t_sim, y_sim=y_sim, y_obs=y_obs, Sigy=use_Sigy)

    dat.standardize_y()
    dat.transform_xt()

    mod=SepiaModel(dat)

    # Change default sepia model parameter based
    # on flags/values passed in.
    if fix_lamOs:
        # Scaling of observation error
        mod.set_param('lamOs',fix=fix_lamOs)

    #fix other parameters to allow a highly-certain linear response emulator
    # Process precision - so it's not forcing N(0,1) marginal response locally
    mod.set_param('lamUz',fix=0.01)
    # Nugget / noise process variance for simulations = near-zero variance
    mod.set_param('lamWs',fix=1e6)
    # Scaling of simulator error = near-zero variance
    # (which is relevant to allocating relative uncertainty for multivariate response)
    mod.set_param('lamWOs',fix=1e6)

    # "fit" the model by generating posterior samples
    mod.do_mcmc(100)  # go past transient
    mod.clear_samples()
    mod.do_mcmc(5000) # collect samples

    return mod<jupyter_output><empty_output><jupyter_text>And some plotting routines<jupyter_code>def plot_trace(mod):
    #call the trace plot built-in with the samples
    samp=mod.get_samples()
    SepiaPlot.mcmc_trace(samp)
    plt.show()

def plot_cal(mod):
    # plot the calibration setup and results
    psamp=mod.get_samples(numsamples=500)
    t_pred=np.linspace(0,5,50).reshape((-1,1))
    spred=SepiaEmulatorPrediction(t_pred=t_pred, model=mod, samples=psamp)
    ysp=spred.get_y()
    yspm=np.mean(ysp,axis=0).squeeze()
    ysps=np.std(ysp,axis=0).squeeze()

    opred=SepiaFullPrediction(model=mod,samples=psamp)
    yop=opred.get_yobs()

    fig,ax1 = plt.subplots()
    ax1.plot(t_pred,yspm,'k')
    ax1.plot(t_pred,yspm+2*ysps,'k.')
    ax1.plot(t_pred,yspm-2*ysps,'k.')
    ax1.plot(mod.data.sim_data.t,mod.data.sim_data.y,'go')
    ax1.plot(psamp['theta_native'],yop.squeeze(),'r.')
    calx=np.array([0,5]).reshape((-1,1))
    ax1.plot(calx,np.ones((2,1))*mod.data.obs_data.y,'b')
    ax1.plot(calx,np.ones((2,1))*(mod.data.obs_data.y+2*np.sqrt(mod.data.obs_data.Sigy)),'b:')
    ax1.plot(calx,np.ones((2,1))*(mod.data.obs_data.y-2*np.sqrt(mod.data.obs_data.Sigy)),'b:')
    ax2=ax1.twinx()
    sns.histplot(psamp['theta_native'],kde=True,color='r',ax=ax2,legend=False)
    ax2.set_ylabel(''); ax2.set_yticks([]); ax2.set_ybound([i*5 for i in ax2.get_ybound()])
    plt.show()

    samp=mod.get_samples()
    theta_std=np.std(samp['theta_native'],ddof=1)
    print('Theta std = %f, var = %f'%(theta_std,theta_std**2))<jupyter_output><empty_output><jupyter_text>See the linearResponseExample for more detail on the settings.Sepia handles emulation of multivariate responses, but we will notshow that in this example. The Al_5083_calibration example python fileand notebook shows a more full model.The observation error is a covariance matrix, Sigy. Often Sigy willonly have entries on this diagonal of this covariance matrix, indicatingerrors are considered independent. In this simple illustrative examplethere is a scalar response, so Sigy is a scalar value.Generate a model with an observation error of N(0,0.1), and calibrate theparameter:<jupyter_code>mod=run_model(fix_lamOs=1, use_Sigy=0.1)<jupyter_output>MCMC sampling: 100%|██████████| 100/100 [00:00<00:00, 788.61it/s]
MCMC sampling: 100%|██████████| 5000/5000 [00:05<00:00, 836.18it/s]<jupyter_text>"Trace plot" shows model parameter samples. This is not very exciting, since themany of the model parameters were set to a fixed value to make thisexample's function very clear.<jupyter_code>plot_trace(mod)<jupyter_output><empty_output><jupyter_text>The key observation to make about the parameters is that we fixedlamOs to 1. lamOS is a precision that divides the Sigy values. More onthat below, but by fixing this at 1, the Sigy will be used exactly.Plot of the emulator and calibration results:<jupyter_code>plot_cal(mod)<jupyter_output><empty_output><jupyter_text>The green points are 'simulations' used to build the emulator. The emulatorpredictions across the domain are in black, and the uncertainty of theemulator is plotted but not visible. There is one observation, at a value2.5 on the y axis, projected across on the blue line. The dotted blue lineis the +/- 2sd values of the observation's uncertainty. The red points aresamples of the calibrated parameter value. This example uses a linear emulatorwith a $45^\circ$ slope so the distribution of the calibrated variable shouldbe the same distribution as the observation. The distribution of theposterior (calibrated) parameter samples is shown in the inset histogram.The standard deviation and variance of these parameter samples is printed.The variance is (approximately, to sampling uncertainty) 0.1, as expected.Another example, using variance 0.25:<jupyter_code>mod=run_model(fix_lamOs=1, use_Sigy=0.25)
plot_cal(mod)<jupyter_output>MCMC sampling: 100%|██████████| 100/100 [00:00<00:00, 582.23it/s]
MCMC sampling: 100%|██████████| 5000/5000 [00:06<00:00, 723.79it/s]<jupyter_text>Finally, some discussion of the parameter lamOs, which is fixed at 1 inthe models above. Often in modeling the observation error is not a complete accounting of alluncertainty, althoughit is expected to be a relative value for multivariate observations. Giventhe data, the observation error needs to be "inflated" to allow the datato be consistent. And on the other hand, if the data is consistent internallyand with the model, we can be satisfied with a posterior result for parametercalibration that corresponds to a smaller observationvariance. The precision lamOs can be fit to allow the observationvariance to be scaled given the rest of the model. This is function of thelamOs parameter, and the default setup of a sepia model has priors on lamOsthat allow these cases. This default expects that the lamOs, and observationerror, will be adjusted given data.In the example below, we continue with an assumption of linearity for theemulator, but will set the observation error variance to be very large, in thiscase the default variance of 1 for a error distribution N(0,1).The lamOs parameter dividing that observation variance will be calibrated.<jupyter_code>mod=run_model()
plot_cal(mod)<jupyter_output>MCMC sampling: 100%|██████████| 100/100 [00:00<00:00, 674.55it/s]
MCMC sampling: 100%|██████████| 5000/5000 [00:07<00:00, 705.31it/s]<jupyter_text>The result is that Sepia's default priors cause lamOs to be large, reducingthe observation error and the corresponding posterior calibrated distributionof the parameter. The lamOs posterior sampled is much larger than 1:<jupyter_code>print('lamOs sample mean = %f'%(np.mean(mod.get_samples()['lamOs'])))<jupyter_output>lamOs sample mean = 940.664634
{'fraction_non_alphanumeric': 0.10071758072783188, 'fraction_numerical': 0.032162993336750385, 'mean_word_length': 5.254006410256411, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 29, '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': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '522529', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2724, 'n_tokens_neox': 2482, 'n_words': 793, 'path': 'examples/observation_error_example.ipynb', 'repo_name': 'lanl/SEPIA', 'size': '7804'}