diff --git a/app/tutorials/Tutorial 1 - Function Approximation.ipynb b/app/tutorials/Tutorial 1 - Function Approximation.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..41a434f56095fb0df7d30de5543c2718be2492f1
--- /dev/null
+++ b/app/tutorials/Tutorial 1 - Function Approximation.ipynb
@@ -0,0 +1,213 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "00c2cc0d-24ab-4a35-bc11-459feebe2455",
+ "metadata": {},
+ "source": [
+ "# Tutorial 1 - Function Approximation with PyThia"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "7bd31586-38d1-4e41-8963-a0e02a070f7e",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import sys\n",
+ "import numpy as np\n",
+ "import pythia as pt # see doc: www.pythia-uq.com\n",
+ "import matplotlib.pyplot as plt\n",
+ "\n",
+ "sys.path.append(\"../../src\")\n",
+ "from utils import Batman # local imports"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "31a019e4-dcc1-47d9-be36-24b95e24abe9",
+ "metadata": {},
+ "source": [
+ "## Define Target Function"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ceb41891-6c8e-40ed-b55e-979e00713ea9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def target_function(y):\n",
+ " return Batman().eval(y)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d4a36225-6620-4919-95fa-cffe3446fef6",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "y = np.linspace(-7.5, 7.5, 401)\n",
+ "fig, ax = plt.subplot_mosaic([[\"batman\"]], constrained_layout=True, figsize=(10, 5))\n",
+ "ax[\"batman\"].plot(y, target_function(y)[:, 0], label=\"$f_1(y)$\")\n",
+ "ax[\"batman\"].plot(y, target_function(y)[:, 1], label=\"$f_2(y)$\")\n",
+ "ax[\"batman\"].set_xlabel(\"y\")\n",
+ "ax[\"batman\"].legend()\n",
+ "plt.show();"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "fde90e4d-be77-4b89-ae03-9ede7faa2447",
+ "metadata": {},
+ "source": [
+ "## Setup Polynomial Chaos Expansion"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b3f735c8-debd-42ba-b642-d4fb0f1ad62e",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "param1 = pt.parameter.Parameter(index=1, name=\"y_1\", domain=[-7.5, 7.5], distribution='uniform')\n",
+ "params = [param1]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "67d2abd1-2ceb-442e-ae7e-417890a86e4b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "sdim = [75]\n",
+ "indices = pt.index.tensor(sdim)\n",
+ "index_set = pt.index.IndexSet(indices)\n",
+ "print(f\"Multiindex information:\")\n",
+ "print(f\" number of indices: {index_set.shape[0]}\")\n",
+ "print(f\" dimension: {index_set.shape[1]}\")\n",
+ "print(f\" maximum dimension: {index_set.max}\")\n",
+ "print(f\" number of sobol indices: {len(index_set.sobol_tuples)}\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "28be82de-2de1-45e6-8d83-7b61175b2e87",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "s = pt.sampler.ParameterSampler(params)\n",
+ "y_train = s.sample(1000)\n",
+ "w_train = s.weight(y_train)\n",
+ "f_train = target_function(y_train)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "59364c66-9409-4214-8c35-4e7ed1afaa72",
+ "metadata": {},
+ "source": [
+ "## Run Polynomial Chaos Approximation"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5bea2800-edb2-4941-b238-d94825b0db23",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "surrogate = pt.chaos.PolynomialChaos(\n",
+ " params, index_set, y_train, w_train, f_train)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "34011607-9de1-4ce0-9bd7-6b843a662703",
+ "metadata": {},
+ "source": [
+ "## Compute Approximation Error"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "269991dd-bbca-4a64-a9e6-89fe2ee0f845",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "test_sampler = pt.sampler.ParameterSampler(params)\n",
+ "y_test = test_sampler.sample(1000)\n",
+ "f_test = target_function(y_test)\n",
+ "f_approx = surrogate.eval(y_test)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "baa189c4-55d9-4087-8330-7de826406097",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "y = np.linspace(-7.1, 7.1, 200).reshape(-1, 1)\n",
+ "f_target = target_function(y)\n",
+ "f_approx = surrogate.eval(y)\n",
+ "f_error = np.abs(f_target-f_approx)/np.abs(f_target)\n",
+ "\n",
+ "fig, ax = plt.subplot_mosaic([[\"target\", \"target\"],\n",
+ " [\"approx\", \"error\"]],\n",
+ " constrained_layout=True, figsize=(10, 4))\n",
+ "fig.suptitle(\"Approximation of function with PyThia\")\n",
+ "\n",
+ "ax[\"target\"].scatter(y_train[:,0], f_train[:,0], s=1)\n",
+ "ax[\"target\"].scatter(y_train[:,0], f_train[:,1], s=1)\n",
+ "ax[\"target\"].set_title(\"training data points\")\n",
+ "\n",
+ "ax[\"approx\"].plot(y, f_approx[:,0])\n",
+ "ax[\"approx\"].plot(y, f_approx[:,1])\n",
+ "ax[\"approx\"].set_title(\"surrogate\")\n",
+ "\n",
+ "ax[\"error\"].semilogy(y, f_error[:,0])\n",
+ "ax[\"error\"].semilogy(y, f_error[:,1])\n",
+ "ax[\"error\"].set_title(\"pointwise relaive error\")\n",
+ "\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5d1f5311-1a1e-42c7-a20e-2f8a274ce434",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.9.16"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/app/tutorials/Tutorial 2 - Sensitivity Analysis.ipynb b/app/tutorials/Tutorial 2 - Sensitivity Analysis.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..caf26096bf57009020c48bb72df6b25d91e8cc3e
--- /dev/null
+++ b/app/tutorials/Tutorial 2 - Sensitivity Analysis.ipynb
@@ -0,0 +1,210 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "ceb7573d-54c3-417a-b533-4a87969adc59",
+ "metadata": {},
+ "source": [
+ "# Tutorial 2 - Sensitivity Analysis with Sobol Indices"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "42d7cee5-4691-41c5-83d2-4ce887966690",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import sys\n",
+ "import numpy as np\n",
+ "import pythia as pt # see doc: www.pythia-uq.com\n",
+ "\n",
+ "sys.path.append(\"../../src\")\n",
+ "from utils import Sobol # local imports"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "462aae79-a70a-4a9c-90c7-3d73b71dc287",
+ "metadata": {},
+ "source": [
+ "## Define Target Function\n",
+ "\n",
+ "For $y\\sim\\mathcal{U}([0,1]^M)$ define\n",
+ "\n",
+ "$$\n",
+ "f(y) = \\prod_{m=1}^M\\frac{\\vert 4y_m-2\\vert + a_m}{1+a_m}\n",
+ "\\qquad\\mbox{for }a_1,\\dots,a_M >= 0.\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "2a2462f7-4abb-4468-89ea-a33e8268209f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "a = np.array([1, 2, 3, 4])\n",
+ "sobol = Sobol(a)\n",
+ "def target_function(y):\n",
+ " return sobol.eval(y)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ce631590-d241-40a9-8604-9b4035931a31",
+ "metadata": {},
+ "source": [
+ "## Compute Polynomial Chaos Approximation"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "7004ecb6-18b9-4f36-96a7-b36e08a09a2d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# define parameters\n",
+ "params = [pt.parameter.Parameter(index=j, name=f\"y_{j+1}\", domain=[0, 1], distribution='uniform')\n",
+ " for j in range(sobol.a.size)]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4911ed05-8f09-4eb9-8880-7c05e49e0f93",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# set index set\n",
+ "indices = pt.index.simplex(len(params), 6) # limit total polynomial degree to 10\n",
+ "index_set = pt.index.IndexSet(indices)\n",
+ "print(f\"multiindex information:\")\n",
+ "print(f\" number of indices: {index_set.shape[0]}\")\n",
+ "print(f\" dimension: {index_set.shape[1]}\")\n",
+ "print(f\" number of sobol indices: {len(index_set.sobol_tuples)}\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "35c9b837-e3ee-436c-b7d7-c3fca55bf389",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# generate training data\n",
+ "s = pt.sampler.ParameterSampler(params)\n",
+ "y_train = s.sample(10_000)\n",
+ "w_train = s.weight(y_train)\n",
+ "f_train = target_function(y_train)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "212a4b5c-fff6-4d36-bb3f-c60d0304a3c2",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# compute PC expansion\n",
+ "surrogate = pt.chaos.PolynomialChaos(params, index_set, y_train, w_train, f_train)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "3c4429d6-90e6-4cf0-bb73-24e50c198f15",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# test approximation\n",
+ "s_test = pt.sampler.ParameterSampler(params)\n",
+ "y_test = s_test.sample(1_000)\n",
+ "f_test = target_function(y_test)\n",
+ "f_approx = surrogate.eval(y_test)\n",
+ "\n",
+ "error_L2 = np.sqrt(np.sum((f_test-f_approx)**2)/y_test.shape[0])\n",
+ "error_L2_rel = error_L2 / np.sqrt(np.sum((f_test)**2)/y_test.shape[0])\n",
+ "error_max = np.max(np.abs(f_test-f_approx))\n",
+ "error_max_rel = np.max(np.abs(f_test-f_approx)/np.abs(f_test))\n",
+ "\n",
+ "print(f\"test error L2 (abs/rel):\",\n",
+ " f\" {error_L2:4.2e} / {error_L2_rel:4.2e}\")\n",
+ "print(f\"test error max (abs/rel):\",\n",
+ " f\" {error_max:4.2e} / {error_max_rel:4.2e}\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "bdbebf3d-d810-401f-984a-79763d8083e0",
+ "metadata": {},
+ "source": [
+ "## Compare Sobol Indices\n",
+ "\n",
+ "For the Sobol function we have $\\mathbb{E}[f]=1$ and\n",
+ "\n",
+ "$$\n",
+ " \\operatorname{Var}[f] = \\prod_{m=1}^M (1+c_m)-1\n",
+ " \\qquad\\mbox{for}\\qquad\n",
+ " c_m = \\frac{1}{3(1+a_m)^2}.\n",
+ "$$\n",
+ "\n",
+ "The Sobol indices are thus given by\n",
+ "\n",
+ "$$\n",
+ " \\operatorname{Sob}_{m_1,\\dots,m_s}\n",
+ " = \\frac{1}{\\operatorname{Var}[f]} \\prod_{k=1}^s c_{m_k}.\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "3c736cef-fe28-4ef6-b32a-e3da774d483d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# compare Sobol indices\n",
+ "tuples, indices = sobol.indices()\n",
+ "approx = surrogate.sobol.T[0]\n",
+ "difference = np.abs(indices-approx)\n",
+ "\n",
+ "print(f\"{'tuple':<15} {'exact':<10} {'approx':<10} {'difference':<9}\")\n",
+ "print(\"-\"*48)\n",
+ "for tup, ind, app, diff in zip(tuples, indices, approx, difference):\n",
+ " print(f\"{str(tup):<14} {ind:4.2e} {app:4.2e} {diff:4.2e}\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a99601f6-52a1-4709-94f6-0aeb2136689a",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.9.16"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/app/tutorials/Tutorial 3 - MCMC.ipynb b/app/tutorials/Tutorial 3 - MCMC.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..8e232907aa7e298cdaded065f9e77abff47665ab
--- /dev/null
+++ b/app/tutorials/Tutorial 3 - MCMC.ipynb
@@ -0,0 +1,261 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "39584f41-1647-4f93-a3b5-84a39b3e13f9",
+ "metadata": {},
+ "source": [
+ "# Tutorial 3 - Markov-Chain Monte Carlo"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "914c9beb-b5de-4a26-a709-3005371211c9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import sys\n",
+ "import numpy as np\n",
+ "import pythia as pt # see doc: www.pythia-uq.com\n",
+ "import emcee\n",
+ "import matplotlib.pyplot as plt\n",
+ "import matplotlib.cm as cmap\n",
+ "\n",
+ "sys.path.append(\"../../src\")\n",
+ "from utils import Sine, Gaussian # local imports"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "43851283-9de8-42c3-90bc-1020638bdbb2",
+ "metadata": {},
+ "source": [
+ "## Define Target Function\n",
+ "For $y\\sim\\mathcal{U}([0,1]^6)$ define\n",
+ "$$\n",
+ "f(y) = (f_1(y), \\dots, f_{200}(y)) \\in\\mathbb{R}^{200}\n",
+ "\\qquad\\mbox{with}\\qquad\n",
+ "f_j(y) = \\sum_{m=1}^6 e^{-m}y_m \\sin\\Bigl(\\frac{2\\pi(j-1)}{199}\\Bigr).\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "31fdf19b-020e-41e4-b4ab-71c9502271c6",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "sine = Sine(np.exp(-np.arange(1,7)))\n",
+ "def target_function(y):\n",
+ " return sine.eval(y)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a19c32ae-3860-4194-9a13-21741a1da6db",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "y = np.random.uniform(0, 1, (10, 6))\n",
+ "val = target_function(y)\n",
+ "fig, ax = plt.subplot_mosaic([[\"sine\"]], constrained_layout=True, figsize=(10, 5))\n",
+ "for j, valj in enumerate(val):\n",
+ " ax[\"sine\"].plot(sine.points, valj, color=cmap.tab20(j), label=f\"sample {j+1}\")\n",
+ " ax[\"sine\"].legend()\n",
+ "plt.show();"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c7d59cb5-3a68-4750-bf18-9220f06a224c",
+ "metadata": {},
+ "source": [
+ "## Set Measurement"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9507a69e-9811-4be4-b15e-59530989018c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# measurement\n",
+ "y_true = np.random.uniform(0.2, 0.8, (1, sine.scales.size))\n",
+ "f_true = target_function(y_true).reshape(1, -1)\n",
+ "c = 0.03\n",
+ "noise = np.random.normal(0, c, f_true.shape)\n",
+ "f_meas = f_true + noise"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ee946804-feda-4a90-93af-e702a267c86f",
+ "metadata": {},
+ "source": [
+ "## Define Bayesian Posterior"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "59f218b2-2470-4672-9e69-de398956c776",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "domain = np.array([[0, 0, 0, 0, 0, 0, 0.01],\n",
+ " [1, 1, 1, 1, 1, 1, 0.5]]).T\n",
+ "prior = pt.sampler.ProductSampler([pt.sampler.UniformSampler(dom) for dom in domain]) # extended prior\n",
+ "forward_model = lambda y: target_function(y[:,:-1]) # extended forward model\n",
+ "sigma = lambda y: np.abs(y[:, -1]).reshape(-1, 1)*np.ones((1, f_meas.shape[1])) # error model\n",
+ "\n",
+ "# define Bayesian posterior\n",
+ "lh = Gaussian(forward_model, sigma, xdim=domain.shape[0])\n",
+ "def log_posterior(y, meas):\n",
+ " return lh.log_likelihood(y, meas) + prior.log_pdf(y)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "37d82fb1-9de9-4e4e-9f33-b9d10d0f2e76",
+ "metadata": {},
+ "source": [
+ "## Setup MCMC"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9e867a60-7337-4f24-ab14-e28f60285376",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "n_walkers = 15\n",
+ "chain_length = 2000\n",
+ "seed = prior.sample(n_walkers)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "fa098095-592e-4b76-ba86-8515fb9b0c76",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# define runner\n",
+ "runner = emcee.EnsembleSampler(n_walkers, domain.shape[0], log_posterior, args=[f_meas])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "76308649-ecad-4a4c-ade5-8074f7673e32",
+ "metadata": {},
+ "source": [
+ "## Run MCMC"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "06aedd14-aadf-4f14-a6b3-3a32c0286286",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# discard burn-in\n",
+ "state = runner.run_mcmc(seed, 1000, progress=True)\n",
+ "runner.reset()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "2ab5096c-8e4a-4dc3-8edd-02a543c495ce",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# run MCMC\n",
+ "runner.run_mcmc(seed, chain_length, progress=True);"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5b2684a6-1a92-4820-a5da-d9259a074d49",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# simple evaluation\n",
+ "samples = runner.get_chain(flat=True)\n",
+ "mean = np.mean(samples, axis=0)\n",
+ "std = np.std(samples, axis=0)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "14dfbf69-c251-43dd-bfa2-5588b0bcc3cc",
+ "metadata": {},
+ "source": [
+ "## Data Analysis"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b84800c1-350a-4581-975c-58e238c83b91",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# plot posterior distribution\n",
+ "layout = [[0, 1], [2, 3], [4, 5], [6, 6]]\n",
+ "fig, axes = plt.subplot_mosaic(layout, figsize=(12, 7.5))\n",
+ "for j, (dom, ax) in enumerate(zip(domain, axes.values())):\n",
+ " bins = 100\n",
+ " ax.hist(samples[:, j], bins=bins, density=True, range=dom, alpha=0.5)\n",
+ " ax.axvline(x=mean[j], color=cmap.tab10(0), label=\"mean\")\n",
+ " ax.axvline(x=mean[j]-std[j], ls=\"--\", color=cmap.tab10(0), label=\"std\")\n",
+ " ax.axvline(x=mean[j]+std[j], ls=\"--\", color=cmap.tab10(0))\n",
+ " if j == 6:\n",
+ " ax.axvline(x=c, color=\"k\", label=\"true value\")\n",
+ " ax.set_ylabel(f\"$c$\")\n",
+ " else: \n",
+ " ax.axvline(x=y_true[0, j], color=\"k\", label=\"true value\")\n",
+ " ax.set_ylabel(f\"$y_{j+1}$\")\n",
+ " ax.set_yticks([])\n",
+ " ax.legend()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9f3e4661-663e-41ab-beb4-5734a11a4d94",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.9.16"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/app/tutorials/Tutorial 4 - Sampling.ipynb b/app/tutorials/Tutorial 4 - Sampling.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..5214b5e9d830fa63c0a2a03536d608297c346d72
--- /dev/null
+++ b/app/tutorials/Tutorial 4 - Sampling.ipynb
@@ -0,0 +1,240 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "b304f6f3-ad5b-4123-9560-316b6d5cb434",
+ "metadata": {},
+ "source": [
+ "# Tutorial 4 - Sampling"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "78eb2cc1-33e6-4fe3-bb47-70e6aaca324a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import pythia as pt # check doc: www.pythia-uq.com\n",
+ "import matplotlib.pyplot as plt\n",
+ "import matplotlib.cm as cmap"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "9b31d394-d268-4c7e-af46-bff7c38d71e2",
+ "metadata": {
+ "tags": []
+ },
+ "source": [
+ "## Univariate Sampling"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "69a027e1-e887-47ea-882b-04d2781b4a93",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# uniform\n",
+ "s = pt.sampler.UniformSampler([-1,2])\n",
+ "samples = s.sample(100_000)\n",
+ "y = np.linspace(-1.5, 2.5, 200)\n",
+ "\n",
+ "fig, ax = plt.subplot_mosaic([[\"fig\"]], constrained_layout=True)\n",
+ "ax[\"fig\"].hist(samples.flatten(), bins=50, density=True, label=\"samples\")\n",
+ "ax[\"fig\"].plot(y, s.pdf(y), label=\"pdf\")\n",
+ "ax[\"fig\"].legend()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "bda40e48-deb9-4095-a25c-f5fd3b7dd5b5",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# normal\n",
+ "s = pt.sampler.NormalSampler(mean=1, var=1)\n",
+ "samples = s.sample(100_000)\n",
+ "y = np.linspace(-3.5, 5.5, 200)\n",
+ "\n",
+ "fig, ax = plt.subplot_mosaic([[\"fig\"]], constrained_layout=True)\n",
+ "ax[\"fig\"].hist(samples.flatten(), bins=50, density=True, label=\"samples\")\n",
+ "ax[\"fig\"].plot(y, s.pdf(y), label=\"pdf\")\n",
+ "ax[\"fig\"].legend()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "c0dea9d0-1c87-43c7-be74-cebba7efe6c1",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Gamma\n",
+ "s1 = pt.sampler.GammaSampler(domain=[0, np.inf], alpha=1, beta=1)\n",
+ "samples1 = s1.sample(100_000)\n",
+ "s2 = pt.sampler.GammaSampler(domain=[0, np.inf], alpha=5, beta=5)\n",
+ "samples2 = s2.sample(100_000)\n",
+ "y = np.linspace(-0.5, 5.5, 200)\n",
+ "\n",
+ "fig, ax = plt.subplot_mosaic([[\"a\", \"b\"]], constrained_layout=True)\n",
+ "ax[\"a\"].hist(samples1.flatten(), bins=50, density=True, label=\"samples\")\n",
+ "ax[\"a\"].plot(y, s1.pdf(y), label=\"pdf\")\n",
+ "ax[\"a\"].legend()\n",
+ "ax[\"b\"].hist(samples2.flatten(), bins=50, density=True, label=\"samples\")\n",
+ "ax[\"b\"].plot(y, s2.pdf(y), label=\"pdf\")\n",
+ "ax[\"b\"].legend()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "8e9ab462-81d0-4f9c-9fa9-21884ba1a9b1",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Beta\n",
+ "s1 = pt.sampler.BetaSampler(domain=[0, 3], alpha=0.8, beta=0.8)\n",
+ "samples1 = s1.sample(100_000)\n",
+ "s2 = pt.sampler.BetaSampler(domain=[0, 3], alpha=10, beta=2)\n",
+ "samples2 = s2.sample(100_000)\n",
+ "y = np.linspace(-0.5, 3.5, 200)\n",
+ "\n",
+ "fig, ax = plt.subplot_mosaic([[\"a\", \"b\"]], constrained_layout=True)\n",
+ "ax[\"a\"].hist(samples1.flatten(), bins=50, density=True, label=\"samples\")\n",
+ "ax[\"a\"].plot(y, s1.pdf(y), label=\"pdf\")\n",
+ "ax[\"a\"].legend()\n",
+ "ax[\"b\"].hist(samples2.flatten(), bins=50, density=True, label=\"samples\")\n",
+ "ax[\"b\"].plot(y, s2.pdf(y), label=\"pdf\")\n",
+ "ax[\"b\"].legend()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "69e22499-9770-47e1-becb-326e29746947",
+ "metadata": {},
+ "source": [
+ "## Multivariate Sampling"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9c2a8fb8-239f-40a0-8fc4-d9776eb1c4a8",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# product sampler\n",
+ "s1 = pt.sampler.UniformSampler(domain=[0, 1])\n",
+ "s2 = pt.sampler.NormalSampler(mean=0, var=1)\n",
+ "s = pt.sampler.ProductSampler([s1, s2])\n",
+ "samples = s.sample(100_000)\n",
+ "x, y = np.linspace(0, 1, 50), np.linspace(-4, 4, 50)\n",
+ "xy = pt.misc.cartProd([x, y])\n",
+ "pdf = s.pdf(xy).reshape(x.size, y.size).T\n",
+ "\n",
+ "fig, ax = plt.subplot_mosaic([[\"x1\", \"pdf\"], [\"samples\", \"x2\"]], constrained_layout=True)\n",
+ "ax[\"pdf\"].contourf(pdf, 30, cmap=\"Greys\", extent=[0, 1, -4, 4])\n",
+ "ax[\"samples\"].scatter(samples[:, 0], samples[:, 1], alpha=0.05, s=1, color=\"grey\")\n",
+ "ax[\"samples\"].axhline(y=0, c=cmap.tab10(0))\n",
+ "ax[\"samples\"].axvline(x=0.5, c=cmap.tab10(1))\n",
+ "ax[\"x1\"].hist(samples[:, 0], bins=100, density=True, color=cmap.tab10(0))\n",
+ "ax[\"x2\"].hist(samples[:, 1], bins=100, density=True, color=cmap.tab10(1))\n",
+ "plt.show();"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "7903b7f1-36a4-4899-a3b2-2ffbc8abce8b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# parameter sampler\n",
+ "param1 = pt.parameter.Parameter(index=1, name=\"y1\", domain=[0, 1], distribution=\"uniform\")\n",
+ "param2 = pt.parameter.Parameter(index=2, name=\"y2\", domain=[-np.inf, np.inf], distribution=\"normal\",\n",
+ " mean=0, var=1)\n",
+ "params = [param1, param2]\n",
+ "s = pt.sampler.ParameterSampler(params)\n",
+ "samples = s.sample(100_000)\n",
+ "x, y = np.linspace(0, 1, 50), np.linspace(-4, 4, 50)\n",
+ "xy = pt.misc.cartProd([x, y])\n",
+ "pdf = s.pdf(xy).reshape(x.size, y.size).T\n",
+ "\n",
+ "fig, ax = plt.subplot_mosaic([[\"x1\", \"pdf\"], [\"samples\", \"x2\"]], constrained_layout=True)\n",
+ "ax[\"pdf\"].contourf(pdf, 30, cmap=\"Greys\", extent=[0, 1, -4, 4])\n",
+ "ax[\"samples\"].scatter(samples[:, 0], samples[:, 1], alpha=0.05, s=1, color=\"grey\")\n",
+ "ax[\"samples\"].axhline(y=0, c=cmap.tab10(0))\n",
+ "ax[\"samples\"].axvline(x=0.5, c=cmap.tab10(1))\n",
+ "ax[\"x1\"].hist(samples[:, 0], bins=100, density=True, color=cmap.tab10(0))\n",
+ "ax[\"x2\"].hist(samples[:, 1], bins=100, density=True, color=cmap.tab10(1))\n",
+ "plt.show();"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ba9900d0-b2e3-4bf0-9bbb-16c09df1eafc",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# WLS tensor sampler\n",
+ "param1 = pt.parameter.Parameter(index=1, name=\"y1\", domain=[0, 1], distribution=\"uniform\")\n",
+ "param2 = pt.parameter.Parameter(index=2, name=\"y2\", domain=[0, 1], distribution=\"beta\",\n",
+ " alpha=5, beta=5)\n",
+ "params = [param1, param2]\n",
+ "s = pt.sampler.WLSTensorSampler(params, [3, 3], tsa=True)\n",
+ "samples = s.sample(100_000)\n",
+ "x, y = np.linspace(0, 1, 50), np.linspace(0, 1, 50)\n",
+ "xy = pt.misc.cartProd([x, y])\n",
+ "pdf = s.pdf(xy).reshape(x.size, y.size).T\n",
+ "\n",
+ "fig, ax = plt.subplot_mosaic([[\"x1\", \"pdf\"], [\"samples\", \"x2\"]], constrained_layout=True)\n",
+ "ax[\"pdf\"].contourf(pdf, 30, cmap=\"Greys\", extent=[0, 1, 0, 1])\n",
+ "ax[\"samples\"].scatter(samples[:, 0], samples[:, 1], alpha=0.05, s=1, color=\"grey\")\n",
+ "ax[\"samples\"].axhline(y=0.5, c=cmap.tab10(0))\n",
+ "ax[\"samples\"].axvline(x=0.5, c=cmap.tab10(1))\n",
+ "ax[\"x1\"].hist(samples[:, 0], bins=100, density=True, color=cmap.tab10(0))\n",
+ "ax[\"x2\"].hist(samples[:, 1], bins=100, density=True, color=cmap.tab10(1))\n",
+ "plt.show();"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1dc2824d-e495-453b-9eba-f48bfdc03945",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.9.16"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}