diff --git a/examples/pyvbmc_example_3.ipynb b/examples/pyvbmc_example_3.ipynb
new file mode 100644
index 00000000..6a9caa58
--- /dev/null
+++ b/examples/pyvbmc_example_3.ipynb
@@ -0,0 +1,512 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Example 3: Using VBMC with PyMC models"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "In this notebook we show how to define a model in PyMC and the use VBMC to perform inference. For this example to work we need to have PyMC and ArviZ installed"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "import arviz as az\n",
+ "import numpy as np\n",
+ "import pymc as pm\n",
+ "from pyvbmc.vbmc import pymc_fit"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "az.style.use(\"arviz-white\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "For the first example we are going to reproduce the [Rosenbrock's banana](https://en.wikipedia.org/wiki/Rosenbrock_function) in $D = 2$ from example 1, but this time using PyMC.\n",
+ "\n",
+ "As the banana function is not available out-of-the-box from PyMC we have to write it ourselves."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "def log_banana(x):\n",
+ " return -pm.math.sum((x[0] ** 2.0 - x[1]) ** 2.0\n",
+ " + (x[0] - 1) ** 2.0 / 100)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This model do not have actual observed data, so we pass the `log_banana` function as a PyMC potential (an arbitrary term that we add to the log-likelihood term). We also define a prior as two independent $\\mathcal{N}(0, 3)$ distributions, for convenience we define it as vector of size 2 using the `shape=2` argument."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "scrolled": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Beginning variational optimization assuming EXACT observations of the log-joint.\n",
+ " Iteration f-count Mean[ELBO] Std[ELBO] sKL-iter[q] K[q] Convergence Action\n"
+ ]
+ },
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/u/32/martino5/unix/proyectos/00_BM/pyvbmc/pyvbmc/vbmc/pymc_fit.py:192: UserWarning: The effect of Potentials on other parameters is ignored during prior predictive sampling. This is likely to lead to invalid or biased predictive samples.\n",
+ " result = pm.sample_prior_predictive(\n"
+ ]
+ },
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ " 0 10 0.70 5.88 34645.41 2 inf start warm-up\n",
+ " 1 15 -2.80 0.33 5.31 2 inf \n",
+ " 2 20 -2.78 0.00 0.18 2 4.31 \n",
+ " 3 25 -2.81 0.00 0.07 2 1.86 \n",
+ " 4 30 -2.79 0.00 0.03 2 0.894 end warm-up\n",
+ " 5 35 -2.78 0.00 0.04 2 0.986 \n",
+ " 6 40 -2.75 0.00 0.03 2 0.843 \n",
+ " 7 45 -2.68 0.00 0.03 5 1.02 \n",
+ " 8 50 -2.59 0.00 0.02 6 0.875 rotoscale, undo\n",
+ " 9 55 -2.57 0.00 0.01 9 0.256 \n",
+ " 10 60 -2.54 0.00 0.01 12 0.251 \n",
+ " 11 65 -2.53 0.00 0.01 15 0.204 \n",
+ " 12 70 -2.45 0.00 0.05 17 1.51 \n",
+ " 13 75 -2.39 0.00 0.02 18 0.646 \n",
+ " 14 80 -2.39 0.00 0.00 16 0.0109 \n",
+ " 15 85 -2.38 0.00 0.01 15 0.157 stable\n",
+ " inf 85 -2.31 0.00 0.05 50 0.157 finalize\n",
+ "Inference terminated: variational solution stable for options.tolstablecountfcn evaluations.\n",
+ "Estimated ELBO: -2.306 +/-0.001.\n"
+ ]
+ }
+ ],
+ "source": [
+ "with pm.Model() as model:\n",
+ " x = pm.Normal(\"x\", 0, 3, shape=2)\n",
+ " y = pm.Potential(\"y\", log_banana(x))\n",
+ " vp, elbo, elbo_sd, idata = pymc_fit()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Time to inspect the posterior, for a quick check of a few summary statistics we can do:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {
+ "scrolled": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " mean | \n",
+ " sd | \n",
+ " hdi_3% | \n",
+ " hdi_97% | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | x[0] | \n",
+ " 0.043 | \n",
+ " 1.080 | \n",
+ " -1.958 | \n",
+ " 1.822 | \n",
+ "
\n",
+ " \n",
+ " | x[1] | \n",
+ " 1.117 | \n",
+ " 1.343 | \n",
+ " -1.095 | \n",
+ " 3.765 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "az.plot_pair(idata, marginals=True, kind=[\"kde\", \"scatter\"], figsize=(4, 4));"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Let us try with a diffent model, a Gaussian model with unkown mean and stantard deviation."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {
+ "scrolled": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "vbmc:TooCloseBounds: For each variable, hard and plausible bounds should not be too close. Moving plausible bounds.\n",
+ "Beginning variational optimization assuming EXACT observations of the log-joint.\n",
+ " Iteration f-count Mean[ELBO] Std[ELBO] sKL-iter[q] K[q] Convergence Action\n",
+ " 0 10 -156.34 7.14 376216.54 2 inf start warm-up\n",
+ " 1 15 -156.96 0.07 0.75 2 inf \n",
+ " 2 20 -156.87 0.00 0.12 2 3.07 \n",
+ " 3 25 -156.85 0.00 0.01 2 0.186 \n",
+ " 4 30 -156.85 0.00 0.00 2 0.0259 \n",
+ " 5 35 -156.85 0.00 0.01 2 0.315 \n",
+ " 6 40 -156.87 0.00 0.02 2 0.55 end warm-up\n",
+ " 7 45 -156.85 0.00 0.01 2 0.22 \n",
+ " 8 50 -156.84 0.00 0.00 2 0.0508 \n",
+ " 9 55 -156.81 0.00 0.00 5 0.141 \n",
+ " 10 60 -156.80 0.00 0.00 8 0.0599 rotoscale, undo\n",
+ " 11 65 -156.80 0.00 0.00 11 0.0197 \n",
+ " 12 70 -156.79 0.00 0.00 11 0.0453 stable\n",
+ " inf 70 -156.78 0.00 0.00 50 0.0453 finalize\n",
+ "Inference terminated: variational solution stable for options.tolstablecountfcn evaluations.\n",
+ "Estimated ELBO: -156.781 +/-0.001.\n"
+ ]
+ }
+ ],
+ "source": [
+ "data = np.random.normal(0, 1, size=100)\n",
+ "\n",
+ "with pm.Model() as model:\n",
+ " a = pm.Normal(\"a\", 0, 1)\n",
+ " b = pm.HalfNormal(\"b\", 1)\n",
+ " y = pm.Normal(\"y\", a, b, observed=data)\n",
+ " vp, elbo, elbo_sd, idata = pymc_fit()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " mean | \n",
+ " sd | \n",
+ " hdi_3% | \n",
+ " hdi_97% | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | a | \n",
+ " -0.098 | \n",
+ " 0.113 | \n",
+ " -0.317 | \n",
+ " 0.104 | \n",
+ "
\n",
+ " \n",
+ " | b | \n",
+ " 1.129 | \n",
+ " 0.075 | \n",
+ " 1.004 | \n",
+ " 1.254 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "