{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Simulation-based calibration\n", "\n", "\n", "[Simulation-based calibration](https://arxiv.org/abs/1804.06788) (SBC) is a method for visually validating Bayesian inferences. SBC is useful for detection of either misspecified models, inaccurate computation or bugs in the \n", "implementation of a probabilistic program. So SBC is not a validation for the inference itself, but the technical aspect of the program." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Consider a generative model:\n", "\n", "\\begin{align}\n", "\\tilde{\\theta} & \\sim P(\\theta) \\\\\n", "\\tilde{y} & \\sim P(y \\mid \\tilde{\\theta}) \\\\\n", "\\{ \\theta_1, \\dots, \\theta_L \\} & \\sim P(\\theta \\mid \\tilde{y})\n", "\\end{align}\n", "\n", "The rank of a prior sample $\\tilde{\\theta}$ in comparison to an *exact* posterior sample $\\{ \\theta_1, \\dots, \\theta_L \\}$: \n", "\n", "\\begin{align}\n", "r(\\{ \\theta_1, \\dots, \\theta_L \\}, \\tilde{\\theta}) = \\sum_l^L \\mathbb{I}(\\theta_l < \\tilde{\\theta}) + 1\n", "\\end{align}\n", "\n", "is a discrete-uniform random variable in $[1, L + 1 ]$ (see the original paper for a proof). Thus we can use this as a testing procedure if our inferences work." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We follow Algorithm 2 from the original paper and implement SBC in Stan." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Loading required package: ggplot2\n", "Loading required package: StanHeaders\n", "rstan (Version 2.18.2, GitRev: 2e1f913d3ca3)\n", "For execution on a local, multicore CPU with excess RAM we recommend calling\n", "options(mc.cores = parallel::detectCores()).\n", "To avoid recompilation of unchanged Stan programs, we recommend calling\n", "rstan_options(auto_write = TRUE)\n" ] } ], "source": [ "library(rstan)\n", "suppressMessages(library(tidyverse))\n", "options(repr.plot.width = 6, repr.plot.height = 3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "SBC is essentially pretty simple. For $N$ iterations we run a while loop to create a posterior sample that has at least a target effective sample size $n_{eff}$ (which is set by us). We do so by resampling and thinning until the posterior sample for every iteration reaches the target $n_{eff}$ (within the while loop). Then we compute the rank for every parameter using the sum as defined above. That is it. If the inference worked, the ranks are uniformly distributed." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "sbc <- function(model, data)\n", "{\n", " ranks <- matrix(0, N, 2, dimnames = list(NULL, c(\"mu\", \"sigma\")))\n", " for (n in seq(N))\n", " { \n", " thin <- init_thin\n", " while (thin < max_thin) \n", " {\n", " fit <- suppressWarnings(\n", " sampling(model, data = data,\n", " chains = 1, iter = 2 * thin * L,\n", " thin = thin, control = list(adapt_delta = 0.99), refresh = 0)\n", " ) \n", " n_eff <- summary(fit)$summary[\"lp__\", \"n_eff\"]\n", " if (n_eff >= target_neff) break;\n", " thin <- 2 * thin\n", " } \n", " ranks[n,] <- apply(rstan::extract(fit)$idsim, 2, sum) + 1\n", " }\n", " ranks\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We start with a simple example where we set two parameters, generate data from them and then compare the posterior to these parameters. The comparison is done in the generated quantities block." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "data {\n", "\tint n;\n", "\tvector[n] x;\n", "}\n", "\n", "transformed data {\n", "\treal beta_sim = normal_rng(0, 1);\n", "\treal sigma_sim = lognormal_rng(0, 1);\n", "\t\n", "\tvector[n] y_sim;\n", "\tfor (i in 1:n)\n", "\t\ty_sim[i] = normal_rng(x[i] * beta_sim, sigma_sim);\n", "}\n", "\n", "parameters {\n", "\treal beta;\n", "\treal sigma;\n", "}\n", "\n", "model {\n", "\tbeta ~ normal(0, 1);\n", "\tsigma ~ lognormal(0, 1);\n", " \ty_sim ~ normal(x * beta, sigma);\n", "}\n", "\n", "generated quantities {\n", "\tint idsim[2] = { beta < beta_sim, sigma < sigma_sim };\n", "}\n" ] } ], "source": [ "model.file <- \"_models/sbc_1.stan\"\n", "cat(readLines(model.file), sep=\"\\n\")" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "model <- stan_model(model.file)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "?sampling" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We run the loop $5000$ times and sample $100$ times. We also set some other parameters that Stan or SBC needs, such as init_thin which specifies the period for saving samples or target_neff which is the effective sample size we want to have." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "N <- 5000\n", "L <- 100\n", "init_thin <- 1\n", "max_thin <- 64\n", "target_neff <- .8 * L" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We also define a histogram plotting method for the ranks." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "plot_fit <- function(fit)\n", "{\n", " as.data.frame(fit) %>% \n", " tidyr::gather(\"param\", \"value\") %>%\n", " ggplot(aes(value)) +\n", " geom_histogram(bins=30) +\n", " scale_y_continuous(\"Count\") + \n", " scale_x_continuous(\"\") + \n", " facet_grid(. ~ param) +\n", " ggthemes::theme_tufte()\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we run SBC and plot the ranks. We create some artifical data for the model first." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtAAAAFoCAIAAADxRFtOAAAABmJLR0QA/wD/AP+gvaeTAAAg\nAElEQVR4nO3deUCU1f7H8TPMMOw4KCoK5oYQojeXeyu1LJGflqWZmhtKimnmWqTV1Ui9maaW\n1yUtLeVqaWJd+f00UzM196XSUCpR1BYVV5Idhhnm98coYgGOOGeemYf36695zjwz8x3OHPhw\nnjPPo7FYLAIAAEAmN6ULAAAA6kfgAAAA0hE4AACAdAQOAAAgHYEDAABIR+AAAADSETgAAIB0\nBA4AACAdgQMAAEhH4AAAANIROAAAgHQEDgAAIB2BAwAASEfgAAAA0hE4AACAdAQOAAAgHYED\nAABIR+AAANiHKT81ODi426KflS4EzkindAEAAJXQuNcaOnRocMuaShcCZ6SxWCxK1wAAAFSO\nQypQTPsmDd87nfLqkGei2kY+/FifDWdyUj6b3v3xLhH3tnxy+LRsk0UIISzG4ODgBedzSx/V\nutE9L5/JUqxoAEIIITJ2fTK426PNQxu1aNPu+alLi0qut7dudM+YU9eEEKa89KnP9b2/RbP7\nOz0xZcW+2fff23fL78LGgS9E0R8HXxvWo3WL8EbNmj/Upf+Sr35T6I3CbggcUNKqYQsHvPPx\n9u9Snm9wbkznjq+fbLHuy6+O7P7PpW3Lx+w8r3R1AMpXnHOo6+BJ5g4DF3+aPGvis4f+M33w\nypO37mKeEPVk8oXgN99fMz9hxOkP4j7IuPlvgy0D/60eQ7+8ED55zkfJK+b3/XvO9OFdfzea\nHfgWYX+s4YCSgl9KaBXoKYR4bNy9r23Z9fErT7lrhHvttjF1vL747qroHKx0gQDKUXht21WT\n+bnnYx8N9BT/aN2sZsgpnxpld7h6bNJ/L2g2bZ/TwkcnRNu/JReEtX259F5bBn7dZ0bOiRnR\ntZanEKJ5pHbWisHH8oob6LWOfaOwJwIHlBQQ4W+9ofXUafX3GHQa66a3m0aUVPwwAIryqfd8\nj+afDvn73x9+4omHHnjwyae7P+53y1+Tsxu+96z5RAuf640+Qf2DPV4tvdeWgT967Kifjxz6\n76affv75x0N7t8h+R3AADqnApVhMuSUscwYU5qar+f5X369f8XbbumL7qpkdWrYaM+9A2R0s\nRovm1r8v7po7eP4S05VX+7bvOWrarrTMsHbdZ3ww3y5lQ1nMcMDJaYQQmcXXQ0buudX5ZqY+\nAIVd+T5x6daCSa+NavXIk0KIlA+feGrOlPdevDkPUb/bvYWJX54omBHmpRNCFFz54pdCk+2H\nSK+lv7HqwKVDp/bX12uFEPkXV9n9LcDxmOGAc9O4P1LDY/2E+UdPnT3+3baXBi9p4EFKBhSm\nN/yx+L0Z495bte/Ij7u3rFvx+W81Qh8vu0Odv7/bpWZRv0GTt+9PObxn8/h+U5t46jRaW2c5\nPAytLCXFqzbsPH3utx92rRvd610hxM69R/KZ4HRlBA44u4Wrp4dmftE7Onr0tA8f/ff6vh3+\n0dCDhWOAkvybxq+ZMiLt038/+3S3F16bfTkyZu3nY27Zw81z8Y51XbyPvzy0z/jpH/3jXxs6\nGzw9aultfH6foBGJr8YkvzU2uvNTb7y/bdB/vlkc98TeOVN+L+KLKi6ME3/BGZkLsy7medav\n5aF0IQCqwpT/07bdv0d16WpdulFiuto+rM2QvT+NrOdTyaMY+OrGDAeckdazBr90ANdlKSl8\necTw0Us3/no5L+fK78sm97vkef+woMrShmDgqx0zHAAA+8v4ZtnYqYu+Tb9kcfdt0jLqXwtn\nd2zoq3RRUBKBAwAgi8VUWKLzZNUVBIEDAAA4AGs4AACAdAQOAAAgHYEDAABIR+AAAADSETgA\nAIB0BA4AACAdgQMAAEhH4AAAANIROAAAgHQEDgAAIB2BAwAASEfgAAAA0hE4AACAdAQOAAAg\nHYEDAABIR+AAAADS6ZQuAAAAiQYOHFj5DqtXr3ZMJdUcMxwAAEA6AgcAAJCOwAEAAKQjcAAA\nAOkIHAAAQDoCBwAAkI7AAQAApCNwAAAA6TjxF1Slymf44dRAACAVMxwAAEA6AgcAAJCOwAEA\nAKQjcAAAAOkIHAAAQDoCBwAAkI7AAQAApCNwAAAA6QgcAABAOgIHAACQjsABAACk41oqAFSI\ni+MAzoYZDgAAIB0zHFXE/08AANiOwAEohtgKoPrgkAoAAJCOGQ4AgN0wb4eKMMMBAACkY4YD\nAOACKp87YeLE+UkMHCWmK+uXLdnybdrl7JL6jcJ6DBoR/bcg610HkxYm7TxyNkcb3rzN4DHP\nhfm5V94OAACqzBkOdUkMHJumTlyZHjB8fHyzALeUbZ8uTBhjXryia7BPelLCjDWnB48eExFg\n2rhk0RvxBauWxms1oqJ2AEDVMCsA5yErcJiNZz88ltnujTmPtw0UQoTe2zLjUL9PFxzt+nbb\nOZ+lNo15t090EyFE6Gy3Z2JnrTg3Ii5YX357iK+kCp2TM4RQAADsTlbgMBWkN2rc+MkIw40G\nTSt/j/3ZeYXXdmQYzaM617e2ehg6tPbVp2y/UNj9VLntIjZUUoUA4BL4PwTqICtweNR4dN68\nR0s3Cy8dXn4+t2FceHH+WiFEhPfNxRkR3rpNqVnFnVPLbS/dPH/+fFbW9U2tVhsWFiapcgAA\nYHeO+JbK6YMb5ryzvLhhl0mPhZh/zRdCBOpufh030F1rzisqKSq/vXRz8eLFmzdvtt4OCAjY\nunWrAypH5fjHC6okad0DyylQzckNHMbsU4lz3/0yJfOhp0eOHtTF202To/cSQmSaSny0Wus+\nV4vNWoPerYJ2qeUBAADHkBg48jN2vzhurqnZ/8z+MC480NPa6O7TQohdJwpMDTyuB4tfCs3+\nkTUqai99tvj4+BdeeMF6282N85UBFWLmCYATkvaX22KaOWG+R6eRS98aVZo2hBAehqggvXbL\n3kvWTVNB2oEcY6vooIraSx9Ys2bN4Bvq1asnq2wAACCBrBmO/Isfp+QYn23p9/2hgzdfzCus\n7d8CJvSOfDVx2rZ6E5sbTBsWv+0V0mlIiJ9GI8ptl1QeAABwJFmBIzs9XQixYs7bZRv9G0z6\nZNGDYQOmvyLmJy2duSRXFx7Zce7Lw6xn96qoHQAAuDpZgSPoobfWP1Thve0GjG834A7aAQCA\nS2P1JQAAkI7AAQAApCNwAAAA6QgcAABAOgIHAACQjsABAACkI3AAAADpHHG1WGfG9RsBAOVy\nwssSOWFJtmOGAwAASFfdZzggXDwyAwBcAjMcAABAOgIHAACQjsABAACkI3AAAADpCBwAAEA6\nvqUC4Ca+sgRAEmY4AACAdAQOAAAgHYEDAABIxxoOAADsjOVQf8UMBwAAkI4ZDgDAn/EPOuyO\nwAFl8OsMAKoVAgdwtwhPAHBbrOEAAADSETgAAIB0HFKBk6rkOAVHKADA5RA41KPylQT8kQYA\nKIhDKgAAQDpmOADcAQ51weq2X84C/oTAAbgkft3j7vEpgiNxSAUAAEjHDAcAANWdA752wAwH\nAACQjsABAACkI3AAAADpWMNRGZZwQ1l8AgGoBoEDroc/wwDgcjikAgAApCNwAAAA6RxxSGXJ\n0H4+sxIH1fG2bl7cP3n4zGNld4hLXNuzlqcQ4mDSwqSdR87maMObtxk85rkwP3cHlAcALo2D\njHAJkgOHxZiybdnGqwV9y7RdS7nmVav7+OGRpS0N/dyFEOlJCTPWnB48ekxEgGnjkkVvxBes\nWhqv1cgtEAAAOIDEwHFhz5yX5u3NM5b8qf3ST9mG5u3bt4+8pdVinPNZatOYd/tENxFChM52\neyZ21opzI+JCfOVVCAAAHENi4AhsFTtrbr8S44Vx8dPLth/NKgroYTAXZF/OKalbx2Cdwii8\ntiPDaB7Vub51Hw9Dh9a++pTtF0RsqLwKpZJ0mljmTgFlMQadE/3i/CQGDp1v3Xt8hdn453Wp\nh/OKS/Ys6LvweLHFovOu02Poi0O6tijOTxVCRHjfXLQR4a3blJpVurlr164zZ85Yb3t6evbr\n109e5Q7A2ABQPTngmh1wTo4+D4fZeC7TLBobHnhz2eu1PYoObUp8Z/FkryYrH9PmCyECdTfT\nSaC71pxXVLr51Vdfbd682Xo7ICDA1QMHqhXyJQA4OnBo9cHJyck3tvwe7jvxxKbvNr5/tFu8\nlxAi01Tio9Va77tabNYa9KUP9PLy8vf3v/4wPz9H1gwAAO6S8mcabV3X6+vMK+4+LYTYdaLA\n1MDjeuD4pdDsH1mjdLfJkydPnjxZoRoBAPgzJi/viKNP/JV1amnMoGHnjOYbDZbd5/NrRIR5\nGKKC9Notey9ZW00FaQdyjK2igxxcHgAAkMHRgcO/UUwzbc6kKR8cOpaW/lNK0oJXvsn3jx8Z\nrtHoJ/SOPJk4bdvhtIzTPy5PmO4V0mlICIdOAABQA0cfUtFoff45L2H54pWL3k7IE75NQu+b\nunBamJdOCBE2YPorYn7S0plLcnXhkR3nvjyMs34BMjAPDMDxpAcOrT5k/fr1ZVs8Alq+MHnO\nC+Xt3G7A+HYDZFcEAAAcTflFo3AM/qm14ucAAIrgarEAAEA6ZjgAAHC0SmZb1Xq6VWY4AACA\ndLbOcLRr1673Z1sn/OXarRf2jXvm9T92b//Y3oUBAOAIalra5czv5TaBIy0tzXrjwIEDTX7+\nOS3P/5a7LaZD//vNvt2/SSoOAFCtOPPfS4dR6w/hNoHj3nvvLb29usv95R5WqtFknF1LgtNR\n66f/jvBDAIC7cZvA8cEHH1hvjBw58pE3/z2gttefdtC6+3d4po+U0gAAgFrcJnA8//zz1htr\n1qzpGffc8/X/vIYDAADgtmxdNLpjxw4hROa505dzi/96b0izMB83zkMOAADKZ2vgKLyyrc/D\nfTcezyz33h9yjff5uNuvKgCQi0U5gIPZGjiWPjVoy681xiVMbB4c8Nd7I71JGwAAoEK2Bo7p\n316O2/r7/EfqSa0GAACokq1nGvXWavrdV1NqKQAAQK1sDRwJ99f5aEeG1FIAAIBa2Ro4Bn+x\n/o9/Pj5z5bZ8s0VqQQAAQH1sXcPRuduL5hrFk56NnjzUo25wkKf2li/BnjlzRkJtAABAJWwN\nHIGBgUIE9uzZUmo1UBm+eQgAsLI1cCQnJ0utAwAAqJitgSMrK6uiuzRaL39fvZ3qAYDqiOlA\nqJ6tgcNgMFR0V0Do4syTL9ipHgAAoEK2Bo6pU6eW3SwpyjlzMvV/k7caHh07P76L/esCAAAq\nYmvgmDJlyl8bc3/5+v7Ibpvzxj9t15oAAIDK2HoejnL5Nope+2brT8a8bq9qAACAKt1V4BBC\n+DT0Kcz8wi6lAAAAtbqrwFFSfGnu6z/ovMLsVQ0AAFAlW9dwtGvX7i9t5nMnjv6eWfSPhPfs\nWxMAAFAZWwNHebQN74vu1XnQ7EkP2K0cAACgRrYGjv3790utAwAAqNidzXAUXU1d/+WB9PST\nl4t9wsLCHnz8qVZ1vSRVBgAAVOMOAse66cNHT0+8UGQubdG61x7yxgcfvd5LQmEAAEA9bP2W\nypnPBvZO+Ej7QL+PN+85+euFy2dP79u6un87/bKE3oP++4vMCgEAgMuzdYZjzosbfOsPOLbt\nkwCdxtoSGNz4wUcftzRq8H/j5ojei6RVCAAAXJ6tMxxJl/PDRr5SmjasNDrDq2Pvzb+8RkJh\nAABAPWwNHD5uboUXC//aXnSpyE3ra9eSAACA2tgaOMY1rXEyMW731VsyR9G1/cOWHK/RdJyE\nwgAAgHrYuoZj2OdTprV8sfM94YPGDn/w3qb+mtxTaYc+Wrjit0L9vz+Lk1oiAABwdbYGjoCI\nsce/qTX2pQmJsxISbzTWbfPUinnvDYoIkFQcAABQhzs4D0fwQwPXfTvgytnTJ0+evGbxb9as\nWZMGte/2arMAAKAauNNrqWgCQ5oGhjSVUgsAAFApmwJH9m8p+66EPNamlnWzMPOL8TMOR0VF\ndevawU+rqfyxQoglQ/v5zEocVMe7tOVg0sKknUfO5mjDm7cZPOa5MD/3ytsBAIBLu80hEWNW\n6qjoyIBGrSdv+L200Vz069J3p/R/4uGGrXpuPJ5V2eMtxpSv3994taBsW3pSwow1B9r3Gj7l\nxVi/MzvfiF9otlTWDgAAXF1lMxzFeUcfbvrAt5lF7Z6KG9OjQWm7d92477YHf/H5R2+9v6Fn\n61ZfZxx/xODx14df2DPnpXl784wlt7RajHM+S20a826f6CZCiNDZbs/EzlpxbkRcsL789hBO\n8gEAgMurbIZjz/i+h64WPpf43d7kjwa0rlXarnHzatup55RFX5zYNMFU+MvwsXvKfXhgq9hZ\ncxcsmPt62cbCazsyjOYunetbNz0MHVr76lO2X6io/a7eHAAAcA6VzXC8k/ybb73nlz7bpqId\nGnWdPbXJh7O+/LcQnct5at+69/gKs/GWTFOcnyqEiPC+uTgjwlu3KTWruHP57aWba9as+eGH\nH6y3fXx8EhISKn9jAADAeVQWOHZnFdXtNaDyx0c/UvdfH++1/fVKivKFEIG6mykk0F1rziuq\nqL10MzU19euvv7beDgjgzB8AALiSygJHXb2bubikkh2EEMY/jG66Gra/npveSwiRaSrx0Wqt\nLVeLzVqDvqL20ge2aNHCZDJZb/v4+Nj+igAAQHGVBY7etb3f258oRKdK9lm275JnzWdtfz13\nnxZC7DpRYGrgcT1Y/FJo9o+sUVF76QP79+/fv39/218IAAA4j8oWjcZNbZd38ePBK45VtMPx\n1XGrLuW1HD/S9tfzMEQF6bVb9l6ybpoK0g7kGFtFB1XUbvszAwAAp1VZ4Gj27LqYMMOquLY9\nR03//lx+2buKc07Nm9C/bewKr9qPrnsp0vbX02j0E3pHnkyctu1wWsbpH5cnTPcK6TQkxK+i\n9iq+LQAA4EwqO6SicfNOPLxX27XLyvcT1i+Z0TA0tHHjxrW9ik+np6efOH3NaPZv0vnzPeuD\n3O/sgiphA6a/IuYnLZ25JFcXHtlx7svDrGcrragdAAC4Oo3FcrvTeVpMO5IWL/nPf1PSTp76\n7UJxicW/zj3h4eH/0zsufnTfWjrXvnzbwIEDlS4BUInVq1crXcJNDG3Ajuwyum24lopG16n/\nuE79xwkhSopzrxa41/Yv57yiAAAAFbmzq8W6ufvW5npqAADgDrn2AREAAOASCBwAAEA6AgcA\nAJCOwAEAAKQjcAAAAOkIHAAAQDoCBwAAkI7AAQAApCNwAAAA6QgcAABAOgIHAACQjsABAACk\nI3AAAADpCBwAAEA6AgcAAJCOwAEAAKQjcAAAAOkIHAAAQDoCBwAAkI7AAQAApCNwAAAA6Qgc\nAABAOgIHAACQjsABAACkI3AAAADpCBwAAEA6AgcAAJCOwAEAAKQjcAAAAOkIHAAAQDoCBwAA\nkI7AAQAApCNwAAAA6QgcAABAOgIHAACQjsABAACkI3AAAADpCBwAAEA6neNf8uL+ycNnHivb\nEpe4tmctTyHEwaSFSTuPnM3RhjdvM3jMc2F+7o4vDwAA2J0CgeNayjWvWt3HD48sbWno5y6E\nSE9KmLHm9ODRYyICTBuXLHojvmDV0nitxvEFAgAAO1MgcFz6KdvQvH379pG3tFqMcz5LbRrz\nbp/oJkKI0Nluz8TOWnFuRFyIr+MrBAAA9qXAGo6jWUUBrQ3mguwLl65ZbjQWXtuRYTR36Vzf\nuulh6NDaV5+y/YLjywMAAHanwAzH4bzikj0L+i48Xmyx6Lzr9Bj64pCuLYrzU4UQEd43F21E\neOs2pWaVbr711lvbtm2z3jYYDOvWrXNw2QAAoMocHTjMxnOZZtHY8MCby16v7VF0aFPiO4sn\nezVZ+Zg2XwgRqLs54xLorjXnFZVuFhQUZGdnW29rtVrbX3HgwIF2qh0AAFSRowOHVh+cnJx8\nY8vv4b4TT2z6buP7R7vFewkhMk0lPjfCxNVis9agL31gly5dmjVrZr3t6enpyJoBAMBdUuCQ\nyp+0ruv1deYVd58WQuw6UWBq4HE9cPxSaPaPrFG6W8eOHTt27KhQjQAA4K44etFo1qmlMYOG\nnTOabzRYdp/PrxER5mGICtJrt+y9ZG01FaQdyDG2ig5ycHkAAEAGRwcO/0YxzbQ5k6Z8cOhY\nWvpPKUkLXvkm3z9+ZLhGo5/QO/Jk4rRth9MyTv+4PGG6V0inISF+Di4PAADI4OhDKhqtzz/n\nJSxfvHLR2wl5wrdJ6H1TF04L89IJIcIGTH9FzE9aOnNJri48suPcl4dx1i8AANRBY7FYbr+X\nK+NbKoBjrF692sGvyOgGHMMuo5uLtwEAAOkIHAAAQDoCBwAAkI7AAQAApCNwAAAA6QgcAABA\nOgIHAACQjsABAACkI3AAAADpCBwAAEA6AgcAAJCOwAEAAKQjcAAAAOkIHAAAQDoCBwAAkI7A\nAQAApCNwAAAA6QgcAABAOgIHAACQjsABAACkI3AAAADpCBwAAEA6AgcAAJCOwAEAAKQjcAAA\nAOkIHAAAQDoCBwAAkI7AAQAApCNwAAAA6QgcAABAOgIHAACQjsABAACkI3AAAADpCBwAAEA6\nAgcAAJCOwAEAAKQjcAAAAOkIHAAAQDoCBwAAkI7AAQAApNMpXcAtDiYtTNp55GyONrx5m8Fj\nngvzc1e6IgAAYAdONMORnpQwY82B9r2GT3kx1u/MzjfiF5otStcEAADswWkCh8U457PUpjFv\n9oluF9n24fGzx+Zf/GbFuVylywIAAHbgLIdUCq/tyDCaR3Wub930MHRo7atP2X5BxIZaWzIz\nMwsKCqy33dzc6tWrp0yhAADgzjlL4CjOTxVCRHjfXLQR4a3blJpVujl37tzNmzdbbwcEBGzd\nutXBFQIAgCpzlkMqJUX5QohA3c16At215rwi5SoCAAB24ywzHG56LyFEpqnER6u1tlwtNmsN\n+tIdRo0aFRMTY72tvbGPLVavXm2/MgE4EUY34EKcJXC4+7QQYteJAlMDj+th4pdCs39kjdId\n6tevX79+fYWqAwAAd8VZDql4GKKC9Notey9ZN00FaQdyjK2ig5StCgAA2IWzBA6NRj+hd+TJ\nxGnbDqdlnP5xecJ0r5BOQ0L8lK4LAADYgcZicaKza+3/dH7SziPnc3XhkQ+MfnlYkN5Z8hAA\nALgbzhU4AACAKjGFAAAApCNwAAAA6QgcAABAOgIHAACQjsABAACkq16BIzY2NioqatmyZUoX\nYjdr166Niorq1auX0oXYTW5ublRUVFRU1L59+5SuxW6mTZsWFRU1efJkpQuxm2+//dbaTX/8\n8YfStQghxL59+6z15ObmKl2L3fTq1SsqKmrt2rVKF2I3y5Yti4qKio2NVboQu8nIyLB+8I4e\nPap0LXYzYcKEqKiomTNn2v2ZneXU5o6Rl5eXnZ1dVKSea8IVFRVlZ2frdOrpR4vFkp2dLYQw\nmUxK12I3BQUF2dnZBQUFShdiNyaTydpNTvK9emerxy5yc3NV+fsqLy9P6ULspvT3ldlsVroW\nu8nPz8/Ozi4sLLT7M1evGQ4AAKAIAgcAAJBOPVPxtmjfvn2zZs1CQ0OVLsRuGjZsGB0d7evr\nq3QhdqPT6aKjo4UQgYGBStdiNy1bthRCREREKF2I3dSqVcvaTXq9XulahBAiMDDQWo+aDi92\n7NgxNze3YcOGShdiN6GhodHR0bVr11a6ELvx9PS0fvAMBoPStdhNmzZt/P39IyMj7f7MnNoc\nAABIxyEVAAAgHYEDAABIp57jnbd1MGlh0s4jZ3O04c3bDB7zXJifu9IVVUWJ6cr6ZUu2fJt2\nObukfqOwHoNGRP8tSAhxcf/k4TOPld0zLnFtz1qeCpV5Byqp3EW7LOfcuzEv7PxTo97nvs8/\nfdNFu2nJ0H4+sxIH1fEubamoa5TqMhf9qJSlvqEtVDe6Gdp32V/VJXCkJyXMWHN68OgxEQGm\njUsWvRFfsGppvFajdFl3btPUiSvTA4aPj28W4Jay7dOFCWPMi1d0Dfa5lnLNq1b38cNvLvNp\n6AoDWAhRUeWu22XeNbu/9lq7si0Hli84GdlFVPxmnZfFmLJt2carBX3LtFXUNUp1met+VMpS\n39AWqhvdDO277S9LdVBSNKJ3z5c+O2XdKvxjT/fu3Zf9nqNsUVVgKvr9qR493v7u8o2GkoWx\nzzz7yj6LxbJr7KARs1MVrK3Kyq9cLV1msVgyf1zVJybhanGJxdW6KWP37P69n+revXv37t0/\nvph3vbWirlGqy1TxUVHl0LaofXQztO9UtVjDUXhtR4bR3KVzfeumh6FDa199yvYLylZVBaaC\n9EaNGz8ZUfr9K00rf4/i7DwhxNGsooDWBnNB9oVL11zre0flVq6aLrOY/pj5r//2e3NiTZ1G\nuFo3BbaKnTV3wYK5r5dtrKhrlOoydXxUVDm0hapHN0O7Cq9bLQ6pFOenCiEivG9OcEV46zal\nZilXURV51Hh03rxHSzcLLx1efj63YVy4EOJwXnHJngV9Fx4vtlh03nV6DH1xSNcWihV6J8qt\nXDVddjp5Rka9fn0a+1k3XaubdL517/EVZuMt/5ZU1DXFnZXpMnV8VFQ5tIWqRzdDuyqvW8V6\nXUpJUb4QIlB384cb6K4157n2FQpOH9ww553lxQ27THosxGw8l2kWjQ0PvLns9doeRYc2Jb6z\neLJXk5X9mtVQuszbqKjyx7Rq6DKz8dzMpPQBi6aVbrpoN5VV0WhSapSpb3SrY2gLVY9uhnbV\n+qtaBA43vZcQItNU4qPVWluuFpu1Bqc4Q2IVGLNPJc5998uUzIeeHjl6UBdvN43QBicnJ9+4\n3+/hvhNPbPpu4/tH+819WMlCbaDVl195t3g1dNn5r+bn+HXpVvf6CvCK3qzzd1NZFY0mpUaZ\nmka3moa2UPXoZmhXrb+qxRoOd58WQogTBTevPvpLodk/0pWyZ6n8jN1jhk04aGwx+8P/TIzt\n6u1W/kLh1nW9inOvOLg2u7BWroous3y85nTooB6V7OGK3VRR1yjVZar4qJD2dQAAAAJZSURB\nVAhRDYa2UM/oZmhXsb+qReDwMEQF6bVb9l6ybpoK0g7kGFtFBylbVVVYTDMnzPfoNHLpW6PC\nA29+wzvr1NKYQcPOGUuvj2zZfT6/RkSYIjXekYoqV0GXFVxJPphTPLTDzZpdt5vKqqhrlOoy\nFXxUhFDh0BbqHd0M7Sr3V7U4pKLR6Cf0jnw1cdq2ehObG0wbFr/tFdJpSIif0nXdsfyLH6fk\nGJ9t6ff9oYOljTqvsDaRMc20X0+a8sHogdE1tYXff/3JN/n+s0aGK1iqjfwblV+5RqNz9S47\nv2mXu2/bUE9taUtFb1bBIqugotGk0QhFukwdo1t9Q1uod3QztKvcX9Xo4m37P52ftPPI+Vxd\neOQDo18eFqR3vdmdC3smj5h97E+N/g0mfbLowaI/ji1fvPLAT7/mCd8moffFjBx+Xz3vcp/E\n2VRSuUt32fK4fnuCJyx/8x9lG12xm8zGs0/3GdX3ozVlT0dYUdco1WUu/VERKh3aQqWjm6Fd\n5f6qRoEDAAAoxWVCJQAAcF0EDgAAIB2BAwAASEfgAAAA0hE4AACAdAQOAAAgHYEDAABIR+AA\nAADSETgAAIB0BA4AACAdgQMAAEhH4AAAANIROAAAgHQEDgAAIB2BAwAASEfgAAAA0hE4AACA\ndAQOAAAgHYEDAABIR+AAAADSETgAAIB0BA4AACAdgQMAAEhH4AAAANIROAAAgHQEDgAAIB2B\nAwAASEfgAAAA0hE4AACAdAQOAAAgHYEDAABIR+AAAADSETgAAIB0BA4AACAdgQMAAEhH4AAA\nANL9P9AtLaC42WydAAAAAElFTkSuQmCC", "text/plain": [ "plot without title" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "n <- 10\n", "x <- rnorm(n)\n", "\n", "fit <- sbc(model, list(x=x, n=n))\n", "plot_fit(fit)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Given the low number of trials, the histogram looks fairly uniform (as expected since the model was specified correctly)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we test some pathological cases to see if the ranks change with mis-specified models." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "data {\n", "\tint n;\n", "\tvector[n] x;\n", "}\n", "\n", "transformed data {\n", "\treal beta_sim = normal_rng(0, 1);\n", "\treal sigma_sim = lognormal_rng(0, 1);\n", "\t\n", "\tvector[n] y_sim;\n", "\tfor (i in 1:n)\n", "\t\ty_sim[i] = normal_rng(x[i] * beta_sim, sigma_sim);\n", "}\n", "\n", "parameters {\n", "\treal beta;\n", "\treal sigma;\n", "}\n", "\n", "model {\n", "\tmu ~ normal(0, 10);\n", "\tsigma ~ lognormal(0, 5);\n", " \ty_sim ~ student_t(10, x * beta, sigma);\n", "}\n", "\n", "generated quantities {\n", "\tint idsim[2] = { beta < beta_sim, sigma < sigma_sim };\n", "}\n" ] } ], "source": [ "model.file <- \"_models/sbc_2.stan\"\n", "cat(readLines(model.file), sep=\"\\n\")" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtAAAAFoCAIAAADxRFtOAAAABmJLR0QA/wD/AP+gvaeTAAAd\n30lEQVR4nO3deVxU9f7H8e8ww7DKopgoGi4IInpzuS1mmo5crMwlNU1RUlyvaxFWD43UNE0t\nykxTS0nLBevG/Wml1lVzX1q8GN1CUSsXXFHZlxnm98cUYjE44Hw5M4fX869zvnMYP+P3fOE9\n37NpzGazAAAAkMlF6QIAAID6ETgAAIB0BA4AACAdgQMAAEhH4AAAANIROAAAgHQEDgAAIB2B\nAwAASEfgAAAA0hE4AACAdAQOAAAgHYEDAABIR+AAAADSETgAAIB0BA4AACAdgQMAAEhH4AAA\nANIROAAA9mHMTwsKCnps6U9KFwJHpFO6AACASmhc640cOTKobV2lC4Ej0pjNZqVrAAAAKsch\nFSjmwebB75xKfWHEk4aOEV0eGbjldE7qx3N7PxoV3qrt42NmZxvNQghhLg4KCnr7fG7ZT7Vv\nevdzp28oVjQAIYQQmXs+Gv5Yt9YhTdt06DRu1sqi0t/b2ze9e9LJ60IIY17GrNGD7mvT8r7u\nvWauObDwvlaDtp8RNg58IYquHX5xVJ/2bcKatmz9UNRTK778TaEPCrshcEBJ60YtGfL6hzu/\nTR3X5NykHl1fOtHm0y++PLr3g0s7Vk/afV7p6gBUrCTnSM/h002dhy7bkLJg2tNHPpg7fO2J\nWzcxxRseT7kQNOfdjYsTxp5aHrs88+bXBlsG/qt9Rn5xIWzGovdT1iwe9PecuWN6nik21eBH\nhP1xDgeUFPRsQrsAdyHEI1Navbh9z4fP93XVCNf6HaPv8vjs26uiR5DSBQKoQOH1HVeNptHj\nYroFuIt727es2/ikl2/5Da7+MP1fFzRbdy5q46UTouPfUgpCOz5X9qotA7/Bk+MXRY/tWc9d\nCNE6QrtgzfAf8kqa6LU1+0FhTwQOKMk/3MeyoHXXafV3++k0llVPF40otf5jABTl1XBcn9Yb\nRvz971169Xro/gcef6L3o3Vu+Wtydst37nV7tfH6vdEr8KkgtxfKXrVl4E+cPOGno0f+tfV/\nP/3045H922V/ItQADqnAqZiNuaWc5gwozEVX990vv9u85rWODcTOdfM7t2036a1D5TcwF5s1\nt/59cdVU4f1LjVdeGPRgvwmz96RnhXbqPW/5YruUDWUxwwEHpxFCZJX8HjJyz63PNzH1ASjs\nyndJK78qmP7ihHYPPy6ESH2vV99FM9955uY8RKPHWhUmfXG8YF6oh04IUXDls18KjbYfIr2e\n8fK6Q5eOnDzYSK8VQuRfXGf3j4CaxwwHHJvG9WFft83xi4+dPPvztzueHb6iiRspGVCY3u/a\nsnfmTXln3YGjP+7d/umaT37zDXm0/AZ3/f2NqLpFg4fN2Hkw9ft926YOntXcXafR2jrL4ebX\nzlxasm7L7lPnfvvvnk8n9n9DCLF7/9F8JjidGYEDjm7J+rkhWZ8NiIycOPu9bm9uHtT53mA3\nThwDlOTTIm7jzLHpG958+onH/vniwssR0Zs+mXTLFi7uy3Z9GuX583MjB06d+/69r2zp4efu\nVk9v4/t7BY5NeiE65dXJkT36vvzujmEffL0sttf+RTPPFHGhihPjxl9wRKbCGxfz3BvVc1O6\nEADVYcz/3469ZwxRPS2nbpQarz4Y2mHE/v+Nb+hVyU8x8NWNGQ44Iq27L790AOdlLi18buyY\niSs///VyXs6VM6tmDL7kft+owMrShmDgqx0zHAAA+8v8etXkWUu/ybhkdvVu3tbwypKFXYO9\nlS4KSiJwAABkMRsLS3XunHUFQeAAAAA1gHM4AACAdAQOAAAgHYEDAABIR+AAAADSETgAAIB0\nBA4AACAdgQMAAEhH4AAAANIROAAAgHQEDgAAIB2BAwAASEfgAAAA0hE4AACAdAQOAAAgHYED\nAABIR+AAAADS6ZQuAAAAyDV06NDKN1i/fr3sGpjhAAAA0hE4AACAdAQOAAAgHYEDAABIR+AA\nAADSETgAAIB0BA4AACAdgQMAAEhH4AAAANIROAAAgHQEDgAAIB2BAwAASEfgAAAA0vG0WAAA\nnIMjPPS12pjhAAAA0hE4AACAdAQOAAAgHYEDAABIR+AAAADSETgAAIB0BA4AACAdgQMAAEjH\njb8AAKiAU99lywExwwEAAKQjcAAAAOkIHAAAQDoCBwAAkI7AAQAApCNwAAAA6bgsFgCgZlzd\n6iCY4QAAANIxwwEAgJ0xrfJXtT1wVL5P1MIdAgAAGTikAgAApCNwAAAA6QgcAABAOgIHAACQ\njsABAACkI3AAAADpCBwAAEA6AgcAAJCOwAEAAKSr7XcaBQA4BW4M7ewIHAAA1HY1kOc4pAIA\nAKQjcAAAAOkIHAAAQDrO4QAA2E3lpwIIzu6sxZjhAAAA0jHDAQjB1zIAkIwZDgAAIB2BAwAA\nSMchFagKR0YAwDExwwEAAKQjcAAAAOkIHAAAQDoCBwAAkI7AAQAApJN4lUqp8crmVSu2f5N+\nObu0UdPQPsPGRv4t0PLS4eQlybuPns3RhrXuMHzS6NA6rpW3Q324nAQAahWJMxxbZ01bu+tq\nn1Fxr73y/MPBBUsSJm0/lyeEyEhOmLfx0IP9x8x8JqbO6d0vxy0xmUUl7QAAwNnJmuEwFZ99\n74esTi8verRjgBAipFXbzCODN7x9rOdrHRd9nNYi+o2Bkc2FECELXZ6MWbDm3NjYIH3F7Y29\nJVUIqBWzR7hz7EWwO1mBw1iQ0bRZs8fD/f5o0LTzcTuYnVd4fVdmsWlCj0aWVje/zu299ak7\nLxT2Pllhu4gJkVQhAAB34raxDOXJChxuvt3eeqtb2Wrhpe9Xn88Njg0ryd8khAj3vHlyRrin\nbmvajZIeaRW2l61++eWX6enplmUPD4/Ro0dLqhwAANhdTdza/NThLYteX10SHDX9kcamX/OF\nEAG6m+eOBLhqTXlFpUUVt5et7tmzZ9u2bZZlf39/AgcAAE5EbuAozj6ZlPjGF6lZDz0xfuKw\nKE8XTY7eQwiRZSz10mot21wtMWn99C5W2sveqm7dukFBQZZlX19fqWUDAAD7khg48jP3PjMl\n0djyHwvfiw0LcLc0unq1EWLP8QJjE7ffg8UvhSafCF9r7WXvFhcXFxcXJ69aAAAgj7TLYs3G\n+fGL3bqPX/nqhLK0IYRw8zME6rXb91+yrBoL0g/lFLeLDLTWLqs8AABQg2TNcORf/DA1p/jp\ntnW+O3L45j/mEdrxb/7xAyJeSJq9o+G01n7GLcte82jcfUTjOhqNqLBdUnmAI+DKQ6DWqmT4\nq3Xgywoc2RkZQog1i14r3+jTZPpHSx8IHTL3ebE4eeX8Fbm6sIiuic+N0mqEEMJaOwAA8nB1\na82QFTgCH3p180NWX+00ZGqnIVVoB1AzmHQBIElNXBYLAFBE5QmS+IiaROBwLHy/BABUmyMf\nHuLx9AAAQDoCBwAAkI5DKurBwVqlcCAMAG6LGQ4AACAdMxyoPr7ZA3AQjnyyJCwIHAAAOBC1\nhicCRzXx5R4AANtxDgcAAJCOwAEAAKQjcAAAAOkIHAAAQDoCBwAAkI6rVMAVNwAA6ZjhAAAA\n0hE4AACAdAQOAAAgHYEDAABIR+AAAADSETgAAIB0BA4AACAd9+GQpfKbW3BnCwBArcIMBwAA\nkI4ZDsAmt70fKwCgEsxwAAAA6WwNHJ06dXr9bO5f2y8cmNLFMNyuJQEAALW5zSGV9PR0y8Kh\nQ4ea//RTep7PLS+bjUf+/fWBvb9JKg4AAKjDbQJHq1atypbXR91X4ZUVvs2n2LUkAIBqcQVf\nrXWbwLF8+XLLwvjx4x+e8+aQ+h5/2kDr6tP5yYFSSgMqwsmbAOCMbhM4xo0bZ1nYuHFjv9jR\n4xp5yy8JAACoja2Xxe7atUsIkXXu1OXckr++2rhlqJeLxp51qR1f0wEAtYqtgaPwyo6BXQZ9\n/nNWha/+N7f4Hi9X+1UFAKiN+DKmYrYGjpV9h23/1XdKwrTWQf5/fTXCk7QB1AqV/D3gdD8A\nlbA1cMz95nLsV2cWP9xQajUAAECVbA0cnlrN4HvqSi0FAFCT7uT4Bcc+UFW2Bo6E++56f1em\n4YmmMosB1IlfzQBga+AY/tnmT+59dH7OO1OjDZ5aLkipXRT5eynpH+VvPwAowtbA0eOxZ0y+\nJdOfjpwx0q1BUKD7rZnj9OnTEmoDAAAqYWvgCAgIECKgX7+2UqsBAACqZGvgSElJkVoHgCrh\n2BAA52Jr4Lhx44a1lzRaDx9vvZ3qAQAAKmRr4PDz87P2kn/IsqwT/7RTPY6FL5FQFnsgANWw\nNXDMmjWr/GppUc7pE2n/TvnKr9vkxXFR9q8LACCEsCF3co9XOAVbA8fMmTP/2pj7y3/ui3hs\nW97UJ+xaEwAAUBmXO/lh76aRm+a0/2jSS/aqBgAAqNIdBQ4hhFewV2HWZ3YpBQAAqNUdBY7S\nkkuJL/1X5xFqr2oAAIAq2XoOR6dOnf7SZjp3/NiZrKJ7E96xb02oBJctAACcka2BoyLa4Hsi\n+/cYtnD6/XYrBwBUqvJvC1xpAtWzNXAcPHhQah0AAEDFqjbDUXQ1bfMXhzIyTlwu8QoNDX3g\n0b7tGnhIqgwqwAEgAIBFFQLHp3PHTJybdKHIVNaida0/4uXl77/UX0JhAABAPWwNHKc/Hjog\nYUNQ16EfTp/wQHiInzb/xE+Hls6ZtiphQGH46Y8GNJVZJOyAyQYAgIJsDRyLntni3WjIDzs+\n8tdpLC0BQc0e6PaouWmT/5uySAxYKq1CAADg9GwNHMmX80MTni9LGxYand8Lk1ttTNgoBIED\ncBpMdwGoebbe+MvLxaXwYuFf24suFblove1aEgAAUBtbA8eUFr4nkmL3Xr0lcxRdPzhqxc++\nLaZIKAwAAKiHrYdURn0yc3bbZ3rcHTZs8pgHWrXw0eSeTD/y/pI1vxXq3/w4VmqJAADA2dka\nOPzDJ//8db3Jz8YnLUhI+qOxQYe+a956Z1i4v6TiAACAOlThPhxBDw399JshV86eOnHixHWz\nT8uWLZs3qX+nT5sFAAC1QFWfpaIJaNwioHELKbUAAKqOy47gFGyaocj+LXXb91fLVguzPhsX\n/0ryF/tyTGZphQEAAPW4TeAovpE2ITLCv2n7GVvOlDWain5d+cbMp3p1CW7X7/Ofb0iuEAAA\nOL3KAkdJ3rEuLe5dvvOnB/rGxvdpUtbu2SD2250psyb0yv1xS7/27XZfL5JfJwAAcGKVBY59\nUwcduVo4Ounb/SnvD2lfr6xd4+LRsXu/mUs/O7413lj4y5jJ++TXCQAAnFhlgeP1lN+8G45b\n+XQHaxs07blwVnO/s1+8KaEwAACgHpUFjr03iho8NKTyn498uEFR9n67lgQAANSmssDRQO9S\nWlJa+c8XXyt20fnatSQAAKA2lQWOAfU9Lx1MqmQDIcSqA5fc6/aya0kAAEBtKgscsbM65V38\ncPiaH6xt8PP62HWX8tpOHS+hMAAAoB6VBY6WT38aHeq3LrZjvwlzvzuXX/6lkpyTb8U/1TFm\njUf9bp8+GyG5SAAA4Nwqu7W5xsUz6fv92p5Ra99N2LxiXnBISLNmzep7lJzKyMg4fup6scmn\neY9P9m0OdOWBKgAAoDK3eZaKq1frNXt/GZG8bMUH/0pNP7Hvq7SSUrPPXXeH3W/4x4DYuImD\n6ulIGwAA4DZseHibRtf9qSndn5oihCgtyb1a4Frfx016XQAAQEWq9rRYF1fv+q6SKgGA2osn\nvkL1OCACAACkq9oMR/WsGDnYa0HSsLs8y1oOJy9J3n30bI42rHWH4ZNGh9ZxrbwdAAA4NcmB\nw1ycumPV51cLBpVry0hOmLfx1PCJk8L9jZ+vWPpyXMG6lXFajdV2AKhJlR/dWL9+fY1VAqiJ\nxMBxYd+iZ9/an1d8683RzcWLPk5rEf3GwMjmQoiQhS5PxixYc25sbJC+4vbG3vIqBAAANUPi\nORwB7WIWJL79duJL5RsLr+/KLDZF9WhkWXXz69zeW5+684K1dnnlAQCAGiNxhkPn3eBub2Eq\nviXTlOSnCSHCPW+enBHuqduadqOkR8XtZauJiYm7d++2LPv6+q5du1Ze5QAAwL5q4qTR8kqL\n8oUQAeVuFxbgqjXlFVlrL1vNyso6d+6cZTk//5b7rAMAAAdX04HDRe8hhMgylnpptZaWqyUm\nrZ/eWnvZD3bt2rVBgwaWZQ8PjxotGgAA3JmaDhyuXm2E2HO8wNjE7fdg8UuhySfC11p72Q9G\nRUVFRUXVcLUAAMAuavrGX25+hkC9dvv+S5ZVY0H6oZzidpGB1tpruDwAACBDTQcOjUYfPyDi\nRNLsHd+nZ576cXXCXI/G3Uc0rmOtvYbLAwAAMtT0IRUhROiQuc+Lxckr56/I1YVFdE18bpTl\n7l7W2gEAgLOTHji0+sabN2/+U2OnIVM7DalgY2vtAADAqfHwNgAAIB2BAwAASEfgAAAA0ilw\n0igAqFLlj5kFajlmOAAAgHQEDgAAIB2BAwAASEfgAAAA0hE4AACAdAQOAAAgHYEDAABIx304\nAKAKuNkGUD3McAAAAOkIHAAAQDoCBwAAkI7AAQAApCNwAAAA6QgcAABAOgIHAACQjsABAACk\nI3AAAADpCBwAAEA6AgcAAJCOwAEAAKQjcAAAAOkIHAAAQDoeTw+g1uER80DNY4YDAABIxwwH\nABViDgNwNMxwAAAA6QgcAABAOgIHAACQjsABAACkI3AAAADpCBwAAEA6AgcAAJCOwAEAAKQj\ncAAAAOkIHAAAQDr139qcOxwDAKA4ZjgAAIB0BA4AACAdgQMAAEhH4AAAANIROAAAgHQEDgAA\nIB2BAwAASEfgAAAA0hE4AACAdAQOAAAgHYEDAABIR+AAAADSETgAAIB0BA4AACAdgQMAAEhH\n4AAAANIROAAAgHQEDgAAIB2BAwAASEfgAAAA0hE4AACAdAQOAAAgHYEDAABIR+AAAADSETgA\nAIB0BA4AACAdgQMAAEhH4AAAANIROAAAgHQEDgAAIB2BAwAASEfgAAAA0hE4AACAdAQOAAAg\nHYEDAABIp1O6AACopqFDhypdAgBbMcMBAACkc6wZjsPJS5J3Hz2bow1r3WH4pNGhdVyVrggA\nANiBA81wZCQnzNt46MH+Y2Y+E1Pn9O6X45aYzErXBAAA7MFhAoe5eNHHaS2i5wyM7BTRscvU\nhZPzL3695lyu0mUBAAA7cJRDKoXXd2UWmyb0aGRZdfPr3N5bn7rzgogJsbRkZWUVFBRYll1c\nXBo2bKhMoQAAoOocJXCU5KcJIcI9b560Ee6p25p2o2w1MTFx27ZtlmV/f/+vvvqqhisEAADV\n5iiHVEqL8oUQAbqb9QS4ak15RcpVBAAA7MZRZjhc9B5CiCxjqZdWa2m5WmLS+unLNpgwYUJ0\ndLRlWfvHNrZYv369/coE4EAY3YATcZTA4erVRog9xwuMTdx+DxO/FJp8InzLNmjUqFGjRo0U\nqg4AANwRRzmk4uZnCNRrt++/ZFk1FqQfyiluFxmobFUAAMAuHCVwaDT6+AERJ5Jm7/g+PfPU\nj6sT5no07j6icR2l6wIAAHagMZsd6O5aBzcsTt599HyuLizi/onPjQrUO0oeAgAAd8KxAgcA\nAFAlphAAAIB0BA4AACAdgQMAAEhH4AAAANIROAAAgHS1K3DExMQYDIZVq1YpXYjdbNq0yWAw\n9O/fX+lC7CY3N9dgMBgMhgMHDihdi93Mnj3bYDDMmDFD6ULs5ptvvrF007Vr15SuRQghDhw4\nYKknNzdX6Vrspn///gaDYdOmTUoXYjerVq0yGAwxMTFKF2I3mZmZlh3v2LFjStdiN/Hx8QaD\nYf78+XZ/Z0e5tXnNyMvLy87OLipSzzPhioqKsrOzdTr19KPZbM7OzhZCGI1GpWuxm4KCguzs\n7IKCAqULsRuj0WjpJge5rt7R6rGL3NxcVf6+ysvLU7oQuyn7fWUymZSuxW7y8/Ozs7MLCwvt\n/s61a4YDAAAogsABAACkU89UvC0efPDBli1bhoSEKF2I3QQHB0dGRnp7eytdiN3odLrIyEgh\nREBAgNK12E3btm2FEOHh4UoXYjf16tWzdJNer1e6FiGECAgIsNSjpsOLXbt2zc3NDQ4OVroQ\nuwkJCYmMjKxfv77ShdiNu7u7Zcfz8/NTuha76dChg4+PT0REhN3fmVubAwAA6TikAgAApCNw\nAAAA6dRzvPO2DicvSd599GyONqx1h+GTRofWcVW6ouooNV7ZvGrF9m/SL2eXNmoa2mfY2Mi/\nBQohLh6cMWb+D+W3jE3a1K+eu0JlVkEllTtpl+WceyP6n7v/1Kj3uueTDXOctJtWjBzstSBp\n2F2eZS3WukapLnPSXaU89Q1tobrRzdC+w/6qLYEjIzlh3sZTwydOCvc3fr5i6ctxBetWxmk1\nSpdVdVtnTVub4T9malxLf5fUHRuWJEwyLVvTM8jreup1j3q9p465eZpPsDMMYCGEtcqdt8s8\n6/Z+8cVO5VsOrX77RESUsP5hHZe5OHXHqs+vFgwq12ata5TqMufdVcpT39AWqhvdDO077S9z\nbVBaNHZAv2c/PmlZK7y2r3fv3qvO5ChbVDUYi8707dPntW8v/9FQuiTmyaefP2A2m/dMHjZ2\nYZqCtVVbxZWrpcvMZnPWj+sGRidcLSk1O1s3Ze5d+NSAvr179+7du/eHF/N+b7XWNUp1mSp2\nFVUObbPaRzdDu6pqxTkchdd3ZRabono0sqy6+XVu761P3XlB2aqqwViQ0bRZs8fDy66/0rTz\ncSvJzhNCHLtR5N/ez1SQfeHSdee67qjCylXTZWbjtfmv/GvwnGl1dRrhbN0U0C5mQeLbbye+\nVL7RWtco1WXq2FVUObSFqkc3Q7sa/26tOKRSkp8mhAj3vDnBFe6p25p2Q7mKqsnNt9tbb3Ur\nWy289P3q87nBsWFCiO/zSkr3vT1oyc8lZrPO864+I58Z0bONYoVWRYWVq6bLTqXMy2w4eGCz\nOpZV5+omnXeDu72FqfiWryXWuqakhzJdpo5dRZVDW6h6dDO0q/PvVrNep1JalC+ECNDd/M8N\ncNWa8pz7CQWnDm9Z9PrqkuCo6Y80NhWfyzKJZn73z1n1Un23oiNbk15fNsOj+drBLX2VLvM2\nrFX+iFYNXWYqPjc/OWPI0tllq07aTeVZG01KjTL1jW51DG2h6tHN0K5ef9WKwOGi9xBCZBlL\nvbRaS8vVEpPWzyHukFgNxdknkxLf+CI166Enxk8cFuXpohHaoJSUlD9er9Nl0LTjW7/9/N1j\ngxO7KFmoDbT6iit/LE4NXXb+y8U5daIea/D7GeDWPqzjd1N51kaTUqNMTaNbTUNbqHp0M7Sr\n11+14hwOV682QojjBTefPvpLocknwpmyZ5n8zL2TRsUfLm6z8L0PpsX09HSp+ETh9g08SnKv\n1HBtdmGpXBVdZv5w46mQYX0q2cIZu8la1yjVZarYVYSoBUNbqGd0M7Sr2V+1InC4+RkC9drt\n+y9ZVo0F6YdyittFBipbVXWYjfPjF7t1H7/y1QlhATev8L5xcmX0sFHnisuej2zeez7fNzxU\nkRqrxFrlKuiygisph3NKRna+WbPzdlN51rpGqS5Twa4ihAqHtlDv6GZoV7u/asUhFY1GHz8g\n4oWk2TsaTmvtZ9yy7DWPxt1HNK6jdF1Vln/xw9Sc4qfb1vnuyOGyRp1HaIeI6Jba/0yfuXzi\n0Mi62sLv/vPR1/k+C8aHKViqjXyaVly5RqNz9i47v3WPq3fHEHdtWYu1D6tgkdVgbTRpNEKR\nLlPH6Fbf0BbqHd0M7Wr3Vy16eNvBDYuTdx89n6sLi7h/4nOjAvXON7tzYd+MsQt/+FOjT5Pp\nHy19oOjaD6uXrT30v1/zhHfzkHuix4+5p6FnhW/iaCqp3Km7bHXs4H1B8avn3Fu+0Rm7yVR8\n9omBEwa9v7H87QitdY1SXebUu4pQ6dAWKh3dDO1q91ctChwAAEApThMqAQCA8yJwAAAA6Qgc\nAABAOgIHAACQjsABAACkI3AAAADpCBwAAEA6AgcAAJCOwAEAAKQjcAAAAOkIHAAAQDoCBwAA\nkI7AAQAApCNwAAAA6QgcAABAOgIHAACQjsABAACkI3AAAADpCBwAAEA6AgcAAJCOwAEAAKQj\ncAAAAOkIHAAAQDoCBwAAkI7AAQAApCNwAAAA6QgcAABAOgIHAACQjsABAACkI3AAAADpCBwA\nAEA6AgcAAJCOwAEAAKQjcAAAAOkIHAAAQDoCBwAAkO7/AUCod7IPlLWDAAAAAElFTkSuQmCC\n", "text/plain": [ "plot without title" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "model <- stan_model(model.file)\n", "fit <- sbc(model, list(x=x, n=n))\n", "plot_fit(fit)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And another pathological example:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "data {\n", "\tint n;\n", "\tvector[n] x;\n", "}\n", "\n", "transformed data {\n", "\treal beta_sim = normal_rng(0, 11);\n", "\treal sigma_sim = lognormal_rng(0, 10);\n", "\t\n", "\tvector[n] y_sim;\n", "\tfor (i in 1:n)\n", "\t\ty_sim[i] = student_t_rng(10, beta_sim, sigma_sim);\n", "}\n", "\n", "parameters {\n", "\treal beta;\n", "\treal sigma;\n", "}\n", "\n", "model {\n", "\tbeta ~ normal(0, 1);\n", "\tsigma ~ lognormal(0, 1);\n", " \ty_sim ~ normal(x * beta, sigma);\n", "}\n", "\n", "generated quantities {\n", "\tint idsim[2] = { beta < beta_sim, sigma < sigma_sim };\n", "}\n" ] } ], "source": [ "model.file <- \"_models/sbc_3.stan\"\n", "cat(readLines(model.file), sep=\"\\n\")" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtAAAAFoCAIAAADxRFtOAAAABmJLR0QA/wD/AP+gvaeTAAAf\nrklEQVR4nO3deUCU1f7H8TPMMOyLikqCqaggogV6s9S86UBmi2a55YIprte1ULOfZVqauZSl\npqkVZClJm13NzLpq7kuLl6KbCKLljkLKzqy/P8YQS3BG5/Aww/v118x3Hma+D2cOfObMM/Oo\nLBaLAAAAkMlN6QYAAIDrI3AAAADpCBwAAEA6AgcAAJCOwAEAAKQjcAAAAOkIHAAAQDoCBwAA\nkI7AAQAApCNwAAAA6QgcAABAOgIHAACQjsABAACkI3AAAADpCBwAAEA6AgcAAJCOwAEAAKQj\ncAAAHMZYnB4SEvLQ8l+VbgQ1jkbpBgAArkPlXm/48OEhbesq3QhqHJXFYlG6BwAA4OJ4SwU1\nRaewJm9mp00f1k/XPqpLj76bjhekfTy354PdI1u1fWTUi/lGixBCWPQhISFLzxSW/1RM09un\nHL+sWNNALXZ219r4h7q2btG0TbuOY2avLjNfqcc0vX3CsUtCCGNR1uyR/Tu0admh28Oz1uxb\n2KFV/60nhY2TXYiyPw4+O6JXTJuIpi1b39v9iVVf/67QjsIxCByoQdaNWDbw1Q+2f582pvHp\nCbH/fD6zzWdffn1493s525Im7DyjdHcArjIUHHogfoap86AVH25YMO3JQ+/NjX8/89pNTFN1\nj2w4FzLnrfVLZo7OXpmw8uzVlwq2TPaXew3/8lzEc4ve2bBmSf9/FMwd9cBJvakadxEOxjEc\nqEFCnp4ZHeQphOgxqdWzW3d98Myj7irhXr/94AZeX3yfK2JDlG4QwBWll7blGk0jxwztGuQp\n7oppWTf0mE9AxQ1yf57x6TnVlu2L2vhohGh/x4aS8PZTym+1ZbI37Dd20eDRD9TzFEK0jlIv\nWBP/c5GhsVZdvTsKhyFwoAapE+lvvaD21Ki1twdqVNar3m4qYa78xwBUO5/bxvRq/eGwf/yj\ny8MP33v3PY881vNBv2v+oZza9INn3Yfb+Fwp+gQ/EeIxvfxWWyb7+Injfj186NMt//v1118O\n7d0qe48gG2+pwJlZjIVmjnoGFOCmqfvW1z9sXDO/fUOxfd0rndtGT3jjQMUNLHqL6tp/Me4q\nO+7fbLw4vX+n3uNe3JWRF96x57yVSxzSNhTECgeci0oIkWe4EjIKT6cUm1j6ABRw8Yfk1d+U\nzHh2XPR9jwgh0t5++NFFs9586uo6RKOHWpUmf3m0ZF64l0YIUXLxixOlRtvfFr2U9cK6AzmH\nju1vpFULIYrPr3P4LqCascIBp6Jyvy/AY+PUJT8dO3Xk+21Px69q7EFoBhSgDfxjxZvzJr25\nbt/hX3Zv/WzNJ78HtHiw4gYN/vFa97plA4Y8t31/2o97vpo8YHaYp0altnWVwyMw2mI2rNu0\nM/v07//d9dn4x18TQuzce7iYRU2nReCAk1mWMrdF3hd94uLGv/h219c39u98VxMPDiIDqpt/\n88T1s0ZnfPj6k4899K9nF16IGvzRJxOu2cLNc8WOz7p7H5kyvO/kue/c9dKm2EBPj3paG+/f\nJ3h08vTBG16eGBf76AtvbRvy3rcrEh7eu2jWyTI+qOKs+OIvOAFT6eXzRZ6N6nko3QgAWxmL\n/7dt90ld9wesh26YjbmdwtsN2/u/sbf5VPFTTHYXxgoHnIDaM4A/QIBzsZhLp4weNX715t8u\nFBVcPPnucwNyPDuMCK4qbQgmu0tjhQMAIMXZb9+dOHv5d1k5FnffsLa6l5Yt/GcTX6WbgmII\nHAAAiSzGUrPGkyOtQOAAAADSSflIobHoxNo3397zv+xLBu3tzWMGjx3dPsTbetPB1GWpOw+f\nKlBHtG4XP2FkuJ/7zdUBAIATkbLCsXbykM36qImje9XXln2zdsn2E43eWTsvUK3KSp05ZX12\n/PgJkXWMm1ct/9HSYd3qRLVK2FsHAADOxfGfUtHn7/voeP7jsyZ3iolqGdVu5IyJ+qJfUi8U\nC4t+0cfpzQfP6RvXMap9l8kLJxaf/3bN6UK76wAAwNk4PnBYLKVdunTpVs/TelXt0UgIYTBb\nSi/tOKs3dY9tZK17BHaO8dWmbT9nb93hDQMAANkcfwyHR4Bu2jSdEEKfl3M278LBL1dq/VsP\naehjOJcuhIj0vnoQRqS3Zkv6ZUOsffWKjzVs2DCT6cq3zsXGxg4bNszhuwMAAG6dxPNQpL00\nZU72ZZXKvfeUVwPVqstlxUKIIM3VNZUgd7WpqMxsZ73iQxw5csRoNFovt2nTRt6+AACAWyEx\ncNz1xgcbhcjJ+Dbx/xItQUn9/L2EEHlGs4/6yuexcw0mdaDWTWtfveJDDB06tHyFo23btvL2\nBQAA3ArHB47LGdt3HPXs3bOT9WqDiK696q78KuXEoClthNh1tMTY+M9TbZ0oNflHBbj72Fev\n+Fjjxo1zeP8AAMDhHH/QqFG/PzlpWY7BfOW6xfhzsdEz2MsjUBesVW/dm3Nls5KMAwX66Lhg\ne+sObxgAAMimnj17tmPv0bNe1E8bP//3kfzQer7Fuae/Tnl1e7YpcdbQ2zy1Eab01JTN9VtE\neJWcT10477R3p5eeuM9Npbaz7th+AQCAdFK++Kv41A+rV6X8N/tkkfBu0uyOvgkj7wnzt960\n/8MlqTsPnynURETdPX7KiGCt283VAQCAE+FcKgAAQDoWDAAAgHQEDgAAIB2BAwAASEfgAAAA\n0kn8plEAqB6DBg2q4taUlJRq6wRAZVjhAAAA0hE4AACAdAQOAAAgHYEDAABIR+AAAADSETgA\nAIB0BA4AACAdgQMAAEhH4AAAANIROAAAgHQEDgAAIB2BAwAASEfgAAAA0hE4AACAdAQOAAAg\nHYEDAABIR+AAAADSETgAAIB0BA4AACAdgQMAAEhH4AAAANIROAAAgHQEDgAAIB2BAwAASEfg\nAAAA0hE4AACAdAQOAAAgHYEDAABIR+AAAADSETgAAIB0BA4AACAdgQMAAEhH4AAAANIROAAA\ngHQEDgAAIB2BAwAASEfgAAAA0hE4AACAdAQOAAAgHYEDAABIR+AAAADSETgAAIB0BA4AACAd\ngQMAAEhH4AAAANIROAAAgHQEDgAAIB2BAwAASEfgAAAA0hE4AACAdAQOAAAgHYEDAABIp1G6\nAQDOZ9CgQVVvkJKSUj2dAHAWrHAAAADpCBwAAEA6AgcAAJCOwAEAAKQjcAAAAOkIHAAAQDoC\nBwAAkI7AAQAApCNwAAAA6QgcAABAutr41eZ8KzMAANWMFQ4AACAdgQMAAEhH4AAAANJJOYbD\nbLy48d1VW7/LuJBvbtQ0vNeQ0XF3BFtvOpi6LHXn4VMF6ojW7eInjAz3c7+5OgAAcCJSVji2\nzJ72/o7cXiMS57/0zH1NSpbNnLD1dJEQIit15rz1Bzo9PmrWU0P9ju98IXGZySJuog4AAJyL\n4wOHSX/q7Z/z7p4248GO0S1a3dFn/Lz7A90+XPqTsOgXfZzefPCcvnEdo9p3mbxwYvH5b9ec\nLrS7DgAAnI3jA4exJKtps2aPRAb+WVBF+3sY8otKL+04qzd1j21krXoEdo7x1aZtP2dv3eEN\nAwAA2Rx/DIdHQNc33uhafrU058ekM4VNEiIMxR8JISK9rx6EEemt2ZJ+2RCbble94mPNnz/f\nbDZbL7dr165Hjx4O3x0AAHDr5H7xV/bBTYteTTI06T6jR6jpt2IhRJDm6ppKkLvaVFRmLrOv\nXvH+P//8c6PRaL2sVqsJHAAA1EyyAoc+/1jy4te+TMu797Gx44d093ZTFWi9hBB5RrOPWm3d\nJtdgUgdq3eysV3yUDh06lAeOJk2aSNoXAABwi6QEjuKzu5+atNjY8v6FbydEBHlai+4+bYTY\ndbTE2NjjSoA4UWryjwqwt17xgZYuXSqjfwAA4FgSPhZrMb4ydYlHt7GrXx5XnjaEEB6BumCt\neuveHOtVY0nGgQJ9dFywvXXHNwwAACRz/ApH8fkP0gr0T7b1++HQwasP4xXe/o46U/tETU9+\ncdtt01oHGjetmO8V2m1YqJ9KJeyqO7xhAAAgm+MDR35WlhBizaL5FYv+jWesXX5P+MC5z4gl\nqatfWVWoiYj65+IpI9QqIYSwtw4AAJyL4wNH8L0vb7y30ls7DpzccaAD6gAAwIlw8jYAACAd\ngQMAAEhH4AAAANIROAAAgHQEDgAAIB2BAwAASEfgAAAA0hE4AACAdAQOAAAgHYEDAABIR+AA\nAADSETgAAIB0BA4AACAdgQMAAEhH4AAAANIROAAAgHQEDgAAIJ2tgaNjx46vnir8e/3cvkld\ndPEObQkAALgaTdU3Z2RkWC8cOHAg7NdfM4r8r7nZYjz0+bf7dv8uqTkAAOAabhA4WrVqVX45\npXuHlOttExA2yaEtAQAAV3ODwLFy5UrrhbFjx9435/WB9b3+soHa3b9zv75SWgMAAK7iBoFj\nzJgx1gvr16/vnTByTCNf+S0BAABXc4PAUW7Hjh1CiLzT2RcKDX+/NbRluI+bypF9AQAAF2Jr\n4Ci9uK1vl/6bj+Rd99b/Furv9HF3XFcAAMCl2Bo4Vj86ZOtvAZNmTmsdUufvt0Z5kzYAAECl\nbA0cc7+7kPDNySX33Sa1GwAA4JJs/eIvb7VqwJ11pbYCAABcla2BY2aHBu/sOCu1FQAA4Kps\nDRzxX2z84/8efOX9bcUmi9SGAACA67H1GI7Yh54yBRhmPBn33HCPhiHBnuprPgR7/PhxCb0B\nAAAXYWvgCAoKEiKod++2UrsBAAAuydbAsWHDBql9AAAAF2Zr4Lh8+XJlN6nUXv6+Wgf1AwAA\nXJCtgSMwMLCym+q0WJGX+S8H9QMAAFyQrYFj9uzZFa+aywqOZ6Z/vuGbwK4TlyR2d3xfAADA\nhdgaOGbNmvX3YuGJ/3SIeuirosmPObQnAADgYmz9Ho7r8m0a99GcmLUTnndUNwAAwCXdUuAQ\nQvg08SnN+8IhrQAAAFd1S4HDbMhZ/Px/NV7hjuoGAAC4JFuP4ejYsePfaqbTR386mVd218w3\nHdsTAABwMbYGjutRN7kz7vHYIQtn3O2wdgAAgCuyNXDs379fah8AAMCF2bfCUZabvvHLA1lZ\nmRcMPuHh4fc8+Gh0Qy9JnQEAAJdhR+D4bO6o8XOTz5WZyitq9/rDXlj5zvOPS2gMAAC4Dls/\npXL840F9Zr6jvnvAB1/tyfzt3IVT2fu+SXmio/bdmX2GfHpCZocAAMDp2brCseipTb6NBv68\nbW0djcpaCQppdk/XBy1NG/970iLRZ7m0DgEAgNOzdYUj9UJx+NhnytOGlUoTOH1iq+IL6yU0\nBgAAXIetgcPHza30fOnf62U5ZW5qX4e2BAAAXI2tgWNS84DM5ITduddkjrJL+0esOhLQfJKE\nxgAAgOuw9RiOEZ/MerHtU7G3RwyZOOqeVs39VYXHMg69s2zN76Xa1z9OkNoiAABwdrYGjjqR\nE498W2/i01OTF8xM/rPYsN2ja954c0hkHUnNAQAA12DH93CE3Dvos+8GXjyVnZmZecni37Jl\ny7DG9W/1bLMAAKAWsPdcKqqg0OZBoc2l9AIAAFyUTSsU+b+nffVjbvnV0rwvxkx9KfXLPQUm\ni7TGAACA67hB4NBfTh8XF1Wnacxzm06WF01lv61+bdYTD3dpEt1785HLkjsEAABOr6rAYSj6\nqUvzu1Zu//WeRxOm9mpcXvdumPD99g2zxz1c+Mum3jHROy+Vye8TAAA4saoCx57J/Q/llo5M\n/n7vhncGxtQrr6vcvNp36z1r+RdHt0w1lp4YNXGP/D4BAIATqypwvLrhd9/bxqx+sl1lGzR9\nYOHssMBTX74uoTEAAOA6qgocuy+XNbx3YNU/H3dfw7L8vQ5tCQAAuJqqAkdDrZvZYK765/V/\n6N00AQ5tCQAAuJqqAkef+t45+5Or2EAI8e6+HM+6Dzu0JQAA4GqqChwJszsWnf8gfs3PlW1w\nJCVhXU5R28ljJTQGAABcR1WBo+WTnw0OD1yX0L73uLk/nC6ueJOh4NgbU59oP3SNV/2unz0d\nJblJAADg3Kr6anOVm3fyj3vVD3R//62ZG1fNa9KiRbNmzep7GbKzsrKOZl/Sm/zDYj/ZszHY\nnROqAACAqtzgXCruPq3X7D4xLHXFqvc+TcvI3PNNusFs8W9we8Tduvv7JCSO719PQ9oAAAA3\nYMPJ21Sabk9M6vbEJCGE2VCYW+Je399Del8AAMCF2He2WDd33/rukjoBAAAuizdEAACAdAQO\nAAAgHYEDAABIR+AAAADSETgAAIB09n1KxV6rhg/wWZA8pIF3eeVg6rLUnYdPFagjWreLnzAy\n3M/95uoAAMCJSFvhsOjT/vPW5tySirWs1Jnz1h/o9PioWU8N9Tu+84XEZSbLzdQBAIBzkbLC\ncW7Poqff2Fukv/bU9hb9oo/Tmw9+rW9cmBCixUK3fkMXrDk9OiFEa1891FdGzwAAQB4pKxxB\n0UMXLF66dPHzFYull3ac1Zu6xzayXvUI7Bzjq03bfs7euoyGAQCAVFJWODS+DW/3FSb9NWnG\nUJwuhIj0vnoQRqS3Zkv6ZUOsffWK9/n555+bzVfWUcLCwqKjox2/MwAA4JbJPWi0InNZsRAi\nqMLJ3oLc1aaiMnvrFe9z/vz5RqPRerlfv34EDgAAaqbq+1ism9ZLCJFnvHpgR67BpPbS2luv\ntoYBAICjVN8Kh7tPGyF2HS0xNvZQWysnSk3+UQH21ive54EDB6qtfwAAcNOqb4XDI1AXrFVv\n3ZtjvWosyThQoI+OC7a3Xm0NAwAAR6m+wKFSaaf2icpMfnHbjxlns39JmjnXK7TbsFA/e+vV\n1jAAAHCU6ntLRQgRPnDuM2JJ6upXVhVqIqL+uXjKCLXqZuoAAMC5SAwcam3oxo0b/1LsOHBy\nx4HX2djeOgAAcCKcvA0AAEhH4AAAANIROAAAgHQEDgAAIB2BAwAASEfgAAAA0hE4AACAdAQO\nAAAgHYEDAABIR+AAAADSETgAAIB0BA4AACAdgQMAAEhH4AAAANIROAAAgHQEDgAAIB2BAwAA\nSEfgAAAA0hE4AACAdAQOAAAgHYEDAABIR+AAAADSETgAAIB0GqUbAAAASho0aFDVG6SkpNz6\no7DCAQAApCNwAAAA6QgcAABAOgIHAACQjsABAACkI3AAAADpCBwAAEA6AgcAAJCOwAEAAKQj\ncAAAAOkIHAAAQDoCBwAAkI7AAQAApCNwAAAA6QgcAABAOgIHAACQjsABAACkI3AAAADpCBwA\nAEA6AgcAAJCOwAEAAKQjcAAAAOkIHAAAQDoCBwAAkI7AAQAApCNwAAAA6QgcAABAOgIHAACQ\njsABAACkI3AAAADpCBwAAEA6AgcAAJCOwAEAAKQjcAAAAOkIHAAAQDoCBwAAkI7AAQAApNMo\n3QAAAJBu0KBByjbACgcAAJCOwAEAAKQjcAAAAOkIHAAAQDoOGr2Oqo+sSUlJqbZOAABwDaxw\nAAAA6QgcAABAOt5SAXB9in9q31F4kxS1RA2fswQOAAAcT1LSreGpogq8pQIAAKRjhcNuVaRL\n1mYBALfIedcwqlbTA8fB1GWpOw+fKlBHtG4XP2FkuJ+70h1VhbeKAUCeG/4nruzPbPX/4A25\naqqoQo0OHFmpM+etz44fPyGyjnHzquUvJJasW52oVindFuA85P25dBk39zrhVv5bVP/vvPpf\nCzndf9Obbtjp9lRBNThwWPSLPk5vPvi1vnFhQogWC936DV2w5vTohFBfpTsDapZb+ZPHn8uq\nyfj93PS/f3LDDTldw7VKzQ0cpZd2nNWbxsU2sl71COwc46tN235ODG1Rvs2ZM2csFov1so+P\nT2BgoAKN2kzSCylneYXqMjtStVqym5Cn+l9q808a1aPmBg5DcboQItL76kEbkd6aLemXK27z\n+OOPG41G6+V+/fpNnz69Ojt0LP5YuMyOVK2W7CYA/EXN/VisuaxYCBGkudphkLvaVFSmXEcA\nAOAm1dwVDjetlxAiz2j2UautlVyDSR2orbjNe++9V/6WSp06dWy8Z9a0ARfDpAZqvpobONx9\n2gix62iJsbHHlcBxotTkHxVQcZtWrVop0RoAALBPzX1LxSNQF6xVb92bY71qLMk4UKCPjgtW\ntisAAHATam7gUKm0U/tEZSa/uO3HjLPZvyTNnOsV2m1YqJ/SfQEAALupyo+BqJn2f7gkdefh\nM4WaiKi7x08ZEaytuQkJAABUpqYHDgAA4AJYMAAAANIROAAAgHQEDgAAIB2BAwAASEfgAAAA\n0tX2wLFlyxadTnf//fcr3YjDPPPMMzqdbt68eUo34hiFhYU6nU6n0+3bt0/pXhxj06ZNOp2u\nR48eSjfiMImJiTqdbsGCBUo3IoQQ2dnZ1idMZmam0r04RlJSkk6nGzp0qNKNOEx8fLxOp0tK\nSlK6EcfIzMy0PuWys7OV7sUx3n77bZ1ON3z4cIffc839avPqYTAY8vPz1X+ersUFlJSU5Ofn\nl5SUKN2IY1gslvz8fCFE+WmBnZ31Kefh4aF0Iw5jfcqVlpYq3YgQQpjNZusTxmw2K92LY5SV\nleXn5xcWFirdiMMUFRXl5+fr9XqlG3EMV33KFRUVOfyea/sKBwAAqAYEDgAAIF1tf0ulUaNG\ncXFxbm6uE7xiYmJ8fHyioqKUbsQxNBpNXFycECIoKEjpXhwjNDQ0Li7O3d1d6UYcpl27dgEB\nAa1bt1a6ESGE8PX1tT5h/Pxc5LxLzZs3j4uLc5nnvxCiU6dO4eHhYWFhSjfiGH5+ftannK+v\nr9K9OEaLFi3i4uIaNGjg8Hvmq80BAIB0rvPKHgAA1FgEDgAAIF1tP4bjYOqy1J2HTxWoI1q3\ni58wMtzPyd5ZNxsvbnx31dbvMi7kmxs1De81ZHTcHcFCiPP7nxv1ys8Vt0xI/qh3PU+F2rRD\nFZ0742AVnH5t8L92/qWo9bnzkw/nOOMYrRo+wGdB8pAG3uWVygZFqcFyxifJX7jYpGZGV2N3\nN6M6J3WtDhxZqTPnrc+OHz8hso5x86rlLySWrFudqFYp3ZY9tsye9n5WnVGTE1vWcUvb9uGy\nmRNMK9Y8EOJzKe2SV72ek0ddPXS0iTNMZiFEZZ076WB51+357LMdK1YOJC3NjOouKt/TGsqi\nT9v27ubckv4VapUNilKD5aRPkr9wsUnNjK65qn9SW2otc9noPr2f/viY9VrpH3t69uz57skC\nZZuyi7Hs5KO9es3//sKfBfOyof2efGafxWLZNXHI6IXpCvZ2067fufMPllXeL+v6Dp6ZazBb\nnGqMzu5e+ESfR3v27NmzZ88PzhddqVY2KEoNlks8SVxvUjOjayZFJnXtPYaj9NKOs3pT99hG\n1qsegZ1jfLVp288p25VdjCVZTZs1eyQy8M+CKtrfw5BfJIT46XJZnZhAU0n+uZxLzvUxpOt2\n7gKDJYSwGP945aVPB8yZVlejEk41RkHRQxcsXrp08fMVi5UNilKD5RpPEteb1MzomkmRSV17\n31IxFKcLISK9ry55RXprtqRfVq4ju3kEdH3jja7lV0tzfkw6U9gkIUII8WORwbxnaf9lRwwW\ni8a7Qa/hTw17oI1ijdrjup27wGAJIbI3zDt724C+za58P4QTjZHGt+HtvsKkv+b1SWWDYohV\nZrBc40niepOaGV0zKTKpa2/gMJcVCyGCNFd/3UHualNRmXId3ZLsg5sWvZpkaNJ9Ro9Qk/50\nnkk0C7x7zrvP1/coO7Ql+dUVz3mFvT+gZYDSbd5AZZ33UDv9YJn0p19JzRq4/MXyq046RuUq\nm0FKzSwXm9HCJSY1M7qGD9BfyJ7UtTdwuGm9hBB5RrPPn2duyzWY1IFaRZu6Gfr8Y8mLX/sy\nLe/ex8aOH9Ld200l1CEbNmz483a/Lv2nHd3y/ea3fhqwuIuSjdpArb1+5w8lOv1gnfl6SYFf\n94caXjkUvLI9rfljVK6yGaTUzHKZGS1caFIzo2v4AP2F7Elde4/hcPdpI4Q4WnL1HKQnSk3+\nUc6URoUQxWd3Txgx9aC+zcK335s29AFvt+sfNBzT0MtQeLGae3MIa+fOP1iWD9ZntxjSq4ot\nnG6MKhsUpQbL+Z8kV7j2pGZG12SyJ3XtDRwegbpgrXrr3hzrVWNJxoECfXRcsLJd2cdifGXq\nEo9uY1e/PC4i6OqnvS8fWz14yIjTelP5drvPFAdEhivSo10q69zZB6vk4oaDBYbhna827Lxj\nVK6yQVFqsJz9SXKFa01qZrQi7d002ZO69r6lolJpp/aJmp784rbbprUONG5aMd8rtNuwUGc6\n4VPx+Q/SCvRPtvX74dDB8qLGK7xd1OCW6v/MmLVy/KC4uurSH/6z9tti/wVjIxRs1Ub+Ta/f\nuUqlcerBOrNll7tv+xae6vJKZXuqYJP2qmwGqVRCkcFygRktXG5SM6MVbPImyJ7Utf3kbfs/\nXJK68/CZQk1E1N3jp4wI1jrTks+5Pc+NXvjzX4r+jWesXX5P2R8/J614/8D/fisSvmEt7hw8\ndtSdt3lf905qmio6d97BSkoYsCdkatKcuyoWnW6MTPpTj/Ud1/+d9RW/lLCyQVFqsJz3SWLl\nepOaGV2TVfOkru2BAwAAVAOnSZQAAMB5ETgAAIB0BA4AACAdgQMAAEhH4AAAANIROAAAgHQE\nDgAAIB2BAwAASEfgAAAA0hE4AACAdAQOAAAgHYEDAABIR+AAAADSETgAAIB0BA4AACAdgQMA\nAEhH4AAAANIROAAAgHQEDgAAIB2BAwAASEfgAAAA0hE4AACAdAQOAAAgHYEDAABIR+AAAADS\nETgAAIB0BA4AACAdgQMAAEhH4AAAANIROAAAgHQEDgAAIB2BAwAASEfgAAAA0hE4AACAdAQO\nAAAgHYEDAABI9/8Y+LH99OVmRgAAAABJRU5ErkJggg==", "text/plain": [ "plot without title" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "model <- stan_model(model.file)\n", "fit <- sbc(model, list(x=x, n=n))\n", "plot_fit(fit)" ] } ], "metadata": { "kernelspec": { "display_name": "R 3.5", "language": "R", "name": "ir35" }, "language_info": { "codemirror_mode": "r", "file_extension": ".r", "mimetype": "text/x-r-source", "name": "R", "pygments_lexer": "r", "version": "3.5.1" } }, "nbformat": 4, "nbformat_minor": 2 }