KRIGING METAMODELING: NOISY MODEL RESPONSE (GAUSSIAN PROCESS REGRESSION)

This example illustrates how to perform Kriging regression on noisy data. The example is based on a simple one-dimensional function similar to the one in uq_Example_Kriging_01_1D, now with additive Gaussian noise added to the model response. Three cases of noise variance are considered: unknown homogeneous (homoscedastic), known homoscedastic, and known heterogeneous (heteroscedastic).

Package imports

from uqpylab import sessions
import numpy as np 
import plotly
plotly.offline.init_notebook_mode()

Start a remote UQCloud session

# Start the session
mySession = sessions.cloud()
# (Optional) Get a convenient handle to the command line interface
uq = mySession.cli
# Reset the session
mySession.reset()

 uqpylab.sessions :: INFO     :: A new session (f460c39c10764c04893085457bb9e179) started.
 uqpylab.sessions :: INFO     :: Reset successful.

Set the random seed for reproducibility

uq.rng(0,'twister');
np.random.seed(0)

Computational model

The computational model is a simple analytical function defined by:

$$y(x) = x \sin(x), \; x \in [-3\pi, 3\pi]$$

Specify this model using a string and create a MODEL object:

ModelOpts = {
    'Type': 'Model',
    'ModelFun': 'xsinx.model',
    'isVectorized': 'true'
}

myModel = uq.createModel(ModelOpts)

Probabilistic input model

The probabilistic input model consists of a single uniform random variable:

$X \sim \mathcal{U}(-3\pi, 3\pi)$

Specify its marginal distribution and create a INPUT object:

InputOpts = {
    'Marginals': {
        'Type': 'Uniform',
        'Parameters': [-3*np.pi, 3*np.pi]
    }
}

myInput = uq.createInput(InputOpts)

Kriging regression models

Three cases of noise variance are considered below: unknown homogeneous (homoscedastic), known homoscedastic, and known heterogeneous (heteroscedastic).

Unknown homogeneous (homoscedastic) noise

In the first case, a Kriging model is built on a noisy data set, while also estimating the unknown homogeneous (homoscedastic) noise variance.

# Create an experimental design
Ntrain = 50
X = uq.getSample(myInput, Ntrain)

Evaluate the corresponding model responses:

Y = uq.evalModel(myModel, X)

Add random Gaussian noise with $\sigma_\epsilon = 0.2 \times \sigma_Y$ to the model response:

Y = Y + 0.2*np.std(Y)*np.random.randn(Y.shape[0],1)

Select the metamodeling tool and the Kriging module:

MetaOpts = {
    'Type': 'Metamodel',
    'MetaType': 'Kriging'
}

Use the experimental design and corresponding model responses generated earlier:

MetaOpts['ExpDesign'] = {
    'X': X.tolist(),
    'Y': Y.tolist()
}

Estimate the homogeneous noise variance:

MetaOpts['Regression'] = {
    'SigmaNSQ': 'auto'
}

Create the Kriging metamodel:

myKrigingRegression1 = uq.createModel(MetaOpts)

Print out a report on the resulting Kriging object:

uq.print(myKrigingRegression1)

%---------------- Kriging metamodel ----------------%
    Object Name:            Model 2      
    Input Dimension:        1            
    Output Dimension:       1            

    Experimental Design    
       X size:              [50x1]       
       Y size:              [50x1]       
       Sampling:            User         

    Trend                  
       Type:                ordinary     
       Degree:              0            
       Beta:                [-2.69816	]  

    Gaussian Process (GP)  
       Corr. type:          ellipsoidal  
       Corr. isotropy:      anisotropic  
       Corr. family:        matern-5_2   
       sigma^2:             6.68357e+02  
       Estimation method:   Cross-validation (CV)

    Hyperparameters        
       theta:               [ 0.97298	]  
       Optim. method:       Hybrid Genetic Alg.

    GP Regression          
       Mode:                regression   
       Est. noise:          true         
       sigmaN^2:            8.25276e-01  

    Error estimates        
       Leave-one-out:       6.58586e-02  
%---------------------------------------------------%

Visualize the result:

fig = uq.display(myKrigingRegression1)
# Set manually the y axis range for comparison with the other models
fig.update_yaxes(range=[-10, 15])
fig.show()

Known homogeneous (homoscedastic) noise

The second case illustrates how to build a Kriging model for a noisy data assuming the homogeneous noise variance is known a priori.

Create an experimental design:

Ntrain = 50
X = uq.getSample(myInput, Ntrain)

Evaluate the corresponding model responses:

Y = uq.evalModel(myModel, X)

Add a random noise to the model responses:

noiseVar = 1.0
Y = Y + np.sqrt(noiseVar)*np.random.randn(Y.shape[0],1)

Select the metamodeling tool and the Kriging module:

MetaOpts = {
    'Type': 'Metamodel',
    'MetaType': 'Kriging'
}

Use the experimental design and corresponding model responses generated earlier:

MetaOpts['ExpDesign'] = {
    'X': X.tolist(),
    'Y': Y.tolist()
}

Impose a known homoscedastic noise variance:

MetaOpts['Regression'] = {
    'SigmaNSQ': noiseVar
}

Create the Kriging metamodel:

myKrigingRegression2 = uq.createModel(MetaOpts);

Print out a report on the resulting Kriging object:

uq.print(myKrigingRegression2)

%---------------- Kriging metamodel ----------------%
    Object Name:            Model 3      
    Input Dimension:        1            
    Output Dimension:       1            

    Experimental Design    
       X size:              [50x1]       
       Y size:              [50x1]       
       Sampling:            User         

    Trend                  
       Type:                ordinary     
       Degree:              0            
       Beta:                [-5.52271	]  

    Gaussian Process (GP)  
       Corr. type:          ellipsoidal  
       Corr. isotropy:      anisotropic  
       Corr. family:        matern-5_2   
       sigma^2:             9.94108e+02  
       Estimation method:   Cross-validation (CV)

    Hyperparameters        
       theta:               [ 1.14666	]  
       Optim. method:       Hybrid Genetic Alg.

    GP Regression          
       Mode:                regression   
       Est. noise:          false        
       sigmaN^2:            1            

    Error estimates        
       Leave-one-out:       7.58969e-02  
%---------------------------------------------------%

Plot a representation of the mean and the 95% confidence bounds of the Kriging regression model:

fig = uq.display(myKrigingRegression2)
fig.update_yaxes(range=[-10, 15])
fig.show()

Known non-homogeneous (heteroscedastic) noise

In the third case, Kriging regression model is built from a data set with non-homogeneous (heteroscedastic) noise variance, known at individual data points. The variances differ but remain independent.

Create an experimental design:

Ntrain = 50
X = uq.getSample(myInput, Ntrain, 'grid')

Evaluate the corresponding model responses:

Y = uq.evalModel(myModel, X)

Define a vector of noise variance at individual data locations:

noiseVar = np.power((0.3*np.abs(Y)), 2)

Add a random noise to the model responses:

Y = Y + np.sqrt(noiseVar) * np.random.randn(Y.shape[0],1)

Select the metamodeling tool and the Kriging module:

MetaOpts = {
    'Type': 'Metamodel',
    'MetaType': 'Kriging'
}

Use the experimental design and corresponding model responses generated earlier:

MetaOpts['ExpDesign'] = {
    'X': X.tolist(),
    'Y': Y.tolist()
}

Impose the known heteroscedastic noise variance:

MetaOpts['Regression'] = {
    'SigmaNSQ': noiseVar.tolist()
}

Create the Kriging regression model:

myKrigingRegression3 = uq.createModel(MetaOpts);

Visualize the results and add an error bar to each data point:

fig = uq.display(myKrigingRegression3)
fig.update_yaxes(range=[-10, 15])
for t in fig.data:
    if t.name.lower() == 'observations':
        t.error_y = {
            'type':'data',
            'array': 2*np.sqrt(noiseVar)*np.ones(Y.shape),
            'color': 'lightgrey'
        }
fig.show()

Terminate the remote UQCloud session

mySession.quit()