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
+}