Evaluating the conditional covariance per point

conditional_covariance_marginal.py

import openturns as ot

def getConditionalMarginalCovariance(krigingResult, x):

    def _getConditionalMarginalCovariancePoint(x):
        cov = krigingResult.getConditionalCovariance(x)
        return cov

    def _getConditionalMarginalCovarianceSample(x):
        cov_coll = []
        for k in range(len(x)):
            cov_coll.append( _getConditionalMarginalCovariancePoint(x[k]) )
        return cov_coll
    
    try: 
        data = ot.Point(x)
        return _getConditionalMarginalCovariancePoint(data)
    except TypeError:
        pass

    data = ot.Sample(x)
    return _getConditionalMarginalCovarianceSample(data)


if not hasattr(ot.KrigingResult, "getConditionalMarginalCovariance"):
    ot.KrigingResult.getConditionalMarginalCovariance = getConditionalMarginalCovariance    

if __name__ == "__main__":
    import numpy as np
    import time

    g = ot.SymbolicFunction(['x'], ['sin(x)'])

    x_train = ot.Sample([1.,3.,4.,6.,7.9,11., 11.5],1)
    y_train = g(x_train)
    n_train = len(x_train)


    # Learning kriging model
    xmin = 0.
    xmax = 12.
    n_test = 10000
    step = (xmax-xmin)/(n_test-1)
    myRegularGrid = ot.RegularGrid(xmin, step, n_test)
    x_test_coord = myRegularGrid.getValues()
    x_test = ot.Sample(x_test_coord, 1)
    y_test = g(x_test)



    # Matern covariance model & constant basis
    dimension = 1
    basis = ot.ConstantBasisFactory(dimension).build()
    covarianceModel = ot.MaternModel([1.]*dimension, 1.5)
    algo = ot.KrigingAlgorithm(x_train, y_train, covarianceModel, basis)
    algo.run()
    result = algo.getResult()


    # Surrogate model
    krigeageMM = result.getMetaModel()
    y_test_MM = krigeageMM(x_test)


    # Conditional covariance : (n_test x n_test) matrix
    tic = time.time()
    covGrid = result.getConditionalCovariance(x_test)
    toc = time.time()
    conditionalStDev = ot.Sample([np.sqrt(covGrid[i, i]) for i in range(n_test)],1)
    print("Evaluating conditional full cov matrix in %1.3f s"%(toc - tic))


    # Covariance without accounting interaction between x data:
    tic = time.time()
    coll_cov_conditional_cov = result.getConditionalMarginalCovariance(x_test)
    toc = time.time()
    print("Evaluating conditional diag cov only matrix in %1.3f s"%(toc - tic))

    # Compare the two methods:
    print((conditionalStDev - ot.Sample([np.sqrt(cov[0,0]) for cov in coll_cov_conditional_cov],1)).getMax())