diff --git a/examples/model-comparison/L2_model.py b/examples/model-comparison/L2_model.py index e945c6888593d257d68d5af258c00a2020ca4c83..6b28c818101e25859bdb222b82cfd9bee741d381 100644 --- a/examples/model-comparison/L2_model.py +++ b/examples/model-comparison/L2_model.py @@ -3,14 +3,14 @@ """ This is a simple linear model. -The code for this numerical experiments is available at https://github.com/MichaelSinsbeck/ -paper sequential-design-model-selection. +The code for this numerical experiments is available at +https://github.com/MichaelSinsbeck/paper_sequential-design-model-selection. Author: Farid Mohammadi, M.Sc. E-Mail: farid.mohammadi@iws.uni-stuttgart.de Department of Hydromechanics and Modelling of Hydrosystems (LH2) -Institute for Modelling Hydraulic and Environmental Systems (IWS), University of Stuttgart -www.iws.uni-stuttgart.de/lh2/ +Institute for Modelling Hydraulic and Environmental Systems (IWS), +University of Stuttgart, www.iws.uni-stuttgart.de/lh2/ Pfaffenwaldring 61 70569 Stuttgart @@ -24,7 +24,11 @@ def L2_model(xx): """ Linear model y = a*x+b - Models adapted from Anneli Guthke + Models adapted from Anneli Guthke's paper: + ch€oniger, A., T. W€ohling, L. Samaniego,and W. Nowak (2014), Model + selection on solid ground: Rigorous comparison ofnine ways to evaluate + Bayesian modelevidence,Water Resour. Res.,50,9484–9513, + doi:10.1002/2014WR016062 Parameters ---------- diff --git a/examples/model-comparison/NL2_model.py b/examples/model-comparison/NL2_model.py new file mode 100644 index 0000000000000000000000000000000000000000..5fd4820e76a9756b85b891b0d8272404e81d3361 --- /dev/null +++ b/examples/model-comparison/NL2_model.py @@ -0,0 +1,57 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +This is a nonlinear cosine model. + +The code for this numerical experiments is available at +https://github.com/MichaelSinsbeck/paper_sequential-design-model-selection. + +Author: Farid Mohammadi, M.Sc. +E-Mail: farid.mohammadi@iws.uni-stuttgart.de +Department of Hydromechanics and Modelling of Hydrosystems (LH2) +Institute for Modelling Hydraulic and Environmental Systems (IWS), +University of Stuttgart, www.iws.uni-stuttgart.de/lh2/ +Pfaffenwaldring 61 +70569 Stuttgart + +Created on Fri Oct 8 2021 + +""" +import numpy as np + + +def NL2_model(xx): + """ + Nonlinear model y = exp(a*x) + b + + Models adapted from Anneli Guthke's paper: + ch€oniger, A., T. W€ohling, L. Samaniego,and W. Nowak (2014), Model + selection on solid ground: Rigorous comparison ofnine ways to evaluate + Bayesian modelevidence,Water Resour. Res.,50,9484–9513, + doi:10.1002/2014WR016062 + + Parameters + ---------- + xx : array + Parameters a and b. + + Returns + ------- + 2D-array + The first row contains the measurement locations. + The second row contains the model outputs. + + """ + n_output = 15 + meas_loc = np.linspace(0.25, 4.75, n_output) + + # NL2_model + NL2_model = np.exp(xx[:, 0] * meas_loc) + xx[:, 1] + + # Output + output = { + 'x_values': meas_loc, + 'Z': NL2_model + } + + return output diff --git a/examples/model-comparison/NL4_model.py b/examples/model-comparison/NL4_model.py index 2a1cc14e154e0baf34c7515ae2d65e6522391dfb..5ca495306d9a6d277ec654a3efbdfb84bfc28ce1 100644 --- a/examples/model-comparison/NL4_model.py +++ b/examples/model-comparison/NL4_model.py @@ -3,14 +3,14 @@ """ This is a nonlinear cosine model. -The code for this numerical experiments is available at https://github.com/MichaelSinsbeck/ -paper sequential-design-model-selection. +The code for this numerical experiments is available at +https://github.com/MichaelSinsbeck/paper_sequential-design-model-selection. Author: Farid Mohammadi, M.Sc. E-Mail: farid.mohammadi@iws.uni-stuttgart.de Department of Hydromechanics and Modelling of Hydrosystems (LH2) -Institute for Modelling Hydraulic and Environmental Systems (IWS), University of Stuttgart -www.iws.uni-stuttgart.de/lh2/ +Institute for Modelling Hydraulic and Environmental Systems (IWS), +University of Stuttgart, www.iws.uni-stuttgart.de/lh2/ Pfaffenwaldring 61 70569 Stuttgart @@ -24,7 +24,11 @@ def NL4_model(xx): """ Nonlinear model y = a*cos(b*x+c)+d - Models adapted from Anneli Guthke + Models adapted from Anneli Guthke's paper: + ch€oniger, A., T. W€ohling, L. Samaniego,and W. Nowak (2014), Model + selection on solid ground: Rigorous comparison ofnine ways to evaluate + Bayesian modelevidence,Water Resour. Res.,50,9484–9513, + doi:10.1002/2014WR016062 Parameters ---------- diff --git a/examples/model-comparison/test_model_comparison.py b/examples/model-comparison/test_model_comparison.py index a55d86afc8227d9dcf77a6e658ed21a4ea314235..4bdba4562f1b0fa5ca88970075f1697872211783 100644 --- a/examples/model-comparison/test_model_comparison.py +++ b/examples/model-comparison/test_model_comparison.py @@ -34,6 +34,7 @@ from bayes_inference.discrepancy import Discrepancy if __name__ == "__main__": + # Read data sigma = 0.6 data = { 'x [m]': np.linspace(0.25, 4.75, 15), @@ -49,16 +50,26 @@ if __name__ == "__main__": myL2Model.link_type = 'Function' myL2Model.py_file = 'L2_model' - myL2Model.name = 'L2_model' + myL2Model.name = 'linear' myL2Model.Output.names = ['Z'] myL2Model.observations = data + # -- Nonlinear exponential model ------- + myNL2Model = PyLinkForwardModel() + + myNL2Model.link_type = 'Function' + myNL2Model.py_file = 'NL2_model' + myNL2Model.name = 'exponential' + myNL2Model.Output.names = ['Z'] + myNL2Model.observations = data + # ------ Nonlinear cosine model --------- + # Data generating process myNL4Model = PyLinkForwardModel() myNL4Model.link_type = 'Function' myNL4Model.py_file = 'NL4_model' - myNL4Model.name = 'NL4_model' + myNL4Model.name = 'cosine' myNL4Model.Output.names = ['Z'] myNL4Model.observations = data @@ -83,14 +94,30 @@ if __name__ == "__main__": L2_Inputs.Marginals[i].name = f'$X_{i+1}$' L2_Inputs.Marginals[i].input_data = L2_input_params[:, i] + # ------ Nonlinear exponential model --------- + NL2_Inputs = Input() + NL2_prior_mean = np.array([0.4, -0.3]) + NL2_prior_cov = np.array( + [[0.003, -0.0001], + [-0.0001, 0.03]] + ) + NL2_input_params = np.random.multivariate_normal( + NL2_prior_mean, NL2_prior_cov, size=n_sample + ) + + for i in range(NL2_input_params.shape[1]): + NL2_Inputs.add_marginals() + NL2_Inputs.Marginals[i].name = f'$X_{i+1}$' + NL2_Inputs.Marginals[i].input_data = NL2_input_params[:, i] + # ------ Nonlinear cosine model --------- NL4_Inputs = Input() NL4_prior_mean = np.array([2.6, 0.5, -2.8, 2.3]) NL4_prior_cov = np.array( - [[0.46, -0.07, 0.24, -0.14], - [-0.07, 0.04, -0.05, 0.02], - [0.24, -0.05, 0.30, -0.16], - [-0.14, 0.02, -0.16, 0.30]] + [[0.44, -0.07, 0.24, -0.14], + [-0.07, 0.02, -0.05, 0.02], + [0.24, -0.05, 0.21, -0.16], + [-0.14, 0.02, -0.16, 0.28]] ) NL4_input_params = np.random.multivariate_normal( NL4_prior_mean, NL4_prior_cov, size=n_sample @@ -151,18 +178,22 @@ if __name__ == "__main__": # 6) chebyshev(FT) 7) grid(FT) 8)user L2_MetaModelOpts.ExpDesign.sampling_method = 'latin_hypercube' + # ------ Nonlinear cosine model --------- + NL2_MetaModelOpts = L2_MetaModelOpts.copy_meta_model_opts(NL2_Inputs) + # ------ Nonlinear cosine model --------- NL4_MetaModelOpts = L2_MetaModelOpts.copy_meta_model_opts(NL4_Inputs) # >>>>>> Train the Surrogates <<<<<<<<<<< L2_MetaModel = L2_MetaModelOpts.create_metamodel(myL2Model) + NL2_MetaModel = NL2_MetaModelOpts.create_metamodel(myNL2Model) NL4_MetaModel = NL4_MetaModelOpts.create_metamodel(myNL4Model) # ===================================================== # ========= POST PROCESSING OF METAMODELS =========== # ===================================================== # ---------- Linear model ------------- - L2_PostPCE = PostProcessing(L2_MetaModel) + L2_PostPCE = PostProcessing(L2_MetaModel, name=myL2Model.name) # Plot to check validation visually. L2_PostPCE.valid_metamodel(n_samples=3) @@ -176,8 +207,23 @@ if __name__ == "__main__": # Plot the sobol indices sobol_cell, total_sobol = L2_PostPCE.sobol_indices() + # ---------- Linear model ------------- + NL2_PostPCE = PostProcessing(NL2_MetaModel, name=myNL2Model.name) + + # Plot to check validation visually. + NL2_PostPCE.valid_metamodel(n_samples=3) + + # Plot moments + NL2_PostPCE.plot_moments() + + # Compute and print RMSE error + NL2_PostPCE.check_accuracy(n_samples=3000) + + # Plot the sobol indices + sobol_cell, total_sobol = NL2_PostPCE.sobol_indices() + # ------ Nonlinear cosine model --------- - NL4_PostPCE = PostProcessing(NL4_MetaModel) + NL4_PostPCE = PostProcessing(NL4_MetaModel, name=myNL4Model.name) # Plot to check validation visually. NL4_PostPCE.valid_metamodel(n_samples=3) @@ -202,8 +248,9 @@ if __name__ == "__main__": # ----- Define the options model ------- metaModels = { - "L2_model": L2_MetaModel, - "NL4_model": NL4_MetaModel + "linear": L2_MetaModel, + "exponential": NL4_MetaModel, + "cosine": NL4_MetaModel } # MCMC inference method @@ -224,7 +271,7 @@ if __name__ == "__main__": # BME Bootstrap OptsDict_Bootstrap = { "bootstrap": True, - "n_samples": 20000, + "n_samples": 500000, "Discrepancy": DiscrepancyOpts, "emulator": False, "plot_post_pred": False