first draft of custom priors

docs/source/custom_priors.ipynb

@@ -0,0 +1,247 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Custom priors\n",
+    "\n",
+    "The prior probability is a critical element of Bayes theorem.\n",
+    "However, to keep `uravu` straightforward to use, by default, a broad uniform prior probability is assigned to the `Relationship` object. \n",
+    "\n",
+    "Of course this may be ignored and custom priors may be used (*and sometimes it may be necessary that this is done*).\n",
+    "This tutorial will show **how** custom priors may be used with `uravu`. \n",
+    "\n",
+    "Let's start as always be producing some synthetic data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.random.seed(2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "x = np.linspace(10, 50, 20)\n",
+    "y = .3 * x ** 2 - 1.4 * x + .2\n",
+    "y += y * np.random.randn(20) * 0.05\n",
+    "dy = 3 * x"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAASTklEQVR4nO3db4xcV3nH8e/TzR8MFDlunMisndpUVkqCC4tWJtQVahOo3YDw1hKKVVFZVSTzIrShQolsVSrtC8uRaCm8aFDcQLEKSrBCmlgIAZEDqoqiGAeHJI6JbBFw/KfxUpRCqyhO7Kcv9joaOzO7s3tn9s6c/X6k1cycvTvz+Hj3t2fPPXNPZCaSpLL8RtMFSJJ6z3CXpAIZ7pJUIMNdkgpkuEtSgS5pugCAK6+8MleuXNl0GZI0VJ544olfZObSdp8biHBfuXIlBw4caLoMSRoqEfHzTp9zWkaSCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkhpyyz2Pccs9j/XluQ13SSqQ4S5JBTLcJalAXYV7RPx1RByKiGci4r6IeFNELImIRyLiSHV7Rcvx2yPiaEQ8FxHr+1e+JKmdGcM9IkaBvwLGM/NdwAiwGdgG7MvM1cC+6jERcV31+euBDcDdETHSn/IlSe10Oy1zCbAoIi4B3gycBDYCu6vP7wYmqvsbgfsz85XMfB44CqztXcmSpJnMGO6ZeQL4B+AYcAr4n8z8LnB1Zp6qjjkFXFV9ySjwQstTHK/aLhARWyPiQEQcmJycrPevkCRdoJtpmSuYGo2vAt4OvCUiPj7dl7Rpyzc0ZO7KzPHMHF+6tO1GIpKkOepmWuaDwPOZOZmZrwIPAr8PvBgRywCq29PV8ceBFS1fv5ypaRxJ0jzpJtyPATdExJsjIoCbgMPAXmBLdcwW4OHq/l5gc0RcHhGrgNXA/t6WLUmazox7qGbm4xHxAPAj4DXgILALeCuwJyJuZeoXwMeq4w9FxB7g2er42zLzbJ/qlyS10dUG2Zn5GeAzFzW/wtQovt3xO4Ad9UqTJM2V71CVpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6Q56ucG13UZ7pIWrEEO57oMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkNeOjgCQ4ee4nHn/8l6+56lIcOnujp8xvukjTPHjp4gu0PPs2Zs+cAOPHSy2x/8OmeBrzhLknz7LPfeY6XXz17QdvLr57ls995rmevYbhL0jw7+dLLs2qfC8NdkubZ2xcvmlX7XBjukjTP7lh/LYsuHbmgbdGlI9yx/tqevcYlPXsmSVJXJsZGAbjzgac4c/Yco4sXccf6a19v7wXDXZIaMDE2yn37jwHw9U+8v+fP77SMJBXIcJekAhnuklQgw12SCmS4Sxpat9zzGLfc81jTZQwkw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVqKtwj4jFEfFARPwkIg5HxPsjYklEPBIRR6rbK1qO3x4RRyPiuYhY37/yJUntdDty/wLw7cz8XeDdwGFgG7AvM1cD+6rHRMR1wGbgemADcHdEjLR9VklqSN09TPu9B2pdM4Z7RLwN+ADwJYDMPJOZLwEbgd3VYbuBier+RuD+zHwlM58HjgJre124JM1V3T1M52MP1Lq6Gbm/A5gE/jUiDkbEvRHxFuDqzDwFUN1eVR0/CrzQ8vXHq7YLRMTWiDgQEQcmJydr/SMkaTbq7mE6H3ug1tVNuF8CvBf4YmaOAf9HNQXTQbRpyzc0ZO7KzPHMHF+6dGlXxUpSL9Tdw3Q+9kCtq5twPw4cz8zHq8cPMBX2L0bEMoDq9nTL8Stavn45cLI35UpSfXX3MJ2PPVDrmjHcM/O/gBci4vzmfjcBzwJ7gS1V2xbg4er+XmBzRFweEauA1cD+nlYtSTXU3cN0PvZAravbbfb+EvhaRFwG/BT4C6Z+MeyJiFuBY8DHADLzUETsYeoXwGvAbZl5tv3TStL8q7uH6XzsgVpXV+GemU8C420+dVOH43cAO2rUJUl9VXcP037vgVqX71CVpAIZ7pIa42Yb/WO4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pLmzKWMg8twl6QCdXttGUlSj/XzsgWO3CWpQIa7pKE06HuYNs1wlzR0hmEP06YZ7pKGzjDsYdo0w13S0BmGPUybZrhLGjrDsIdp0wx3SUNnGPYwbZrr3CUNnWHYw7RphrukoTToe5g2zWkZSSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIK5Dp3aQE7v0VeE+vEz1+y98zZc6y769FG3oRU8vp4w13SvOt0yV5gqN5lOsi/HJyWkTTvvGRv/xnukuadl+ztP8Nd0rzzkr39Z7hLmpM6e5h6yd7+84SqpFmre0LUS/b2n+EuDbGmljJOd0K024D2kr395bSMpFnzhOjgM9wlzZonRAef4S416JZ7Hnt9amWYeEJ08HUd7hExEhEHI+Kb1eMlEfFIRBypbq9oOXZ7RByNiOciYn0/CpfUnImxUXZuWsNlI1MRMrp4ETs3rfGE6ACZzcj9duBwy+NtwL7MXA3sqx4TEdcBm4HrgQ3A3RExgqSiTIyNMnbNYt63agk/2HajwT5gugr3iFgOfBi4t6V5I7C7ur8bmGhpvz8zX8nM54GjwNrelCtJ6ka3I/fPA3cC51rars7MUwDV7VVV+yjwQstxx6u2C0TE1og4EBEHJicnZ124JKmzGcM9Ij4CnM7MJ7p8zmjTlm9oyNyVmeOZOb506dIun1qS1I1u3sS0DvhoRNwMvAl4W0R8FXgxIpZl5qmIWAacro4/Dqxo+frlwMleFi1Jmt6MI/fM3J6ZyzNzJVMnSh/NzI8De4Et1WFbgIer+3uBzRFxeUSsAlYD+3teuSSpozqXH7gL2BMRtwLHgI8BZOahiNgDPAu8BtyWmWc7P40kqddmFe6Z+X3g+9X9/wZu6nDcDmBHzdokaVpek6Yz36EqSQUy3CWpQIa7JBXIcJekAhnu0gJVZ5s8DT7DXVqAOm2TZ8CXw3CXFqDptslTGQx3qYZh3WzDbfLKZ7hLC5Db5JXPcJeGVJ0Tom6TV74615aR1JBOJ0SBrnZEOn/MnQ88xZmz5xhdvIg71l/rbkoFceQuNaTOyLsXJ0TdJq9sjtylBtQdeZdyQtQLf/WPI3epAXVH3p4Q1UwMd6kBdUfenhDVTAx3qQF1R94TY6Ps3LSGy0amfoRHFy9i56Y1zpvrdc65Sw24Y/21bH/w6QumZmY78p4YG+W+/ceA5uaunTMfXIa71ACXIqrfDHepIYMw8la5nHOXpAIZ7lrQhvXCX9JMDHdJKpDhLkkFMtwlqUCGuyQVyHCX5sgNpjXIDHdpDtxgWoPOcJfmwA2mNegMd2kOSrmeuspluGuoNfUmJK+nrkFnuEtz4PXUNei8cJg0B17VUYPOcJfmyKs6apAZ7tIC5i+lcjnnLkkFMtwlqUBOy0hDzGkVdWK4a8E6f22YM2fPse6uRxtZ7WI4q19mnJaJiBUR8b2IOBwRhyLi9qp9SUQ8EhFHqtsrWr5me0QcjYjnImJ9P/8BGm5NvQnJa8OodN3Mub8GfDoz3wncANwWEdcB24B9mbka2Fc9pvrcZuB6YANwd0SMtH1mqSFeG0almzHcM/NUZv6ouv9r4DAwCmwEdleH7QYmqvsbgfsz85XMfB44CqztdeFSHV4bRqWb1WqZiFgJjAGPA1dn5imY+gUAXFUdNgq80PJlx6u2i59ra0QciIgDk5OTs69cqsFrw6h0XYd7RLwV+Abwqcz81XSHtmnLNzRk7srM8cwcX7p0abdlaMA0NWdel9eGUem6CveIuJSpYP9aZj5YNb8YEcuqzy8DTlftx4EVLV++HDjZm3Kl3pgYG2XnpjVcNjL1IzC6eBE7N63x2jAqRjerZQL4EnA4Mz/X8qm9wJbq/hbg4Zb2zRFxeUSsAlYD+3tXstQbE2OjjF2zmPetWsIPtt1osKso3Yzc1wF/DtwYEU9WHzcDdwEfiogjwIeqx2TmIWAP8CzwbeC2zDzb/qmluXMPU6mzGd/ElJn/Sft5dICbOnzNDmBHjbqkaXVapw44Apfw2jIaUq5Tl6ZnuGsouU5dmp7hrqHkOnVpeoa7htKgrFP/+ife78W/NJC8KuQCd/4NSMMWUO5hKk3PcFdj6l5y1z1Mpc6cllEjvOSu1F+GuxrhUkapvwx3zVmdd4i6lFHqL8Ndc1J3WsWljFJ/Ge6ak7rTKoOylFEqlatlhlxTSxnrTqu4lFHqL8Ndc/L2xYs40SbIZzOt4lJGqX+cltGcOK0iDTZH7pqTUqZV/ItBpTLcGzasb/8Hp1WkQea0zALmTkZSuQz3Bcq3/0tlM9wXKN/+L5XNcF+gfPu/VDZPqC5QvVinPgg8kSu158h9iNU5Ieo6dalsjtxramopY6cTokBXa81LWacuqT1H7g2qM/LuxQnRibFRxq5ZzPtWLeEH22402KWCGO4NqbsU0ROikqbjtExDpht5dzOC9oSopOks+JH7Lfc89vq8+XyqO/L2hKik6Thyb0jdkfegnBB15C0NpgU/cm9KL0benhCV1InhXkOd1S4TY6Ps3LSGy0am/gtGFy9i56Y1BrSknhj6aZlhXWd+/jgvmSupHxb0yL3pdeaS1C8LNtxdZy6pZEM/LTNXrjOf4nSQVKahHrnXmVZxnbmkkg1tuNedVuk0wp7NOnNXu0gaVEMb7nVPaLrOXFLJ+jbnHhEbgC8AI8C9mXlXL5+/7rTKoLzDsy7nzCW105dwj4gR4J+BDwHHgR9GxN7MfLZXr9GLE5qDsM7ccJbUD/2allkLHM3Mn2bmGeB+YGMvX8ATmpLUWb/CfRR4oeXx8aqtZzyhKUmd9WvOPdq05QUHRGwFtgJcc801c3qRQZhWkaRB1K9wPw6saHm8HDjZekBm7gJ2AYyPj18Q/MPEXyqSBlG/wv2HwOqIWAWcADYDf9an16rFcJZUor6Ee2a+FhGfBL7D1FLIL2fmoX68luEsSW/Ut3Xumfkt4Fv9en5JUmdD+w5VSVJnhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQWKzOYv6xIRk8DPazzFlcAvelROL1nX7FjX7FjX7JRY129n5tJ2nxiIcK8rIg5k5njTdVzMumbHumbHumZnodXltIwkFchwl6QClRLuu5ouoAPrmh3rmh3rmp0FVVcRc+6SpAuVMnKXJLUw3CWpQEMV7hHx5Yg4HRHPtLQtiYhHIuJIdXvFgNT1dxFxIiKerD5ubqCuFRHxvYg4HBGHIuL2qr3RPpumrkb7LCLeFBH7I+LHVV1/X7U33V+d6mr8e6yqYyQiDkbEN6vHjf9Mdqir8f6KiJ9FxNPV6x+o2vrSX0MV7sBXgA0XtW0D9mXmamBf9Xi+fYU31gXwT5n5nuqjiV2pXgM+nZnvBG4AbouI62i+zzrVBc322SvAjZn5buA9wIaIuIHm+6tTXdD89xjA7cDhlsdN99d5F9cFg9Fff1S9/vm17X3pr6EK98z8D+CXFzVvBHZX93cDE/NaFB3ralxmnsrMH1X3f83UN/ooDffZNHU1Kqf8b/Xw0uojab6/OtXVuIhYDnwYuLelufGfyQ51Daq+9NdQhXsHV2fmKZgKDeCqhutp9cmIeKqatmnkT9PzImIlMAY8zgD12UV1QcN9Vv0p/yRwGngkMweivzrUBc1/j30euBM419LWeH91qAua768EvhsRT0TE1qqtL/1VQrgPqi8Cv8PUn9GngH9sqpCIeCvwDeBTmfmrpuq4WJu6Gu+zzDybme8BlgNrI+Jd811DOx3qarS/IuIjwOnMfGI+X3cm09TV+PcXsC4z3wv8CVPTkR/o1wuVEO4vRsQygOr2dMP1AJCZL1Y/kOeAfwHWNlFHRFzKVIB+LTMfrJob77N2dQ1Kn1W1vAR8n6lzKY33V7u6BqC/1gEfjYifAfcDN0bEV2m+v9rWNQD9RWaerG5PA/9e1dCX/ioh3PcCW6r7W4CHG6zldef/syp/CjzT6dg+1hDAl4DDmfm5lk812med6mq6zyJiaUQsru4vAj4I/ITm+6ttXU33V2Zuz8zlmbkS2Aw8mpkfp+H+6lRX0/0VEW+JiN88fx/446qG/vRXZg7NB3AfU39OvQocB24FfoupM8xHqtslA1LXvwFPA09V/3nLGqjrD5ia43sKeLL6uLnpPpumrkb7DPg94GD1+s8Af1u1N91fnepq/HuspcY/BL45CP01TV1Nf3+9A/hx9XEI+Jt+9peXH5CkApUwLSNJuojhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgr0/5yBfCRpydWtAAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plt.errorbar(x, y, dy, marker='o', ls='')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The model for this data is a two degree polynomial, below is a function that defines this. \n",
+    "The `Relationship` object is also created."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def two_degree(x, a, b, c):\n",
+    "    return c * x ** 2 + b * x + a"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from uravu.relationship import Relationship"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "modeller = Relationship(two_degree, x, y, dy)\n",
+    "modeller.max_likelihood()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The max likelihood (which makes no consideration of the prior) is found, "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[ 1.9562903  -1.70809576  0.30668596]\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(modeller.variable_medians)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The default prior probabilities for these variables with `uravu` are uniform in the range $[x - 10x, x + 10x)$, where $x$ is the current value of the variable.\n",
+    "\n",
+    "However, if you wanted the prior probability to be a normal distribution, centred on the current value of the varible with a width of 1, it would be necessary to create a custom prior function. \n",
+    "This function is shown below. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from scipy.stats import norm\n",
+    "\n",
+    "def custom_prior():\n",
+    "    priors = []\n",
+    "    for var in modeller.variable_medians:\n",
+    "        priors.append(norm(loc=var, scale=1))\n",
+    "    return priors"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Note that the function returns a list of 'frozen' `scipy` RV objects that describe the shape of the priors. \n",
+    "\n",
+    "To make use of these priors, they must be passed to the `mcmc` or `nested_sampling` functions as the `prior_function` keyword argument. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "100%|██████████| 1000/1000 [01:18<00:00, 12.73it/s]\n"
+     ]
+    }
+   ],
+   "source": [
+    "modeller.mcmc(prior_function=custom_prior)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "2422it [00:21, 110.52it/s, +500 | bound: 3 | nc: 1 | ncall: 18962 | eff(%): 15.410 | loglstar:   -inf < -88.203 <    inf | logz: -92.165 +/-  0.100 | dlogz:  0.001 >  0.509]\n"
+     ]
+    }
+   ],
+   "source": [
+    "modeller.nested_sampling(prior_function=custom_prior)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "-92.16+/-0.10\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(modeller.ln_evidence)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Any `scipy` [statistical function](https://docs.scipy.org/doc/scipy/reference/stats.html) that has a `logpdf` class method may be used in the definition of priors. "
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "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.7.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

docs/source/tutorials.rst

@@ -9,6 +9,7 @@ If you want to go even deeper, then check out `Information Theory, Inference, an
 1. `Maximum likelihood`_
 2. `Markov chain Monte Carlo`_
 3. `Nested Sampling`_
+4. `Custom priors`_
 4. `Unaccounted uncertainty`_
 
 
@@ -19,6 +20,7 @@ If you want to go even deeper, then check out `Information Theory, Inference, an
    max_likelihood
    mcmc
    nested_sampling
+   custom_priors
    unaccounted_uncertainty
 
 
@@ -27,4 +29,5 @@ If you want to go even deeper, then check out `Information Theory, Inference, an
 .. _Maximum likelihood: ./max_likelihood.html
 .. _Markov chain Monte Carlo: ./mcmc.html
 .. _Nested Sampling: ./nested_sampling.html
+.. _Custom priors: ./custom_priors.html
 .. _Unaccounted uncertainty: ./unaccounted_uncertainty.html