Skip to content
Snippets Groups Projects
Commit b827cffa authored by Hu Zhao's avatar Hu Zhao
Browse files

docs: track built auto_examples

parent 1a5b9ec0
Branches
No related tags found
No related merge requests found
Showing
with 3431 additions and 0 deletions
File added
File added
File added
:orphan:
Emulator
========
.. raw:: html
<div class="sphx-glr-thumbnails">
.. raw:: html
<div class="sphx-glr-thumbcontainer" tooltip="PPGaSP: GP emulation for multi-output functions">
.. only:: html
.. image:: /auto_examples/emulator/images/thumb/sphx_glr_plot_ppgasp_thumb.png
:alt: PPGaSP: GP emulation for multi-output functions
:ref:`sphx_glr_auto_examples_emulator_plot_ppgasp.py`
.. raw:: html
<div class="sphx-glr-thumbnail-title">PPGaSP: GP emulation for multi-output functions</div>
</div>
.. raw:: html
<div class="sphx-glr-thumbcontainer" tooltip="ScalarGaSP: GP emulation for a single-output function">
.. only:: html
.. image:: /auto_examples/emulator/images/thumb/sphx_glr_plot_scalargasp_thumb.png
:alt: ScalarGaSP: GP emulation for a single-output function
:ref:`sphx_glr_auto_examples_emulator_plot_scalargasp.py`
.. raw:: html
<div class="sphx-glr-thumbnail-title">ScalarGaSP: GP emulation for a single-output function</div>
</div>
.. raw:: html
</div>
.. toctree::
:hidden:
/auto_examples/emulator/plot_ppgasp
/auto_examples/emulator/plot_scalargasp
.. only:: html
.. container:: sphx-glr-footer sphx-glr-footer-gallery
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download all examples in Python source code: emulator_python.zip </auto_examples/emulator/emulator_python.zip>`
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download all examples in Jupyter notebooks: emulator_jupyter.zip </auto_examples/emulator/emulator_jupyter.zip>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/emulator/plot_ppgasp.py"
.. LINE NUMBERS ARE GIVEN BELOW.
.. only:: html
.. note::
:class: sphx-glr-download-link-note
Click :ref:`here <sphx_glr_download_auto_examples_emulator_plot_ppgasp.py>`
to download the full example code
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_auto_examples_emulator_plot_ppgasp.py:
PPGaSP: GP emulation for multi-output functions
===============================================
.. GENERATED FROM PYTHON SOURCE LINES 9-23
This example shows how to apply Gaussian process emulation to a multi-output
function using class :class:`.PPGaSP`.
The multi-output function that we are going to look at is the `DIAMOND`
(diplomatic and military operations in a non-warfighting domain) computer model.
It is used as a testbed to illustrate the `PP GaSP` emulator in the `R` package
`RobustGaSP`, see :cite:t:`Gu2019` for more detail.
The simulator has :math:`13` input parameters and :math:`5` outputs. Namely,
:math:`\mathbf{y}=f(\mathbf{x})` where :math:`\mathbf{x}=(x_1,\ldots,x_{13})^T`
and :math:`\mathbf{y}=(y_1,\ldots,y_5)^T`.
The training and testing data are provided in the folder '.../tests/data/'.
We first load the training and testing data.
.. GENERATED FROM PYTHON SOURCE LINES 24-40
.. code-block:: default
import numpy as np
import os
dir_data = os.path.abspath('../../../tests/data/')
humanityX = np.genfromtxt(os.path.join(dir_data, 'humanityX.csv'), delimiter=',')
humanityY = np.genfromtxt(os.path.join(dir_data, 'humanityY.csv'), delimiter=',')
print(f"Number of training data points: ", humanityX.shape[0])
print(f"Input dimension: ", humanityX.shape[1])
print(f"Output dimension: ", humanityY.shape[1])
humanityXt = np.genfromtxt(os.path.join(dir_data, 'humanityXt.csv'), delimiter=',')
humanityYt = np.genfromtxt(os.path.join(dir_data, 'humanityYt.csv'), delimiter=',')
print(f"Number of testing data points: ", humanityXt.shape[0])
.. rst-class:: sphx-glr-script-out
.. code-block:: none
Number of training data points: 120
Input dimension: 13
Output dimension: 5
Number of testing data points: 120
.. GENERATED FROM PYTHON SOURCE LINES 41-43
.. note:: You may need to modify ``dir_data`` according to where you save them
on your local machine.
.. GENERATED FROM PYTHON SOURCE LINES 45-50
``humanityX`` and ``humanitY`` are the training data, corresponding to ``design``
and ``response`` respectively. ``humanityXt`` are testing input data, at which
we are going to make predictions once the emulator is trained. ``humanityYt``
are the true outputs at ``humanityXt``, which is then used to validate the
performance of the trained emulator.
.. GENERATED FROM PYTHON SOURCE LINES 52-54
To build a `PP GaSP` emulator for the above simulator, first import class
:class:`.PPGaSP` by
.. GENERATED FROM PYTHON SOURCE LINES 55-58
.. code-block:: default
from psimpy.emulator import PPGaSP
.. GENERATED FROM PYTHON SOURCE LINES 59-63
Then, create an instance of :class:`.PPGaSP`. The parameter ``ndim``
(dimension of function input ``x``) must be specified. Optional parameters, such
as ``method``, ``kernel_type``, etc., can be set up if desired. Here, we leave
all the optional parameters to their default values.
.. GENERATED FROM PYTHON SOURCE LINES 64-68
.. code-block:: default
emulator = PPGaSP(ndim=humanityX.shape[1])
.. GENERATED FROM PYTHON SOURCE LINES 69-71
Next, we train the `PP GaSP` emulator based on the training data using
:py:meth:`.PPGaSP.train` method.
.. GENERATED FROM PYTHON SOURCE LINES 72-75
.. code-block:: default
emulator.train(design=humanityX, response=humanityY)
.. rst-class:: sphx-glr-script-out
.. code-block:: none
The upper bounds of the range parameters are 296,974 297,1814 294,7672 295,9345 296,563 295,935 297,3198 296,4444 296,9323 298,2373 298,0619 298,7249 298,7249
The initial values of range parameters are 5,93948 5,943628 5,895344 5,918691 5,93126 5,918701 5,946395 5,928888 5,938645 5,964745 5,961238 5,974498 5,974498
Start of the optimization 1 :
The number of iterations is 44
The value of the marginal posterior function is -5279,534
Optimized range parameters are 25,46132 2,891428 5,30812 26,95519 291,1828 48,17274 84,56133 2,874385 39,12121 51,81304 0,5407651 1,728696 1,020605
Optimized nugget parameter is 0
Convergence: TRUE
The initial values of range parameters are 12,37503 12,38367 12,28307 12,33172 12,3579 12,33174 12,38944 12,35296 12,37329 12,42767 12,42037 12,44799 12,44799
Start of the optimization 2 :
The number of iterations is 45
The value of the marginal posterior function is -5279,534
Optimized range parameters are 25,46132 2,891428 5,308121 26,95519 291,183 48,17274 84,56133 2,874385 39,12121 51,81304 0,5407651 1,728696 1,020605
Optimized nugget parameter is 0
Convergence: TRUE
.. GENERATED FROM PYTHON SOURCE LINES 76-79
With the trained emulator, we can make predictions for any arbitrary set of
input points using :py:meth:`PPGaSP.predict` method.
Here, we make predictions at testing input points ``humanityXt``.
.. GENERATED FROM PYTHON SOURCE LINES 80-83
.. code-block:: default
predictions = emulator.predict(humanityXt)
.. GENERATED FROM PYTHON SOURCE LINES 84-86
We can validate the performance of the trained emulator based on the true outputs
``humanityYt``.
.. GENERATED FROM PYTHON SOURCE LINES 87-106
.. code-block:: default
import matplotlib.pyplot as plt
fig, ax = plt.subplots(5, 1, figsize=(6,15))
for i in range(humanityY.shape[1]):
ax[i].set_xlabel(f'Actual $y_{i+1}$')
ax[i].set_ylabel(f'Emulator-predicted $y_{i+1}$')
ax[i].set_xlim(np.min(humanityYt[:,i]),np.max(humanityYt[:,i]))
ax[i].set_ylim(np.min(humanityYt[:,i]),np.max(humanityYt[:,i]))
_ = ax[i].plot([np.min(humanityYt[:,i]),np.max(humanityYt[:,i])], [np.min(humanityYt[:,i]),np.max(humanityYt[:,i])])
_ = ax[i].errorbar(humanityYt[:,i], predictions[:,i,0], predictions[:,i,3], fmt='.', linestyle='', label='prediction and std')
_ = ax[i].legend()
plt.tight_layout()
.. image-sg:: /auto_examples/emulator/images/sphx_glr_plot_ppgasp_001.png
:alt: plot ppgasp
:srcset: /auto_examples/emulator/images/sphx_glr_plot_ppgasp_001.png
:class: sphx-glr-single-img
.. GENERATED FROM PYTHON SOURCE LINES 107-109
We can also draw any number of samples at testing input ``humanityXt`` using
:py:meth:`.PPGaSPsample()` method.
.. GENERATED FROM PYTHON SOURCE LINES 110-112
.. code-block:: default
samples = emulator.sample(humanityXt, nsamples=10)
print("Shape of samples: ", samples.shape)
.. rst-class:: sphx-glr-script-out
.. code-block:: none
Shape of samples: (120, 5, 10)
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 0 minutes 5.333 seconds)
.. _sphx_glr_download_auto_examples_emulator_plot_ppgasp.py:
.. only:: html
.. container:: sphx-glr-footer sphx-glr-footer-example
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_ppgasp.py <plot_ppgasp.py>`
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_ppgasp.ipynb <plot_ppgasp.ipynb>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/emulator/plot_scalargasp.py"
.. LINE NUMBERS ARE GIVEN BELOW.
.. only:: html
.. note::
:class: sphx-glr-download-link-note
Click :ref:`here <sphx_glr_download_auto_examples_emulator_plot_scalargasp.py>`
to download the full example code
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_auto_examples_emulator_plot_scalargasp.py:
ScalarGaSP: GP emulation for a single-output function
=====================================================
.. GENERATED FROM PYTHON SOURCE LINES 8-13
This example shows how to apply Gaussian process emulation to a single-output
function using class :class:`.ScalarGaSP`.
The task is to build a GP emulator for the function :math:`y = x * sin(x)`
based on a few number of training data.
.. GENERATED FROM PYTHON SOURCE LINES 16-17
First, import the class :class:`.ScalarGaSP` by
.. GENERATED FROM PYTHON SOURCE LINES 18-21
.. code-block:: default
from psimpy.emulator import ScalarGaSP
.. GENERATED FROM PYTHON SOURCE LINES 22-26
Then, create an instance of :class:`.ScalarGaSP`. The parameter ``ndim``
(dimension of function input ``x``) must be specified. Optional parameters, such
as ``method``, ``kernel_type``, etc., can be set up if desired. Here, we leave
all the optional parameters to their default values.
.. GENERATED FROM PYTHON SOURCE LINES 27-30
.. code-block:: default
emulator = ScalarGaSP(ndim=1)
.. GENERATED FROM PYTHON SOURCE LINES 31-34
Given training input points ``design`` and corresponding output values ``response``,
the emulator can be trained using :py:meth:`.ScalarGaSP.train`. Below we train
an emulator using :math:`8` selected points.
.. GENERATED FROM PYTHON SOURCE LINES 35-47
.. code-block:: default
import numpy as np
def f(x):
#return x + 3*np.sin(x/2)
return x*np.sin(x)
x = np.arange(2,10,1)
y = f(x)
emulator.train(design=x, response=y)
.. rst-class:: sphx-glr-script-out
.. code-block:: none
The upper bounds of the range parameters are 1881,404
The initial values of range parameters are 37,62809
Start of the optimization 1 :
The number of iterations is 9
The value of the marginal posterior function is -13,79116
Optimized range parameters are 2,237467
Optimized nugget parameter is 0
Convergence: TRUE
The initial values of range parameters are 0,21875
Start of the optimization 2 :
The number of iterations is 9
The value of the marginal posterior function is -13,79116
Optimized range parameters are 2,237467
Optimized nugget parameter is 0
Convergence: TRUE
.. GENERATED FROM PYTHON SOURCE LINES 48-50
We can validate the performance of the trained emulator using the leave-one-out
cross validation method :py:meth:`loo_validate()`.
.. GENERATED FROM PYTHON SOURCE LINES 51-54
.. code-block:: default
validation = emulator.loo_validate()
.. GENERATED FROM PYTHON SOURCE LINES 55-57
Let's plot emulator predictions vs actual outputs. The error bar indicates the
standard deviation.
.. GENERATED FROM PYTHON SOURCE LINES 58-74
.. code-block:: default
import matplotlib.pyplot as plt
fig , ax = plt.subplots(figsize=(4,4))
ax.set_xlabel('Actual y')
ax.set_ylabel('Emulator-predicted y')
ax.set_xlim(np.min(y)-1,np.max(y)+1)
ax.set_ylim(np.min(y)-1,np.max(y)+1)
_ = ax.plot([np.min(y)-1,np.max(y)+1], [np.min(y)-1,np.max(y)+1])
_ = ax.errorbar(y, validation[:,0], validation[:,1], fmt='.', linestyle='', label='prediction and std')
_ = plt.legend()
plt.tight_layout()
.. image-sg:: /auto_examples/emulator/images/sphx_glr_plot_scalargasp_001.png
:alt: plot scalargasp
:srcset: /auto_examples/emulator/images/sphx_glr_plot_scalargasp_001.png
:class: sphx-glr-single-img
.. GENERATED FROM PYTHON SOURCE LINES 75-80
With the trained emulator at our deposit, we can use the
:py:meth:`.ScalarGaSP.predict()` method to make predictions at
any arbitrary set of input points (``testing_input``). It should be noted that,
``testing_trend`` should be set according to ``trend`` used during emulator
training.
.. GENERATED FROM PYTHON SOURCE LINES 81-95
.. code-block:: default
testing_input = np.arange(0,10,0.1)
predictions = emulator.predict(testing_input)
plt.plot(testing_input, predictions[:, 0], 'r-', label= "mean")
plt.scatter(x, y, s=15, c='k', label="training data", zorder=3)
plt.plot(testing_input, f(testing_input), 'k:', zorder=2, label="true function")
plt.fill_between(testing_input, predictions[:, 1], predictions[:, 2], alpha=0.3, label="95% CI")
plt.xlabel('x')
plt.ylabel('emulator-predicted y')
plt.xlim(testing_input[0], testing_input[-1])
_ = plt.legend()
plt.tight_layout()
.. image-sg:: /auto_examples/emulator/images/sphx_glr_plot_scalargasp_002.png
:alt: plot scalargasp
:srcset: /auto_examples/emulator/images/sphx_glr_plot_scalargasp_002.png
:class: sphx-glr-single-img
.. GENERATED FROM PYTHON SOURCE LINES 96-98
We can also draw any number of samples at ``testing_input``` using
:py:meth:`.ScalarGaSP.sample()` method.
.. GENERATED FROM PYTHON SOURCE LINES 99-117
.. code-block:: default
nsamples = 5
samples = emulator.sample(testing_input, nsamples=nsamples)
# sphinx_gallery_thumbnail_number = 3
for i in range(nsamples):
plt.plot(testing_input, samples[:,i], '--', label=f'sample{i+1}')
plt.scatter(x, y, s=15, c='k', label="training data", zorder=2)
plt.plot(testing_input, f(testing_input), 'k:', zorder=2, label="true function")
plt.fill_between(testing_input, predictions[:, 1], predictions[:, 2], alpha=0.3, label="95% CI")
plt.xlabel('x')
plt.ylabel('emulator-predicted y')
plt.xlim(testing_input[0], testing_input[-1])
_ = plt.legend()
plt.tight_layout()
.. image-sg:: /auto_examples/emulator/images/sphx_glr_plot_scalargasp_003.png
:alt: plot scalargasp
:srcset: /auto_examples/emulator/images/sphx_glr_plot_scalargasp_003.png
:class: sphx-glr-single-img
.. GENERATED FROM PYTHON SOURCE LINES 118-124
.. tip:: Above example shows how to train a GP emulator based on noise-free training data,
which is often the case of emulating a deterministic simulator. If you are dealing
with noisy training data, you can
- set the parameter ``nugget`` to a desired value, or
- set ``nugget`` to :math:`0` and ``nugget_est`` to `True`, meaning that ``nugget``
is estimated from the noisy training data.
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 0 minutes 0.728 seconds)
.. _sphx_glr_download_auto_examples_emulator_plot_scalargasp.py:
.. only:: html
.. container:: sphx-glr-footer sphx-glr-footer-example
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_scalargasp.py <plot_scalargasp.py>`
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_scalargasp.ipynb <plot_scalargasp.ipynb>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
:orphan:
.. _sphx_glr_auto_examples_emulator_sg_execution_times:
Computation times
=================
**00:06.062** total execution time for **auto_examples_emulator** files:
+------------------------------------------------------------------------------------+-----------+--------+
| :ref:`sphx_glr_auto_examples_emulator_plot_ppgasp.py` (``plot_ppgasp.py``) | 00:05.333 | 0.0 MB |
+------------------------------------------------------------------------------------+-----------+--------+
| :ref:`sphx_glr_auto_examples_emulator_plot_scalargasp.py` (``plot_scalargasp.py``) | 00:00.728 | 0.0 MB |
+------------------------------------------------------------------------------------+-----------+--------+
:orphan:
Inference
=========
.. raw:: html
<div class="sphx-glr-thumbnails">
.. raw:: html
<div class="sphx-glr-thumbcontainer" tooltip="Bayesian inference">
.. only:: html
.. image:: /auto_examples/inference/images/thumb/sphx_glr_plot_bayes_inference_thumb.png
:alt: Bayesian inference
:ref:`sphx_glr_auto_examples_inference_plot_bayes_inference.py`
.. raw:: html
<div class="sphx-glr-thumbnail-title">Bayesian inference</div>
</div>
.. raw:: html
<div class="sphx-glr-thumbcontainer" tooltip="Active learning">
.. only:: html
.. image:: /auto_examples/inference/images/thumb/sphx_glr_plot_active_learning_thumb.png
:alt: Active learning
:ref:`sphx_glr_auto_examples_inference_plot_active_learning.py`
.. raw:: html
<div class="sphx-glr-thumbnail-title">Active learning</div>
</div>
.. raw:: html
</div>
.. toctree::
:hidden:
/auto_examples/inference/plot_bayes_inference
/auto_examples/inference/plot_active_learning
.. only:: html
.. container:: sphx-glr-footer sphx-glr-footer-gallery
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download all examples in Python source code: inference_python.zip </auto_examples/inference/inference_python.zip>`
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download all examples in Jupyter notebooks: inference_jupyter.zip </auto_examples/inference/inference_jupyter.zip>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/inference/plot_bayes_inference.py"
.. LINE NUMBERS ARE GIVEN BELOW.
.. only:: html
.. note::
:class: sphx-glr-download-link-note
Click :ref:`here <sphx_glr_download_auto_examples_inference_plot_bayes_inference.py>`
to download the full example code
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_auto_examples_inference_plot_bayes_inference.py:
Bayesian inference
==================
.. GENERATED FROM PYTHON SOURCE LINES 8-16
This example shows how to perform Bayesian inference given the uniform prior
:math:`p(\mathbf{x})=p(x_1,x_2)=0.01`
where :math:`x_i \in [-5,5], i=1,2`, and likelihood
:math:`L(\mathbf{x}|\mathbf{d})=\exp \left(-\frac{1}{100}\left(x_1-1\right)^2-\left(x_1^2-x_2\right)^2\right)`.
.. GENERATED FROM PYTHON SOURCE LINES 17-29
.. code-block:: default
import numpy as np
ndim = 2
bounds = np.array([[-5,5],[-5,5]])
def prior(x):
return 0.01
def likelihood(x):
return np.exp(-(x[0]-1)**2/100 - (x[0]**2-x[1])**2)
.. GENERATED FROM PYTHON SOURCE LINES 30-33
To estimate the posterior using grid estimation, we need to import the
:class:`.GridEstimation` class, create an instance, and call the
:py:meth:`.GridEstimation.run` method.
.. GENERATED FROM PYTHON SOURCE LINES 34-40
.. code-block:: default
from psimpy.inference.bayes_inference import GridEstimation
grid_estimator = GridEstimation(ndim, bounds, prior, likelihood)
posterior, x_ndim = grid_estimator.run(nbins=[50,40])
.. GENERATED FROM PYTHON SOURCE LINES 41-43
The following figure plots the estimated posterior.
.. GENERATED FROM PYTHON SOURCE LINES 44-61
.. code-block:: default
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1,figsize=(6,4))
# mask insignificant values
posterior = np.where(posterior < 1e-10, None, posterior)
contour = ax.contour(x_ndim[0], x_ndim[1], np.transpose(posterior), levels=10)
plt.colorbar(contour, ax=ax)
ax.set_xlim(bounds[0,0], bounds[0,1])
ax.set_ylim(bounds[1,0], bounds[1,1])
ax.set_title('Grid estimation')
ax.set_xlabel('x1')
ax.set_ylabel('x2')
plt.tight_layout()
.. image-sg:: /auto_examples/inference/images/sphx_glr_plot_bayes_inference_001.png
:alt: Grid estimation
:srcset: /auto_examples/inference/images/sphx_glr_plot_bayes_inference_001.png
:class: sphx-glr-single-img
.. GENERATED FROM PYTHON SOURCE LINES 62-67
To estimate the posterior using Metropolis Hastings estimation, we need to import
the :class:`.MetropolisHastingsEstimation` class, create an instance, and call the
:py:meth:`.MetropolisHastingsEstimation.run` method. The
:py:meth:`.MetropolisHastingsEstimation.run` method has a parameter, ``mh_sampler``,
which takes an instance of :class:`.MetropolisHastings` as argument.
.. GENERATED FROM PYTHON SOURCE LINES 68-92
.. code-block:: default
from psimpy.inference.bayes_inference import MetropolisHastingsEstimation
mh_estimator = MetropolisHastingsEstimation(ndim, bounds, prior, likelihood)
# create a mh_sampler
from psimpy.sampler.metropolis_hastings import MetropolisHastings
from scipy.stats import multivariate_normal
init_state = np.array([-4,-4])
f_sample = multivariate_normal.rvs
nburn = 100
nthin = 10
seed = 1
kwgs_f_sample = {'random_state': np.random.default_rng(seed)}
mh_sampler = MetropolisHastings(ndim=ndim, init_state=init_state,
f_sample=f_sample, bounds=bounds, nburn=nburn, nthin=nthin, seed=seed,
kwgs_f_sample=kwgs_f_sample)
nsamples = 5000
mh_samples, mh_accept = mh_estimator.run(nsamples, mh_sampler)
.. GENERATED FROM PYTHON SOURCE LINES 93-96
The following figure plots the samples drawn from the unnormalized posterior,
which can be used to estimate the posterior and its poperties.
.. GENERATED FROM PYTHON SOURCE LINES 97-110
.. code-block:: default
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1,figsize=(5,4))
ax.scatter(mh_samples[:,0], mh_samples[:,1], s=10, c='r', marker='o', alpha=0.1)
ax.set_xlim(bounds[0,0], bounds[0,1])
ax.set_ylim(bounds[1,0], bounds[1,1])
ax.set_title('MH estimation')
ax.set_xlabel('x1')
ax.set_ylabel('x2')
plt.tight_layout()
.. image-sg:: /auto_examples/inference/images/sphx_glr_plot_bayes_inference_002.png
:alt: MH estimation
:srcset: /auto_examples/inference/images/sphx_glr_plot_bayes_inference_002.png
:class: sphx-glr-single-img
.. GENERATED FROM PYTHON SOURCE LINES 111-116
.. note:: Besides ``prior`` and ``likelihood``, one can also instantiate
the :class:`.MetropolisHastingsEstimation` class with
- ``ln_prior`` and ``ln_likelihood``: Natural logarithm of ``prior`` and ``likelihood``.
- ``ln_pxl``: Natural logarithm of the product of ``prior`` and ``likelihood``.
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 0 minutes 9.739 seconds)
.. _sphx_glr_download_auto_examples_inference_plot_bayes_inference.py:
.. only:: html
.. container:: sphx-glr-footer sphx-glr-footer-example
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_bayes_inference.py <plot_bayes_inference.py>`
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_bayes_inference.ipynb <plot_bayes_inference.ipynb>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
:orphan:
.. _sphx_glr_auto_examples_inference_sg_execution_times:
Computation times
=================
**02:22.392** total execution time for **auto_examples_inference** files:
+-----------------------------------------------------------------------------------------------+-----------+--------+
| :ref:`sphx_glr_auto_examples_inference_plot_active_learning.py` (``plot_active_learning.py``) | 02:12.653 | 0.0 MB |
+-----------------------------------------------------------------------------------------------+-----------+--------+
| :ref:`sphx_glr_auto_examples_inference_plot_bayes_inference.py` (``plot_bayes_inference.py``) | 00:09.739 | 0.0 MB |
+-----------------------------------------------------------------------------------------------+-----------+--------+
:orphan:
Sampler
=======
.. raw:: html
<div class="sphx-glr-thumbnails">
.. raw:: html
<div class="sphx-glr-thumbcontainer" tooltip="Saltelli sampling">
.. only:: html
.. image:: /auto_examples/sampler/images/thumb/sphx_glr_plot_saltelli_thumb.png
:alt: Saltelli sampling
:ref:`sphx_glr_auto_examples_sampler_plot_saltelli.py`
.. raw:: html
<div class="sphx-glr-thumbnail-title">Saltelli sampling</div>
</div>
.. raw:: html
<div class="sphx-glr-thumbcontainer" tooltip="Latin hypercube sampling">
.. only:: html
.. image:: /auto_examples/sampler/images/thumb/sphx_glr_plot_latin_thumb.png
:alt: Latin hypercube sampling
:ref:`sphx_glr_auto_examples_sampler_plot_latin.py`
.. raw:: html
<div class="sphx-glr-thumbnail-title">Latin hypercube sampling</div>
</div>
.. raw:: html
<div class="sphx-glr-thumbcontainer" tooltip="Metropolis Hastings sampling">
.. only:: html
.. image:: /auto_examples/sampler/images/thumb/sphx_glr_plot_metropolis_hastings_thumb.png
:alt: Metropolis Hastings sampling
:ref:`sphx_glr_auto_examples_sampler_plot_metropolis_hastings.py`
.. raw:: html
<div class="sphx-glr-thumbnail-title">Metropolis Hastings sampling</div>
</div>
.. raw:: html
</div>
.. toctree::
:hidden:
/auto_examples/sampler/plot_saltelli
/auto_examples/sampler/plot_latin
/auto_examples/sampler/plot_metropolis_hastings
.. only:: html
.. container:: sphx-glr-footer sphx-glr-footer-gallery
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download all examples in Python source code: sampler_python.zip </auto_examples/sampler/sampler_python.zip>`
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download all examples in Jupyter notebooks: sampler_jupyter.zip </auto_examples/sampler/sampler_jupyter.zip>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/sampler/plot_latin.py"
.. LINE NUMBERS ARE GIVEN BELOW.
.. only:: html
.. note::
:class: sphx-glr-download-link-note
Click :ref:`here <sphx_glr_download_auto_examples_sampler_plot_latin.py>`
to download the full example code
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_auto_examples_sampler_plot_latin.py:
Latin hypercube sampling
========================
.. GENERATED FROM PYTHON SOURCE LINES 8-11
This example shows how to draw samples using Latin hypercube sampling.
.. GENERATED FROM PYTHON SOURCE LINES 14-17
For the illustration purpose, let's have a look at a two-dimensional example
where we have two random variables X and Y. Each is uniformly distributed in its
range.
.. GENERATED FROM PYTHON SOURCE LINES 18-25
.. code-block:: default
import numpy as np
ndim = 2
# range of X is 10 to 20, range of Y is -10 to 0
bounds = np.array([[10,20], [-10,0]])
.. GENERATED FROM PYTHON SOURCE LINES 26-28
Given this setting, we can import :class:`.Latin`, create an instance, and
call the :py:meth:`.Latin.sample` method to draw required number of samples
.. GENERATED FROM PYTHON SOURCE LINES 29-36
.. code-block:: default
from psimpy.sampler import LHS
# setting seed leads to same samples every time when the codes are run
lhs_sampler = LHS(ndim, bounds, seed=10)
lhs_samples = lhs_sampler.sample(nsamples=5)
.. GENERATED FROM PYTHON SOURCE LINES 37-38
The samples are plotted in the following figure.
.. GENERATED FROM PYTHON SOURCE LINES 39-54
.. code-block:: default
import matplotlib.pyplot as plt
fig , ax = plt.subplots(figsize=(6,4))
ax.scatter(lhs_samples[:,0], lhs_samples[:,1], s=10, c='blue', marker='o')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_xlim(bounds[0])
ax.set_ylim(bounds[1])
ax.set_title("Latin hypercube samples (criterion='random')")
_ = ax.grid(visible=True, which='major', axis='both')
.. image-sg:: /auto_examples/sampler/images/sphx_glr_plot_latin_001.png
:alt: Latin hypercube samples (criterion='random')
:srcset: /auto_examples/sampler/images/sphx_glr_plot_latin_001.png
:class: sphx-glr-single-img
.. GENERATED FROM PYTHON SOURCE LINES 55-58
There are different criterions to pick samples in each hypercube. The default
is `random`, as used above. Other options are `center` and `maximin`. For instance,
we can use the `center` criterion to draw :math:`5` samples as follows:
.. GENERATED FROM PYTHON SOURCE LINES 59-75
.. code-block:: default
lhs_sampler = LHS(ndim, bounds, criterion='center', seed=10)
lhs_samples = lhs_sampler.sample(nsamples=5)
fig , ax = plt.subplots(figsize=(6,4))
ax.scatter(lhs_samples[:,0], lhs_samples[:,1], s=10, c='blue', marker='o')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_xlim(bounds[0])
ax.set_ylim(bounds[1])
ax.set_title("Latin hypercube samples (criterion='center')")
_ = ax.grid(visible=True, which='major', axis='both')
.. image-sg:: /auto_examples/sampler/images/sphx_glr_plot_latin_002.png
:alt: Latin hypercube samples (criterion='center')
:srcset: /auto_examples/sampler/images/sphx_glr_plot_latin_002.png
:class: sphx-glr-single-img
.. GENERATED FROM PYTHON SOURCE LINES 76-77
And we can use the `maximin` criterion as follows:
.. GENERATED FROM PYTHON SOURCE LINES 78-91
.. code-block:: default
lhs_sampler = LHS(ndim, bounds, criterion='maximin', seed=10, iteration=500)
lhs_samples = lhs_sampler.sample(nsamples=5)
fig , ax = plt.subplots(figsize=(6,4))
ax.scatter(lhs_samples[:,0], lhs_samples[:,1], s=10, c='blue', marker='o')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_xlim(bounds[0])
ax.set_ylim(bounds[1])
ax.set_title("Latin hypercube samples (criterion='maximin')")
_ = ax.grid(visible=True, which='major', axis='both')
.. image-sg:: /auto_examples/sampler/images/sphx_glr_plot_latin_003.png
:alt: Latin hypercube samples (criterion='maximin')
:srcset: /auto_examples/sampler/images/sphx_glr_plot_latin_003.png
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 0 minutes 0.270 seconds)
.. _sphx_glr_download_auto_examples_sampler_plot_latin.py:
.. only:: html
.. container:: sphx-glr-footer sphx-glr-footer-example
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_latin.py <plot_latin.py>`
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_latin.ipynb <plot_latin.ipynb>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/sampler/plot_metropolis_hastings.py"
.. LINE NUMBERS ARE GIVEN BELOW.
.. only:: html
.. note::
:class: sphx-glr-download-link-note
Click :ref:`here <sphx_glr_download_auto_examples_sampler_plot_metropolis_hastings.py>`
to download the full example code
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_auto_examples_sampler_plot_metropolis_hastings.py:
Metropolis Hastings sampling
============================
.. GENERATED FROM PYTHON SOURCE LINES 8-22
This example shows how to draw samples using Metropolis Hastings sampling.
The target probability distribution is
:math:`p(\mathbf{x})=p(x_1, x_2) \propto \exp \left(-\frac{1}{100}\left(x_1-1\right)^2-\left(x_1^2-x_2\right)^2\right)`
where :math:`x_1 \in [-5,5]` and :math:`x_2 \in [-5,5]`.
It should be noted that the right hand side of the equation is an unnormalized
probability density function since its integral is not equal to :math:`1`.
This can happen, for example, when the normalization constant is unknown or difficult
to compute.
We can define the target probability distribution in Python as follows:
.. GENERATED FROM PYTHON SOURCE LINES 23-32
.. code-block:: default
import numpy as np
def target_dist(x):
if (x[0]>=-5 and x[0]<=5) and (x[1]>=-5 and x[1]<=5):
return np.exp(-0.01*(x[0]-1)**2 - (x[0]**2-x[1])**2)
else:
return 0
.. GENERATED FROM PYTHON SOURCE LINES 33-34
The figure below shows how the target distribution looks like.
.. GENERATED FROM PYTHON SOURCE LINES 35-60
.. code-block:: default
import matplotlib.pyplot as plt
import itertools
x1_values = np.linspace(-5,5.1,100)
x2_values = np.linspace(-5,5.1,100)
target_values = np.zeros((100, 100))
for i, j in itertools.product(range(100), range(100)):
x1 = x1_values[i]
x2 = x2_values[j]
target_values[i,j] = target_dist(np.array([x1, x2]))
fig , ax = plt.subplots(figsize=(4,4))
ax.contourf(x1_values, x2_values, np.transpose(target_values), levels=10, cmap='Blues')
ax.set_title('(unnormalized) target distribution')
ax.set_xlabel(r'$x_1$')
ax.set_ylabel(r'$x_2$')
ax.set_xlim(-5,5)
ax.set_ylim(-5,5)
plt.tight_layout()
.. image-sg:: /auto_examples/sampler/images/sphx_glr_plot_metropolis_hastings_001.png
:alt: (unnormalized) target distribution
:srcset: /auto_examples/sampler/images/sphx_glr_plot_metropolis_hastings_001.png
:class: sphx-glr-single-img
.. GENERATED FROM PYTHON SOURCE LINES 61-67
To perform Metropolis Hastings sampling, we need to choose a proposal distribution
which can be used to determine a new state ``x'`` given a current state ``x`` at
each iteration. This is defined by the parameter ``f_sample`` which should be a
function. A usual choice is to choose the new state from a Gaussian distribution
centered at the current state.
.. GENERATED FROM PYTHON SOURCE LINES 68-75
.. code-block:: default
from scipy.stats import multivariate_normal
f_sample = multivariate_normal.rvs
# make the samples reproducible
kwgs_f_sample = {'random_state': np.random.default_rng(seed=1)}
.. GENERATED FROM PYTHON SOURCE LINES 76-77
Then, we create an instance of :class:`.MetropolisHastings` class.
.. GENERATED FROM PYTHON SOURCE LINES 78-85
.. code-block:: default
from psimpy.sampler.metropolis_hastings import MetropolisHastings
mh_sampler = MetropolisHastings(ndim=2, init_state=np.array([-4,-4]),
f_sample=f_sample, target=target_dist, nburn=1000, nthin=5,
seed=1, kwgs_f_sample=kwgs_f_sample)
.. GENERATED FROM PYTHON SOURCE LINES 86-88
Next, we call the :py:meth:`.MetropolisHastings.sample` method to draw required
number of samples.
.. GENERATED FROM PYTHON SOURCE LINES 89-95
.. code-block:: default
mh_samples, mh_accept = mh_sampler.sample(nsamples=5000)
print("Acceptance ratio: ", np.mean(mh_accept))
.. rst-class:: sphx-glr-script-out
.. code-block:: none
Acceptance ratio: 0.3712
.. GENERATED FROM PYTHON SOURCE LINES 96-97
The following figure shows how the samples look like.
.. GENERATED FROM PYTHON SOURCE LINES 98-117
.. code-block:: default
import matplotlib.pyplot as plt
import itertools
# sphinx_gallery_thumbnail_number = 2
fig , ax = plt.subplots(figsize=(4,4))
ax.contourf(x1_values, x2_values, np.transpose(target_values), levels=10, cmap='Blues')
ax.scatter(mh_samples[:,0], mh_samples[:,1], s=5, c='r', marker='o',alpha=0.05)
ax.set_title('Metropolis Hastings samples')
ax.set_xlabel(r'$x_1$')
ax.set_ylabel(r'$x_2$')
ax.set_xlim(-5,5)
ax.set_ylim(-5,5)
plt.tight_layout()
.. image-sg:: /auto_examples/sampler/images/sphx_glr_plot_metropolis_hastings_002.png
:alt: Metropolis Hastings samples
:srcset: /auto_examples/sampler/images/sphx_glr_plot_metropolis_hastings_002.png
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 0 minutes 4.769 seconds)
.. _sphx_glr_download_auto_examples_sampler_plot_metropolis_hastings.py:
.. only:: html
.. container:: sphx-glr-footer sphx-glr-footer-example
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_metropolis_hastings.py <plot_metropolis_hastings.py>`
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_metropolis_hastings.ipynb <plot_metropolis_hastings.ipynb>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/sampler/plot_saltelli.py"
.. LINE NUMBERS ARE GIVEN BELOW.
.. only:: html
.. note::
:class: sphx-glr-download-link-note
Click :ref:`here <sphx_glr_download_auto_examples_sampler_plot_saltelli.py>`
to download the full example code
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_auto_examples_sampler_plot_saltelli.py:
Saltelli sampling
=================
.. GENERATED FROM PYTHON SOURCE LINES 8-11
This example shows how to draw samples using Saltelli sampling.
Assume that there is a three-dimensional problem where X, Y, and Z are the
three random variables.
.. GENERATED FROM PYTHON SOURCE LINES 12-19
.. code-block:: default
import numpy as np
ndim = 3
# range of X is 10 to 20, range of Y is 100 to 200, range of Z is 1000 to 2000
bounds = np.array([[10,20], [100,200], [1000,2000]])
.. GENERATED FROM PYTHON SOURCE LINES 20-22
Given this setting, we can import :class:`.Saltelli`, create an instance, and
call the :py:meth:`.Saltelli.sample` method to draw required number of samples.
.. GENERATED FROM PYTHON SOURCE LINES 23-29
.. code-block:: default
from psimpy.sampler import Saltelli
saltelli_sampler = Saltelli(ndim, bounds, calc_second_order=False)
saltelli_samples = saltelli_sampler.sample(nbase=128)
.. GENERATED FROM PYTHON SOURCE LINES 30-34
In above codes, we set ``calc_second_order`` to `False`. It means that picked
samples can be used in following Sobol' analysis to compute first-order and
total-effect Sobol' indices but not second-order Sobol' indices. It leads to
:math:`nbase*(ndim+2)` samples as shown below
.. GENERATED FROM PYTHON SOURCE LINES 35-51
.. code-block:: default
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(saltelli_samples[:,0], saltelli_samples[:,1], saltelli_samples[:,2], marker='o')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.tight_layout()
print('Number of samples: ', f'{len(saltelli_samples)}')
.. image-sg:: /auto_examples/sampler/images/sphx_glr_plot_saltelli_001.png
:alt: plot saltelli
:srcset: /auto_examples/sampler/images/sphx_glr_plot_saltelli_001.png
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-script-out
.. code-block:: none
Number of samples: 640
.. GENERATED FROM PYTHON SOURCE LINES 52-55
If we want to draw samples which can also be used to compute second-order
Sobol' indices, we need to set ``calc_second_order`` to `True`.
It leads to :math:`nbase*(2*ndim+2)` samples.
.. GENERATED FROM PYTHON SOURCE LINES 56-73
.. code-block:: default
saltelli_sampler = Saltelli(ndim, bounds, calc_second_order=True)
saltelli_samples = saltelli_sampler.sample(nbase=128)
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(saltelli_samples[:,0], saltelli_samples[:,1], saltelli_samples[:,2], marker='o')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.tight_layout()
print('Number of samples: ', f'{len(saltelli_samples)}')
.. image-sg:: /auto_examples/sampler/images/sphx_glr_plot_saltelli_002.png
:alt: plot saltelli
:srcset: /auto_examples/sampler/images/sphx_glr_plot_saltelli_002.png
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-script-out
.. code-block:: none
Number of samples: 1024
.. GENERATED FROM PYTHON SOURCE LINES 74-78
.. note:: If one has a two-dimensional problem, there is no need to set
``calc_second_order`` to `True`. The reason is that the second-order Sobol'
index can be directly computed based on the first-order and total-effect
Sobol' index in that case.
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 0 minutes 0.309 seconds)
.. _sphx_glr_download_auto_examples_sampler_plot_saltelli.py:
.. only:: html
.. container:: sphx-glr-footer sphx-glr-footer-example
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_saltelli.py <plot_saltelli.py>`
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_saltelli.ipynb <plot_saltelli.ipynb>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
:orphan:
.. _sphx_glr_auto_examples_sampler_sg_execution_times:
Computation times
=================
**00:05.348** total execution time for **auto_examples_sampler** files:
+-----------------------------------------------------------------------------------------------------+-----------+--------+
| :ref:`sphx_glr_auto_examples_sampler_plot_metropolis_hastings.py` (``plot_metropolis_hastings.py``) | 00:04.769 | 0.0 MB |
+-----------------------------------------------------------------------------------------------------+-----------+--------+
| :ref:`sphx_glr_auto_examples_sampler_plot_saltelli.py` (``plot_saltelli.py``) | 00:00.309 | 0.0 MB |
+-----------------------------------------------------------------------------------------------------+-----------+--------+
| :ref:`sphx_glr_auto_examples_sampler_plot_latin.py` (``plot_latin.py``) | 00:00.270 | 0.0 MB |
+-----------------------------------------------------------------------------------------------------+-----------+--------+
:orphan:
Sensitivity
===========
.. raw:: html
<div class="sphx-glr-thumbnails">
.. raw:: html
<div class="sphx-glr-thumbcontainer" tooltip="Sobol&#x27; analysis">
.. only:: html
.. image:: /auto_examples/sensitivity/images/thumb/sphx_glr_plot_sobol_analyze_thumb.png
:alt: Sobol' analysis
:ref:`sphx_glr_auto_examples_sensitivity_plot_sobol_analyze.py`
.. raw:: html
<div class="sphx-glr-thumbnail-title">Sobol' analysis</div>
</div>
.. raw:: html
</div>
.. toctree::
:hidden:
/auto_examples/sensitivity/plot_sobol_analyze
.. only:: html
.. container:: sphx-glr-footer sphx-glr-footer-gallery
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download all examples in Python source code: sensitivity_python.zip </auto_examples/sensitivity/sensitivity_python.zip>`
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download all examples in Jupyter notebooks: sensitivity_jupyter.zip </auto_examples/sensitivity/sensitivity_jupyter.zip>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/sensitivity/plot_sobol_analyze.py"
.. LINE NUMBERS ARE GIVEN BELOW.
.. only:: html
.. note::
:class: sphx-glr-download-link-note
Click :ref:`here <sphx_glr_download_auto_examples_sensitivity_plot_sobol_analyze.py>`
to download the full example code
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_auto_examples_sensitivity_plot_sobol_analyze.py:
Sobol' analysis
===============
.. GENERATED FROM PYTHON SOURCE LINES 8-16
This example shows how to compute Sobol' sensitivity indices using
the :class:`.SobolAnalyze` class.
Let us consider the Ishigami function :cite:p:`LeGratiet2014`
:math:`y=f\left(x_1, x_2, x_3\right)=\sin \left(x_1\right)+7 \sin \left(\mathrm{x}_2\right)^2+0.1 x_3^4 \sin \left(x_1\right)`
where :math:`x_i` is uniformly distributed on :math:`[-\pi, \pi]` for :math:`i \in [1,2,3]`. We define the function as follows:
.. GENERATED FROM PYTHON SOURCE LINES 17-24
.. code-block:: default
import numpy as np
def f(x1,x2,x3):
return np.sin(x1) + 7*np.sin(x2)**2 + 0.1*x3**4*np.sin(x1)
.. GENERATED FROM PYTHON SOURCE LINES 25-32
The variation of each :math:`x_i` leads to the variation of :math:`y`. Sobol' analysis is used to quantify the contribution
of each :math:`x_i` and their interactions to the variation of :math:`y`. Steps of a Sobol' analysis are:
#. Draw samples of input.
#. Evaluate the model at each input sample to obtain model output.
#. Compute Sobol' indices based on model outputs.
.. GENERATED FROM PYTHON SOURCE LINES 33-57
.. code-block:: default
from psimpy.sampler import Saltelli
from psimpy.sensitivity import SobolAnalyze
ndim = 3
bounds = np.array([[-np.pi, np.pi], [-np.pi, np.pi], [-np.pi, np.pi]])
calc_second_order = False
nbase = 1024
# draw Saltelli samples
saltelli_sampler = Saltelli(ndim, bounds, calc_second_order)
saltelli_samples = saltelli_sampler.sample(nbase)
# Evaluate the model to obtain model outputs
Y = f(saltelli_samples[:,0], saltelli_samples[:,1], saltelli_samples[:, 2])
# Compute Sobol' indices
sobol_analyzer = SobolAnalyze(ndim, Y, calc_second_order, seed=2)
S_res = sobol_analyzer.run()
print("Estimated first-order Sobol' index: ", S_res['S1'][:,0])
print("Analytical solution: [0.314, 0.442, 0]")
.. rst-class:: sphx-glr-script-out
.. code-block:: none
Estimated first-order Sobol' index: [0.31430871 0.44545644 0.01325799]
Analytical solution: [0.314, 0.442, 0]
.. GENERATED FROM PYTHON SOURCE LINES 58-62
Figure below shows the results. The bar plots represent estimated first-order
and total-effect Sobol' indices, as well as their 95% confidence interval. The
red stars represent true values of the first-order Sobol' indices.
.. GENERATED FROM PYTHON SOURCE LINES 63-103
.. code-block:: default
import matplotlib.pyplot as plt
name_vars = [r'$x_1$',r'$x_2$',r'$x_3$']
x = {'S1':[1,4,7], 'ST':[2,5,8]}
color = {'S1':'#69b3a2', 'ST':'#3399e6'}
width=0.8
capsize=3
fig, ax1 = plt.subplots(figsize=(4,3))
# estimated first-order Sobol' indices and their confidence interval
ax1.bar(x['S1'], S_res['S1'][:,0], width=width, yerr=S_res['S1'][:,2], color=color['S1'], capsize=capsize)
ax1.axis([0,9,0,0.7])
ax1.set_xticks([1.5,4.5,7.5])
ax1.set_xticklabels(name_vars)
ax1.tick_params(axis='both',bottom=False)
ax1.tick_params(axis="y", labelcolor=color['S1'])
ax1.yaxis.grid(True, linestyle='--', which='major', color='grey', alpha=.25)
ax1.set_ylabel('First-order index', color=color['S1'])
# analytical solution for first-order indices
ax1.scatter(x['S1'], [0.314, 0.442, 0], s=30, c='red', marker='*')
# estimated total-effect Sobol' indices and their confidence interval
ax2 = ax1.twinx()
ax2.bar(x['ST'], S_res['ST'][:,0], width=width, yerr=S_res['ST'][:,2], color=color['ST'], capsize=capsize)
ax2.axis([0,9,0,0.7])
ax2.set_xticks([1.5,4.5,7.5])
ax2.tick_params(axis='both',bottom=False)
ax2.tick_params(axis="y", labelcolor=color['ST'])
ax2.yaxis.grid(True, linestyle='--', which='major', color='grey', alpha=.25)
ax2.set_ylabel('Total-effect index', color=color['ST'])
plt.tight_layout()
.. image-sg:: /auto_examples/sensitivity/images/sphx_glr_plot_sobol_analyze_001.png
:alt: plot sobol analyze
:srcset: /auto_examples/sensitivity/images/sphx_glr_plot_sobol_analyze_001.png
:class: sphx-glr-single-img
.. GENERATED FROM PYTHON SOURCE LINES 104-109
.. note:: The interaction between one parameter and all other parameters can be
quantified by the difference between its first-order and total-effet Sobol'
indices. If detailed information on interaction between any two parameters
is needed, one can set ``calc_second_order`` to `True`.
.. GENERATED FROM PYTHON SOURCE LINES 112-128
**Many models of interest in real-world applications are very computationally
expensive. In that case, we resort to emulation techniques in order to efficiently
conduct Sobol' analysis.** See :cite:t:`Zhao2021a` for the detailed theory.
Steps of such an analysis are:
#. Draw input samples for emulator training.
#. Evaluate the model at each training input point to obtain training data.
#. Train an emulator for the computationally expensive model of interest.
#. Draw samples of input for Sobol' analysis.
#. Evaluate the emulator at each input sample at step 3 to obtain approximated model outputs.
#. Compute Sobol' indices based on approximated model outputs.
Assume that the above Ishigami function is computationally expensive, an emulator-based
Sobol' analysis can be conducted as follows:
.. GENERATED FROM PYTHON SOURCE LINES 129-178
.. code-block:: default
import numpy as np
from psimpy.sampler import LHS
from psimpy.emulator import ScalarGaSP
from psimpy.sampler import Saltelli
from psimpy.sensitivity import SobolAnalyze
# define the model
def f(x1,x2,x3):
return np.sin(x1) + 7*np.sin(x2)**2 + 0.1*x3**4*np.sin(x1)
ndim = 3
bounds = np.array([[-np.pi, np.pi], [-np.pi, np.pi], [-np.pi, np.pi]])
seed = 2
# draw 100 input samples for emulator training using Latin hypercube sampling
lhs_sampler = LHS(ndim, bounds, seed, criterion='maximin', iteration=200)
design = lhs_sampler.sample(nsamples=100)
# evaluate the model at each training input point
response = f(design[:,0], design[:,1], design[:,2])
# train an emulator
scalar_gasp = ScalarGaSP(ndim)
scalar_gasp.train(design, response)
# the trained emulator needs to be validated before moving on
# here we leave out this step for simplicity...
# draw samples of input for Sobol' analysis
calc_second_order = False
nbase = 1024
saltelli_sampler = Saltelli(ndim, bounds, calc_second_order)
saltelli_samples = saltelli_sampler.sample(nbase)
# evaluate the emulator to obtain approximated model output.
# in order to take emulator-induced uncertainty into account, we draw 50 realizations
# from the emulator
Y = scalar_gasp.sample(saltelli_samples, 50)
# Compute Sobol' indices based on approximated model outputs
sobol_analyzer = SobolAnalyze(ndim, Y, calc_second_order, seed=seed)
S_res = sobol_analyzer.run(mode='parallel', max_workers=5)
print("\n")
print("Emulator-based first-order Sobol' index: ", S_res['S1'][:,0])
print("Analytical solution: [0.314, 0.442, 0]")
.. rst-class:: sphx-glr-script-out
.. code-block:: none
The upper bounds of the range parameters are 173,6849 174,6991 174,3157
The initial values of range parameters are 3,473697 3,493981 3,486314
Start of the optimization 1 :
The number of iterations is 14
The value of the marginal posterior function is -284,0714
Optimized range parameters are 3,183999 2,786247 4,315419
Optimized nugget parameter is 0
Convergence: TRUE
The initial values of range parameters are 1,721972 1,732027 1,728226
Start of the optimization 2 :
The number of iterations is 14
The value of the marginal posterior function is -284,0714
Optimized range parameters are 3,183999 2,786247 4,315419
Optimized nugget parameter is 0
Convergence: TRUE
Emulator-based first-order Sobol' index: [0.29064086 0.46853437 0.01096656]
Analytical solution: [0.314, 0.442, 0]
.. GENERATED FROM PYTHON SOURCE LINES 179-182
The following figure shows how the results of emulator-based Sobol' analysis
look like.
.. GENERATED FROM PYTHON SOURCE LINES 182-213
.. code-block:: default
import matplotlib.pyplot as plt
name_vars = [r'$x_1$',r'$x_2$',r'$x_3$']
x = {'S1':[1,4,7], 'ST':[2,5,8]}
color = {'S1':'#69b3a2', 'ST':'#3399e6'}
width=0.8
capsize=3
fig, ax1 = plt.subplots(figsize=(4,3))
ax1.bar(x['S1'], S_res['S1'][:,0], width=width, yerr=S_res['S1'][:,2], color=color['S1'], capsize=capsize)
ax1.axis([0,9,0,0.7])
ax1.set_xticks([1.5,4.5,7.5])
ax1.set_xticklabels(name_vars)
ax1.tick_params(axis='both',bottom=False)
ax1.tick_params(axis="y", labelcolor=color['S1'])
ax1.yaxis.grid(True, linestyle='--', which='major', color='grey', alpha=.25)
ax1.set_ylabel('First-order index', color=color['S1'])
ax1.scatter(x['S1'], [0.314, 0.442, 0], s=30, c='red', marker='*')
ax2 = ax1.twinx()
ax2.bar(x['ST'], S_res['ST'][:,0], width=width, yerr=S_res['ST'][:,2], color=color['ST'], capsize=capsize)
ax2.axis([0,9,0,0.7])
ax2.set_xticks([1.5,4.5,7.5])
ax2.tick_params(axis='both',bottom=False)
ax2.tick_params(axis="y", labelcolor=color['ST'])
ax2.yaxis.grid(True, linestyle='--', which='major', color='grey', alpha=.25)
ax2.set_ylabel('Total-effect index', color=color['ST'])
plt.tight_layout()
.. image-sg:: /auto_examples/sensitivity/images/sphx_glr_plot_sobol_analyze_002.png
:alt: plot sobol analyze
:srcset: /auto_examples/sensitivity/images/sphx_glr_plot_sobol_analyze_002.png
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 0 minutes 16.499 seconds)
.. _sphx_glr_download_auto_examples_sensitivity_plot_sobol_analyze.py:
.. only:: html
.. container:: sphx-glr-footer sphx-glr-footer-example
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_sobol_analyze.py <plot_sobol_analyze.py>`
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_sobol_analyze.ipynb <plot_sobol_analyze.ipynb>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
:orphan:
.. _sphx_glr_auto_examples_sensitivity_sg_execution_times:
Computation times
=================
**00:16.499** total execution time for **auto_examples_sensitivity** files:
+---------------------------------------------------------------------------------------------+-----------+--------+
| :ref:`sphx_glr_auto_examples_sensitivity_plot_sobol_analyze.py` (``plot_sobol_analyze.py``) | 00:16.499 | 0.0 MB |
+---------------------------------------------------------------------------------------------+-----------+--------+
:orphan:
Simulator
=========
.. raw:: html
<div class="sphx-glr-thumbnails">
.. raw:: html
<div class="sphx-glr-thumbcontainer" tooltip="Mass Point Model">
.. only:: html
.. image:: /auto_examples/simulator/images/thumb/sphx_glr_plot_mass_point_model_thumb.png
:alt: Mass Point Model
:ref:`sphx_glr_auto_examples_simulator_plot_mass_point_model.py`
.. raw:: html
<div class="sphx-glr-thumbnail-title">Mass Point Model</div>
</div>
.. raw:: html
<div class="sphx-glr-thumbcontainer" tooltip="RunSimulator: Mass Point Model">
.. only:: html
.. image:: /auto_examples/simulator/images/thumb/sphx_glr_plot_run_mass_point_model_thumb.png
:alt: RunSimulator: Mass Point Model
:ref:`sphx_glr_auto_examples_simulator_plot_run_mass_point_model.py`
.. raw:: html
<div class="sphx-glr-thumbnail-title">RunSimulator: Mass Point Model</div>
</div>
.. raw:: html
<div class="sphx-glr-thumbcontainer" tooltip="Ravaflow24 Mixture Model">
.. only:: html
.. image:: /auto_examples/simulator/images/thumb/sphx_glr_plot_ravaflow24_thumb.png
:alt: Ravaflow24 Mixture Model
:ref:`sphx_glr_auto_examples_simulator_plot_ravaflow24.py`
.. raw:: html
<div class="sphx-glr-thumbnail-title">Ravaflow24 Mixture Model</div>
</div>
.. raw:: html
<div class="sphx-glr-thumbcontainer" tooltip="RunSimulator: Ravaflow24 Mixture Model">
.. only:: html
.. image:: /auto_examples/simulator/images/thumb/sphx_glr_plot_run_ravaflow24_thumb.png
:alt: RunSimulator: Ravaflow24 Mixture Model
:ref:`sphx_glr_auto_examples_simulator_plot_run_ravaflow24.py`
.. raw:: html
<div class="sphx-glr-thumbnail-title">RunSimulator: Ravaflow24 Mixture Model</div>
</div>
.. raw:: html
</div>
.. toctree::
:hidden:
/auto_examples/simulator/plot_mass_point_model
/auto_examples/simulator/plot_run_mass_point_model
/auto_examples/simulator/plot_ravaflow24
/auto_examples/simulator/plot_run_ravaflow24
.. only:: html
.. container:: sphx-glr-footer sphx-glr-footer-gallery
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download all examples in Python source code: simulator_python.zip </auto_examples/simulator/simulator_python.zip>`
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download all examples in Jupyter notebooks: simulator_jupyter.zip </auto_examples/simulator/simulator_jupyter.zip>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment