r/rethinking/probability.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Probability Theory\n",
    "Why we need to learn probability theory. What we usually want is a model that allows us to make predictions. We may have a deterministic model of something (a law of physics), it can be a linear model or a non-linear model. In case it is linear, then we need probability to compute the error of the model from the observables. This is the easiest case. In case it is non-linear we need probability to estimate the average state of the model from the many that are allowed by the smallest differences in the initial values of the observables. If we do not have any deterministic model, like in the social sciences but also in physics when we don't know what happens inside a system, then the only thing we can do is to make observations to build a statistical model, a distribution of the outcomes given certain conditions. We can observe that some distributions occur in many cases, like the Gaussian, Poisson etc.\n",
    "\n",
    "This is the main page of my of scripts and notebooks about probability and bayesian inference. I write my notes while I read and study the material and code from several books. My main reference textbooks are\n",
    "\n",
    " - [McElreath - Statistical Rethinking, 2nd edition](https://github.com/rmcelreath/statrethinking_winter2019)\n",
    " - [Gelman et al. - Baysian Data Analysis, 3rd edition](http://www.stat.columbia.edu/~gelman/book/)\n",
    " - [Kruschke - Doing Bayesian Data Analysis, 2nd Edition](https://www.amazon.com/Doing-Bayesian-Data-Analysis-Tutorial/dp/0124058884/)\n",
    " - [DeGroot et al. - Probability and Statistics, 4th edition](https://www.amazon.com/Probability-Statistics-4th-Morris-DeGroot/dp/0321500466)\n",
    " - [Jaynes - Probability Theory: The Logic of Science](https://www.amazon.com/Probability-Theory-Science-T-Jaynes/dp/0521592712)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definitions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sample space\n",
    "A sample space $\\Omega$ is the set of all possible observations of a process or experiment. An event is a subset of the sample space. The probability of an event is a function that maps the event to a real number in the interval $[0,1]$. The probability function follows the Kolmogorov's axioms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Joint probability\n",
    "Given a vector variable y in the sample space and a parameter vector $\\theta$ in the parameter space, $P(y = y_k, \\theta = \\theta_j)$, or simply $P(y_k, \\theta_j)$, represents the joint probability of a sample in which $y = y_k$ and $\\theta = \\theta_j$. The variable $y_k$ may depend on some other variables.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the components of the sample vector x are independent variables, e.g. uncorrelated events, then  \n",
    "$P(x_k, \\theta_j) = \\prod_{i} P(x_{k,i}, \\theta_j)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conditional probability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$P(x_k,\\theta_j) = P(\\theta_j|x_k)P(x_k)$  \n",
    "$P(x_k,\\theta_j) = P(x_k|\\theta_j)P(\\theta_j)$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Marginal probability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$P(\\theta_j) = \\sum_{k} P(\\theta_j|x_k)P(x_k)$  \n",
    "$P(x_k) = \\sum_{j} P(x_k|\\theta_j)P(\\theta_j)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayes rule\n",
    "In Bayesian statistics data are the observed variables and paramenters are the unobserved variables.  \n",
    "\n",
    "$$P(x_{k}|\\theta_{j})P(\\theta_{j}) = P(\\theta_{j}| x_{k})P(x_{k})$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$P(\\theta_{j}| x_{k}) = \\frac{P(x_{k}|\\theta_{j})P(\\theta_{j})}{P(x_{k})}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$P(\\theta_{j}| x_{k}) = \\frac{P(x_{k}|\\theta_{j})P(\\theta_{j})}{\\sum_{i} P(x_k|\\theta_i)P(\\theta_i)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$P(\\theta_j)$ is called the prior probability and represents the probability of $\\theta$ before we see any data.  \n",
    "$P(x_k|\\theta_j)$ is called the likelihood and represents the probability of any mathematically plausible value of $\\theta$.    \n",
    "$P(x_k)$ is called the marginal probability  \n",
    "$P(\\theta_j|x_k)$ is called the posterior probability  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|           |$\\theta_1$|$\\theta_2$|$\\theta_3$|$\\theta_4$|\n",
    "|-----------|-----------|----------|---------|----------|\n",
    "|$x_1$      |$P(x_1,\\theta_1)$|$P(x_1,\\theta_2)$|$P(x_1,\\theta_3)$|$P(x_1,\\theta_4)$|\n",
    "|$x_2$      |$P(x_2,\\theta_1)$|$P(x_2,\\theta_2)$|$P(x_2,\\theta_3)$|$P(x_2,\\theta_4)$|\n",
    "|$x_3$      |$P(x_3,\\theta_1)$|$P(x_3,\\theta_2)$|$P(x_3,\\theta_3)$|$P(x_3,\\theta_4)$|\n",
    "|$x_4$      |$P(x_4,\\theta_1)$|$P(x_4,\\theta_2)$|$P(x_4,\\theta_3)$|$P(x_4,\\theta_4)$|"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Distributions\n",
    "Distributions are used to represent priors and likelihoods when building a model. A first differentiation is between discrete and continuous distributions. Discrete distributions are used to represent proceses whose data are integers numbers. Continuous distributions are used to represent processes whose data are real numbers. Another characteristic that guides us to choose a distribution is the scale type, the type of values we are going to deal with. Some continuous distributions are defined in an closed interval, e.g. $[0, 1]$ or infinite, e.g. $[0, + \\infty)$. The discrete Bernoulli distribution is defined in the set {0, 1}. There exist many standard distributions used for the analysis of experiments or observations that are shown in the next sections. When it is used to represent the likelihood, a distribution is selected also depending on the sampling method or stop rule, that is the rule that applies when we collect a sample. The two main sampling methods, or stop rules, are set or by the number of events, that is the sample size, or by the time used to gather the data of the sample. For example the binomial distribution is used when our sample has a specific size. The same applies for the normal distribution, the bernoulli distribution and others. The Poisson disribution and the exponential distribution are used when the stop rule is the time, or the space, between events. The likelihood represents all the possible states in which our system can be found. The data in our sample tell us which of all those possible states are more or less likely taking into account what we already know about the system before we start to collect the data and that is represented by the prior distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Binomial distribution\n",
    "A sample of the discrete binomial distribution is characterized by  \n",
    "*n* the number of events (e.g. tosses of a coin)  \n",
    "*k* the number of successes (e.g. heads)  \n",
    "*$\\theta$* the probability of success in a toss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$Binomial(k, n, \\theta) = \\binom{n}{k}\\theta^k(1 - \\theta)^{(n - k)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$P(x_k|\\theta_j) =\\binom{n}{k}{\\theta_j}^k(1 - \\theta_j)^{(n - k)}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAM1BMVEUAAABNTU1oaGh8fHyH\nzuuMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD///8LNvnPAAAACXBIWXMAABJ0\nAAASdAHeZh94AAAY/0lEQVR4nO3d7VbbTLKAURmOYwgf5v6v9mADCTOTBGFXd1e19v4x46zF\nm26UemIsy8ryAlxtGb0BmIGQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQ\nIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQ\nIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQ\nIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQ\nIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQ\nIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIICQIECHkBYo5oIpjw9nwBIQSUgQQEgQ\nQEgQQEgQQEgQQEgQQEgQQEgQQEgQQEgQQEgQQEgQQEgQQEgQQEgQQEgQQEgQQEgQQEgQQEgQ\nQEgQQEhT+7/RG9gMIU1NSL0IaWpC6kVIUxNSL0KampB6EdLUhNSLkKYmpF6ENDUh9SKkqQmp\nFyFNTUi9CGlqQupFSFMTUi9CmpqQehHS1ITUi5CmJqRehDQ1IfUipKkJqRchTU1IvQhpakLq\nRUhTE1IvQpqakHoR0tSE1IuQpiakXoQ0NSH1IqSpCakXIU1NSL0IaWpC6kVIUxNSL0KampB6\nEdLUhNSLkKYmpF6ENDUh9SKkqQmpFyFNTUi9CGlqQupFSFMTUi9CmpqQehHS1ITUi5CmJqRe\nhDQ1IfUipKkJqRchTU1IvQhpakLqRUhTE1IvQpqakHoR0tSE1IuQpiakXoQ0NSH1IqSpCakX\nIU1NSL0IaWpC6kVIUxNSL0Iq7OtMvvwKpQURUmFCykNIhQkpDyEVJqQ8hFSYkPIQUmFCykNI\nhQkpDyEVJqQ8hFSYkPIQUmFCykNIhQkpDyEVJqQ8hFSYkPIQUmFCykNIhQkpj/4h3d8sy/6h\n6RJbIaQ8Ooa0nP/D2+Xs0GSJjRFSHr1DOiyH48vL82G5b7HExggpj94h7Zbj6fFxuWmxxMYI\nKY/eIS3Lp1+EL7ExQsqjd0g/PkLatVhiY4SUR9eQ9nf3D8vP14fHw7/PNghpFSHl0TWkN+eH\nu2OLJTZGSHn0fB/p6en+fr8/n3I4/LMjIa0jpDxc2VCYkPIQUmFCykNIhQkpDyEVJqQ8hFSY\nkPIYcPr791nw8CU2Rkh5dAzpXkjBhJRH1/eRdretl9gWIeXR9TXS0xcfQwpYYlOElEffkw33\ny9Pff9u1P/fxQUh5OGtXmJDyEFJhQspDSIUJKQ8hFSakPEaF5H2kAELKQ0iFCSkPP9oVJqQ8\nhFSYkPIQUmFCyqNrSI93+/N1C/vDY6slNkVIeXQM6Xjz6Rqgf1++KqRVhJRHx5AOy+7n26V2\nzw8797ULIKQ8Ooa0+3TF6pM7rQYQUh7d/1mXP/0ibImNEVIenpEKE1IefV8jPTyfH3mNFENI\nefQ8/X376azdjXt/X09IefR9H+lwfh9pt7/zPlIEIeXhyobChJSHkAoTUh5CKkxIeQipMCHl\nIaTChJSHkAoTUh5CKkxIeQipMCHlIaTChJSHkAoTUh5CKkxIeQipMCHlIaTChJSHkAoTUh5C\nKkxIeQipMCHlIaTChJSHkAoTUh5CKkxIeQipMCHlIaTChJSHkAoTUh5CKkxIeQipMCHlIaTC\nhJSHkAoTUh5CKkxIeQipMCHlIaTChJSHkAoTUh5CKkxIeQipMCHlIaTChJSHkAoTUh5CKkxI\neQipMCHlIaTChJSHkAoTUh5CKkxIeQipMCHlIaTChJSHkAoTUh5CKkxIeQipMCHlIaTChJSH\nkAoTUh5CKkxIeQipMCHlIaTChJSHkAoTUh5CKkxIeQipMCHlIaTChJSHkAoTUh5CKkxIeQip\nMCHlIaTChJSHkAoTUh5CKkxIeQipMCHlIaTChJSHkAoTUh5CKkxIeQipMCHlIaTChJSHkAoT\nUh5CKkxIeQipMCHlIaTChJSHkAoTUh5CKkxIeQipMCHlIaTChJSHkAoTUh5CKkxIeQipMCHl\nIaTChJSHkAoTUh5CSisgk4iQpLaKkNISUiVCSktIlQgpLSFVIqS0hFSJkNISUiVCSktIlQgp\nLSFVIqS0hFSJkNISUiVCSktIlQgpLSFVIqS0hFSJkNISUiVCSktIlQgpLSFVIqS0hFSJkNIS\nUiVCSktIlQgpLSFVIqS0hFSJkNISUiVCSktIlQgpLSFVIqS0hFSJkNISUiVCSktIlQgpLSFV\nIqS0hFSJkNISUiVCSktIlQgpLSFVIqS0hFSJkNISUiVCSktIlQgpLSFVIqS0hFSJkNISUiVC\nSktIlQgpLSFVIqS0hFSJkNISUiVCSktIlQgpLSFVIqS0hFSJkNISUiVCSktIlQgpLSFV0jWk\nx7v9crI/PLZaYiJCqqRjSMeb5bfbJktMRUiVdAzpsOx+Pp0fPT/slkOLJaYipEo6hrRbnn49\nflp2LZaYipAq6RjSsvztF2FLTEVIlXhGSktIlfR9jfTwfH7kNdIaQqqk5+nv209n7W6OTZaY\niZAq6fs+0uH8PtJuf+d9pK8JqRJXNqQlpEqElJaQKrkypJu757Ct/GWJzRJSJVeGdDpt0KAl\nIb0IqZYrQzr+/NGiJSG9CKmWgNdIj3c3q1pa/lPwruYjpEpiTjY87V7LuP/iv7sX0rcIqZKQ\nkB5uV3w04pTbV1/x9yU2SEiVXB/S8e716ejm4fha0/6L//Lp3xcGXbWr+QipkmtDejydbDi8\nXY367x/XTu4/Xbf6P7/t2p/7tkJIlVz7PtLrk9H9x2Vz/76i+9IlNktIlVz7PtL+IWwrf1li\ns4RUybXvI4Vt5K9LbJaQKrn6yob3B7uwH+v+e4nNElIlQSE9f/f8gPeRviSkSq4I6eE/zrPd\nfPM3EdJXhFTJNc9In+9Td/PFR/Wa72o+Qqok6jVSLCG9CKkWH+xLS0iVXBHS6dnoexcjuPf3\ndwipko4huff39wipEvf+TktIlbjTalpCquSqH+2+d8G2e39/j5Aq6RiSZ6TvEVIl7v2dlpAq\nce/vtIRUSd/3kdz7+xuEVEnXkFruaj5CqsQlQmkJqRIhpSWkSq4O6efpZU/0nRuE9CKkWq4N\n6eNM3Fe3tLtiia0SUiVXhnR6b+j1/x52y13Ujv57ic0SUiVXhvRxtcIXVypcs8RmCamSqE/I\nOv0dTkiVXP2j3ccz0sq7en9/ic0SUiXXnmzYn18jPe5+BO3nD0tslZAq6Xj1d9tdzUdIlQgp\nLSFV4sqGtIRUiZDSElIlUSE9hl7aIKQXIdVybUgHr5FaEVIlV7+P9CH0slUhvQiplqsvEfr5\ncrs8P98ubqIfTUiVBFwidPf6bPT0xa1Tr1his4RUSUBID8u9a+0aEFIlV4a0f/3R7nm5eXkU\nUjghVXJlSA+ngM4f7gu92E5IL0Kq5drT33enX/1YYi/+FtKJkCpxZUNaQqpESGkJqRJ3EUpL\nSJW4i1BaQqrEXYTSElIl7iKUlpAqcRehtIRUibsIpSWkStxFKC0hVeLmJ2kJqRIhpSWkSlzZ\nkJaQKhFSWkKq5PpLhE7XNux/Bm3nj0tslJAqibpEKPST5kI6EVIlV4Z0/+sSofuoHf33Epsl\npEquDOnm1xuyNzH7+d8lNktIlbhEKC0hVRL2jOSi1WhCqsRrpLSEVImzdmkJqZKgj5p7Hyme\nkCpxZUNaQqrk2jutxt7P7k9LbJaQKok6/R1LSC9CquXq09/HsK38ZYnNElIlV4Z03N+G/sNI\nf1his4RUydU/2vlgXytCqkRIaQmpEqe/0xJSJUJKS0iVxFzZ8MNN9OMJqRI30U9LSJWE3UTf\n1d/RhFRJ2E30fUI2mpAq8QnZtIRUSdhN9ENfJAnpRUi1XP2vmr/fRN8H+8IJqZK4Kxsir24Q\n0ouQahFSWkKqxJUNaQmpEiGlJaRKhJSWkCoR0ihdZjzHIlsgpFFyzLiQgghplBwzLqQgQhol\nx4wLKYiQRskx40IKIqRRcsy4kIIIaZQcMy6kIEIaJceMCymIkEbJMeNCCiKkUXLMuJCCCGmU\nHDMupCBCGiXHjAspiJBGyTHjQgoipFFyzLiQgghplBwzLqQgQholx4wLKYiQRskx40IKIqRR\ncsy4kIIIaZQcMy6kIEIaJceMCymIkEbJMeNCCiKkUXLMuJCCCGmUHDMupCBCGiXHjAspiJBG\nyTHjQgoipFFyzLiQgghplBwzLqQgQholx4wLKYiQRskx40IKIqRRcsy4kIIIaZQcMy6kIEIa\nJceMCymIkEbJMeNCCiKkUXLMuJCCCGmUHDMupCBCGiXHjAspiJBGyTHjQgoipFFyzLiQgghp\nlBwzLqQgQholx4wLKYiQRskx40IKIqRRcsy4kIIIaZQcMy6kIEIaJceMCymIkEbJMeNCCiKk\nUXLMuJCCCGmUHDMupCBCGiXHjAspiJBGyTHjQgoipFFyzLiQgghplBwzLqQgQholx4wLKYiQ\nRskx40IKIqRRcsy4kIIIaZQcMy6kIEIaJceMCymIkEbJMeNCCiKkUXLMuJCCCGmUHDMupCD9\nQ7q/WZb9Q9MlSsgx40IK0jGk5fwf3i5nhyZLVJJjxoUUpHdIh+VwfHl5Piz3LZaoJMeMCylI\n75B2y/H0+LjctFiikhwzLqQgvUNalk+/CF+ikhwzLqQgvUP68RHSrsUSleSYcSEF6RrS/u7+\nYfn5+vB4+PfZBiGt+QIhJdI1pDfnh7tjiyUqyTHjQgrS832kp6f7+/3+fMrh8M+OhLTqC4SU\niCsbRskx40IKIqRRcsy4kIJ0Denxbn9+lbQ/PLZaoo4cMy6kIB1DOt4sv902WaKSHDMupCAd\nQzosu59P50fPDzunv3PMuJCCdAxptzz9evzkDdkcMy6kIN2v/v7TL8KWqCTHjAspiGekUXLM\nuJCC9H2N9PB8fuQ10kuWGRdSkJ6nv28/nbW7+Z9LG5bPLl2ikBwzLqQgfd9HOpzfR9rt77yP\nlGTGhRTElQ2j5JhxIQUR0ig5ZlxIQYQ0So4ZF1KQUSF5HynHjAspiJBGyTHjQgriR7tRcsy4\nkIIIaZQcMy6kIEIaJceMCymID/aNkmPGhRTEB/tGyTHjQgrig32j5JhxIQXxMYpRcsy4kIL4\nYN8oOWZcSEE8I42SY8aFFMQH+0bJMeNCCpLng30hSxSSY8aFFMQH+0bJMeNCCuLKhlFyzLiQ\ngghplBwzLqQgQholx4wLKYiQRskx40IKIqRRcsy4kIIIaZQcMy6kIEIaJceMCymIkEbJMeNC\nCiKkUXLMuJCCCGmUHDMupCBCGiXHjAspiJBGyTHjQgoipFFyzLiQgghplBwzLqQgQholx4wL\nKYiQRskx40IKIqRRcsy4kIIIaZQcMy6kIEIaJceMCymIkEbJMeNCCiKkUXLMuJCCCGmUHDMu\npCBCGiXHjAspiJBGyTHjQgoipDaqzHiSReoTUhvzzLiQVhFSG/PMuJBWEVIb88y4kFYRUhvz\nzLiQVhFSG/PMuJBWEVIb88y4kFYRUhvzzLiQVhFSG/PMuJBWEVIb88y4kFYRUhvzzLiQVhFS\nG/PMuJBWEVIb88y4kFYRUhvzzLiQVhFSG/PMuJBWEVIb88y4kFYRUhvzzLiQVhFSG/PMuJBW\nEVIb88y4kFYRUhvzzLiQVhFSG/PMuJBWEVIb88y4kFYRUhvzzLiQVhFSG/PMuJBWEVIb88y4\nkFYRUhvzzLiQVhFSG/PMuJBWEVIb88y4kFYRUhvzzLiQVhFSG/PMuJBWEVIb88y4kFYRUhvz\nzLiQVhFSG/PMuJBWEVIb88y4kFYRUhvzzLiQVhFSG/PMuJBWEVIb88y4kFYRUhvzzLiQVhFS\nG/PMuJBWEVIb88y4kFYRUhvzzLiQVhFSG/PMuJBWEVIb88y4kFYRUhvzzLiQVhFSG/PMuJBW\nEVIb88y4kFYRUhvzzLiQVhFSG/PMuJBWEVIb88y4kFYRUhvzzLiQVhFSG/PMuJBWEVIb88y4\nkFYRUhvzzLiQVhFSG/PMuJBWEVIb88y4kFYRUhvzzLiQVhFSG/PMuJBWEVIb88y4kFYRUhvz\nzLiQVhFSG/PMuJBWEVIb88y4kFYRUhvzzLiQVhFSG/PMuJBWEVIb88y4kFYRUhvzzLiQVhFS\nG/PMuJBWEVIb88y4kFYRUhvzzLiQVhFSG/PMuJBWEVIb88y4kFYRUhvzzLiQVhFSG/PMuJBW\nEdIlkozflhbJTkiXKDN+8yySnZAuUWb85lkkOyFdosz4zbNIdkK6RJnxm2eR7IR0iTLjN88i\n2QnpEmXGb55FshPSJcqM3zyLZCekS5QZv3kWyU5IlygzfvMskp2QLlFm/OZZJDshXaLM+M2z\nSHZCukSZ8ZtnkeyEdIky4zfPItkJ6RJlxm+eRbIT0iXKjN88i2QnpEuUGb95FslOSJcoM37z\nLJKdkC5RZvzmWSQ7IV2izPjNs0h2QrpEmfGbZ5HshHSJMuM3zyLZCekSZcZvnkWyE9Ilyozf\nPItkJ6RLlBm/eRbJTkiXKDN+8yySnZAuUWb85lkkOyFdosz4zbNIdkK6RJnxm2eR7IR0iTLj\nN88i2XUN6fFuv5zsD4+tluijzPjNs0h2HUM63iy/3TZZopcy4zfPItl1DOmw7H4+nR89P+yW\nQ4sleikzfvMskl3HkHbL06/HT8uuxRJB5hm/DS0yWMeQluVvvwhbIkiOybDI975iLM9If5Bj\nMizyva8Yq+9rpIfn86Psr5FyTIZFvvcVY/U8/X376azdzbHJEjFyTIZFvvcVY/V9H+lwfh9p\nt7/L/T5SjsmwyPe+YixXNvxBjsmwyPe+Yiwh/UGOybDI975iLJcI/UGOybDI975iLJcI/UGO\nybDI975iLJcI/UGOybDI975iLG/I/kGOybDI975irA1eIlRlMiwSvEhTG3xGmmcyLPKt36Kp\nDV4iNM9kWORbv0VTeS4RWj7722/R56ftKpNhkTyLuEQIAriyAQIICQIICQKMCin1R83hu4QE\nAfxoBwGEBAGEBAF8sA8C+GAfBPDBPgiwwY9RQLwNfrAP4nlGggAb/GAfxMvzwb6QJWAMH+yD\nAK5sgABCggBCggBCggBJQ4JiLpjy+HCqKXMIqmx0k/us8k03VOYQVNnoJvdZ5ZtuqMwhqLLR\nTe6zyjfdUJlDUGWjm9xnlW+6oTKHoMpGN7nPKt90Q2UOQZWNbnKfVb7phsocgiob3eQ+q3zT\nDZU5BFU2usl9VvmmGypzCKpsdJP7rPJNN1TmEFTZ6Cb3WeWbbqjMIaiy0U3us8o33VCZQ1Bl\no5vcZ5VvGlITEgQQEgQQEgQQEgQQEgQQEgQQEgQQEgQQEgQQEgQQEgQQEgQQEgQQEgQQEgTY\neEgX3zO9r/uPHR52y+7wz39odKiPfeY+rPc3vw5i3PHM+s328ZT7T/zD08cO3/7h3puxu/m7\nj33mPqyH8952p3wCj2fSb7aTp2U/egsrPO3eZ/Jx2T2dfvXFP9o7yq99pj6sT8uP4+m580fs\n8dx2SPfL3egtfO1+uX0f0MPy8Pq/P5Nu+vc+Ux/W/dseT1uNPJ5bD+l+9Ba+thxe3gd0vzy/\n5P37/vc+KxzW01Yjj+e2Q9ovDz9eX22O3sa/Pb18DOh//l82v/dZ4LAel9vY45nzz6SX/duL\n4tvR+/hKiZBePoWU/rDen36qE1KUZfn5+pfTIf1PIsVCyn9Yn3enH+eEFOuY94Tyu2IhvUl8\nWI+787OlkILlncx37xvclQop8T5v3xKPPJ5pv9ee8v6Jv/uPs3bPSc/avVQJ6fnm9vn8IPJ4\nJv1eO9ktp/e3E0/mu/eRvDu/7/GwpD0f9uuZM/Nhffh1FiTyeG47pMPpGB7f3pfLrMaVDb/2\nmfqwPv8+m+jKhijH3fk8bdq/4T98/JB0k/y08vs+Ux/WH8vvKwEDj+e2Q3r9a3O33OQ9S/vh\nI6Tj+WrlsXv5l8/7zHpYl08hBR7PjYcEMYQEAYQEAYQEAYQEAYQEAYQEAYQEAYQEAYQEAYQE\nAYQEAYQEAYQEAYQEAYQEAYQEAYQEAYQEAYQEAYQEAYQEAYQEAYQEAYQEAYQEAYQEAYQEAYQE\nAYQEAYQEAYQEAYQEAYRUU9Z/MXyz/HnUJKRk/HnUJKRk/HnUJKRk/HnU9BbSYbkbvRHeCKmm\nc0iH5Xb0PngnpJpOIekoESHV9BqSjjIRUk3Lcrssj6N3wS9CqmlZlt1yM3oX/CKkml6fjp6W\n5efobfBBSDWdTjbcLbvR2+CDkGo6n/6+8TZSGkKq6RzS6w93x9Eb4Y2Qanq7suFu2Y/eCG+E\nBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGE\nBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAH+H7AVor/Y8Wv6AAAAAElF\nTkSuQmCC",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N <- 20 # number of flips\n",
    "k <- seq(from = 1, to = N, by = 1) # number of successes (e.g. heads)\n",
    "theta <- 0.5 # probability of success at every flip\n",
    "y_k <- dbinom(k, size = N, prob = theta) # predicted variable\n",
    "plot(k, y_k, type = \"h\", col = \"skyblue\", ylab = \"probability\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bernoulli distribution\n",
    "For N = 1 and $k \\in {0, 1}$ the binomial distribution becomes\n",
    "\n",
    "$$Bernoulli(\\theta, k) = \\theta^k (1 - \\theta)^{(1 - k)}$$  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Beta distribution\n",
    "The beta continuous distribution is defined over the interval [0, 1]. Its shape depends on two positive parameters, a and b. It can be used as conjugate prior distribution of a binomial distribution likelihood."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$Beta(\\theta|a,b) = \\frac{\\theta^{(a - 1)}(1 - \\theta)^{(b - 1)}}{B(a,b)}$$\n",
    "\n",
    "where\n",
    "\n",
    "$$B(a,b) = \\int^1_0 \\theta^{(a - 1)}(1 - \\theta)^{(b - 1)}d\\theta$$  \n",
    "\n",
    "is the beta function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other parameters of the beta distribution can be computed from the shape parameters $a$ and $b$.  \n",
    "\n",
    "$$\\mu = \\frac{a}{a + b}$$  \n",
    "$$\\omega = \\frac{a - 1}{a + b - 2}$$  \n",
    "$$k = a + b$$  \n",
    "\n",
    "where $\\mu$ is the mean of the beta distribution, $\\omega$ the mode and $k$ the concentration. The larger is k the narrower, or more concentrated, is the distribution. We can invert the relationships to get the shape parameters $a$ and $b$.\n",
    "\n",
    "$$a = \\omega(k - 2) + 1 $$  \n",
    "$$b = (1 - \\omega)(k - 2) + 1$$    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAM1BMVEUAAABNTU1oaGh8fHyH\nzuuMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD///8LNvnPAAAACXBIWXMAABJ0\nAAASdAHeZh94AAAfm0lEQVR4nO3d61oqydZt4UxYiIiA93+1S1BneYA8RPQYPTKyvT++7apn\nzwoZMVqpiM7uDUC2zv0OAC0gJECAkAABQgIECAkQICRAgJAAAUICBAgJECAkQICQAAFCAgQI\nCRAgJECAkAABQgIECAkQICRAgJAAAUICBAgJECAkQICQAAFCAgQICRAgJECAkAABQgIECAkQ\nICRAgJAAAUICBAgJECAkQICQAAFCAgQICRAgJECAkAABQgIECAkQICRAgJAAAUICBAgJECAk\nQICQAAFCAgQICRAgJECAkAABQgIECAkQICRAgJAAAUICBAgJECAkQICQAAFCAgQICRAgJECA\nkAABQgIECAkQICRAgJAAAUICBAgJECAkQICQAAFCAgQICRAgJECAkAABQgIECAkQICRAgJAA\nAUICBAgJECAkQICQAAFCAgQICRAgJECAkAABQgIECAkQICRAgJAAAUICBAgJECAkQICQAAFC\nAgQICRDIDOmw6brdceQIYGHml5Aa0sdZ249j90WOAEyiQ9p3+8vb23nfHUocAZhEh9R3l+vb\nl25T4gjAJDqkr88lhz+nJCQsTHRIT18h9SWOAExCQ9o9H47dy/ubl/3wsw2EhIUJDenf84Rd\n119KHAGYBIb0djodDrvd7SmH/WBHhISliQypqiMAJUICBAgJEHCFxPeR0JR6Qsp8BSDgxKd2\ngAAhAQKEBAhEhnR56rrt54/08WQDmhIY0qW/PY+w+/iXEBJaEhjS7Yf5Lod+e/uXEBJaEhhS\n//EHz/3mTEhoTPjvbHj/oLTdEhIaExjSpvt6yfdmS0hoS2BIh+7p861ztyUkzPO/b9zvyx2R\nT3/v/9VzHHkVECHhh//94X6Pfgv9huxp9/XW+YmQMNXfjOpLiVc2oHL3M6otJUJC1R5nVFdK\nhISKDWdUU0mEhHqNdlRPSYSEWk3IqJ6SCAmVmtZRLSUREuo0taNKSiIk1Gh6RpWkREioESFJ\nENLKzeuohpIICdWZm1ENJRESqpMQkr0kQkJlUjIiJNsRqFVaSO6UCAlVSc2IkExHoE7pIXlL\nIiTUJKMjb0mEhIpkdURIjiNQobyOrCUREuqRG5KxJEJCLbIzIiTDEaiOICRfSYSESig6IqTw\nI1AZSUe+kggJdRCF5CqJkFADVUaEFH0EqkJIRRDSyug6cpVESKiAMiRPSYQEP2lHhBR6BOqh\n7chTEiHBjpBKIaQ1UXdkKYmQYKbviJACj0AlCnTkKImQYFWkI0KKOwJ1KBNSfEmEBKdCHRFS\n2BGoQamO4ksiJBgRUlmEtA7lOiKkqCNQgYIhRZdESLAp2REhBR0Bv6IhBZdESHAp2xEhxRwB\nO0IqjpBWoHBHwSUREjyKdxRbEiHBg5AICfkCOiIkQmpeREeEREjNCwkpsiRCgkFMR5ElERIM\nCCntj1R4BIyiOiKk8kfAiJAS/0iFR8AnrKPAkggJ4Qgp9Y9UeARsAjsiJDQrsqO4kggJwQgp\n/Y9UeARMYjsKK4mQECq6I0JCkwgp549UeAQ8CCnnj1R4BCzCO4oqiZAQyNARIaE9jpBiSiIk\nxLF0REhoDSFl/pEKj4CBJ6SQkggJYUwdERKa4uoopCRCQhRCyv4jFR6BcISU/UcqPALRfB0R\nEtph7CiiJEJCDEIS/JEKj0Asa0eEhFZ4QypfEiEhgrkjQkIbCEnyRyo8ApHcHRESmuDuqHxJ\nhIQA7owICS1wV/Q/QkID3BFdlX6MhITi3BHdFH6MhITi3A3dFH6MhITS3Al9KvsgCQmFuQP6\nUvZRhob0+rzrrnb711JHoDrugL6UfZSBIV023X+2RY5AhdwBfSn7KAND2nf9y+n21vnYd/sS\nR6A+7n7+KfswA0Pqu9O/t09dX+IIVMedzzdFH2dgSF336H/IjkB13PV8U/Rx8hEJRbnr+a7k\n44z9Gul4vr3F10ir4W7nh5IPNPLp7+23Z+02lyJHoC7udH4q+Uhjv4+0v30fqd89832kdXCn\n81PJR8orG1CQO52fSj5SQkI57nJ+K/hQCQnluMP5reBDdYXE95FWwN3NX+Ueaz0hdd8pjoCb\nO5u/yj1WPrVDKe5q7ij3YAkJpbiruaPcgyUklOKu5o5yDzYypMtT122Pn/8SnmxonTuau4o9\n2sgf7Os/fjz2419CSI1zJ3NfsYcb+qLVw3tNh/72w7GE1Dp3MvcVe7ihP0Zx+3/O/eZMSO1z\nJ3NfsYdr+MG+y3ZLSM1zF/NAsccbGNKm+/rRic2WkBrnDuaRYg84MKRD9/T51rnbElLb3ME8\nVOoBRz79vf9Xz3HkVUCEtHTuXh4q9YBDvyF72n29dX4ipJa5cxlQ6BHzygYU4K5lQKFHTEgo\nwF3LgEKPmJCg545lSKGHTEiQc7cyqNBjJiTIuVsZVuYxExLk3KkMK/OYCQlq7lJGlHnQhAQx\ndyijijxqQoKYu5NRRR41IUHM3cmoIo+akKDlzmRckYdNSNByZzKuyMMmJEi5K5mgyOMmJEi5\nK5mixOMmJCi5G5mkxAMnJCi5G5mmwAMnJAi5C5mowCMnJAi5C5mowCMnJAi5C5mowCMnJOi4\nA5mqwEMnJOi4A5mqwEMnJOi4A5lM/9AJCTLuPKbTP3ZCgoq7jhn0D56QoOKuYw75gyckqLjj\nmEP+4AkJIu42ZpE/ekKCiLuNWeSPnpCg4U5jHvnDJyRouNOYR/7wCQkS7jLmUj9+QoKEO4y5\n1I+fkKDg7mI29QAICQruLmZTD4CQoODuYjb1AAgJAu4s5lNPgJAg4M4igXgChAQBdxUJxBMg\nJORzR5FCPAJCQjZ3E2m0MyAkZHMnkUY7A0JCNncSabQzICTkcheRSDsEQkIudxGJtEMgJGRy\nB5FKOwVCQiZ3EKm0UyAkZHIHkUw6BUJCHncO6aRjICRkcdeQQToHQkIWdw05lHMgJGRxx5BD\nOQdCQg53C1mUgyAk5HC3kEU5CEJCBncKeZSTICRkcKeQRzkJQkIGdwqZhJMgJKRzh5BLOApC\nQjp3CLmEoyAkpHOHkEs4CkJCMncH+XSzICSkclcgoBsGISGVuwIB3TAICancFQjohkFISOSO\nQEE3DUJCIncECrppEBLSuBvQkI2DkJDGnYCGbByEhDTuBDRk4yAkJHEXICKbByEhhTsAFdlA\nCAkp3AGoyAZCSEjhDkBFNhBCQgL3/uuoJkJISOBefx3VRAgJCdzrr6OaCCFhPvf2C6lGQkiY\nz739SqKREBJmc+++lGgmhITZ3LsvJZoJIWE29+5LiWZCSJjLvfpaoqEQEmZyb76YaCqEhJnc\nmy8mmgohYSb35qtppkJImMe993KasRAS5nHvvZxmLISEWdxrr6eZCyFhFvfaFyCZCyFhFvfW\nFyCZCyFhDvfSlyAZTHxIh03X7Y5Fj0Ap7p0vQjKZwJC62x/cdjf7IkegMPfOFyGZTHRI+25/\neXs777tDiSNQmHvni5BMJjqkvrtc3750mxJHoCz3ypchGU10SF337X/Ij0BZ7pUvRDGa6JCe\nvkLqSxyBstwbX4hiNKEh7Z4Px+7l/c3LfvjZBkKqknvhS1HMJjSkD7c3+0uJI1CUe+FLUcwm\n8vtIp9PhsNvdnnLYD3ZESFVy73s5guHwygZM5V73cgTDISRM5V73cgTDISRM5N72ggTTcYXE\n95GWxr3sJQnGU09I3XeKIyDlXvaSBOPhUztM4172kgTjISRM4t71svLnQ0iYxL3qZeXPh5Aw\nhXvTC8sfECFhCvemF5Y/IELCFO5NLy17QIYXrU54hpuQKuPe8+KyJxQY0oGQFsu958VlTyj0\n1d/9tvQRKMO958VlTyj0a6TTyC8PEhyBEtxrXl72iGKfbDh0p9JHoAD3mpeXPSKetcMo95YH\nyJ4RIWGUe8sDZM+IkDDKveURcmdESBjj3vEQuUMiJIxwr3iM3CkREka4VzxG7pQyQ9o8n3Pf\ng7EjYOZe8Ri5U8oMqeu6Ei0RUj3cGx4kd0yZIV1enkq0REj1cG94kNwxCb5Gen3eqFsipGq4\nFzxK7pw0Tzac+vePS4N/dVj2EfBwL3iYzDlJQjp+/H2WU1/bnXIETNz7HSZzTvkhXZ7fPxxt\njpf3mnaZ78yjI2DjXu84mYPKDen1+mTD/uM13brf60hIlXBvd6DMSeV+H+n9g9Hh669oGf5b\n+FKPgJF7uyPlTSr3+0i7Y97540fAyL3ckfImlft9pLzTpxwBH/duh8obVfYrGz7f6GWf1v0+\nAj7u3Q6VNypRSGftXyBBSHVw73aovFFlhHT88du1NnnvR/Z7BT33asfKm1XOR6TN945e896P\n7PcKcu7NDpY3LNXXSFqEVAP3ZgfLGxY/2IdH3JsdLG9YGSFdPxoV+tsqCakC7sUOlzUtQsID\n7r0OlzUtPrXDfe61jpc1LkLCfe61jpc1rqxP7Sb/NS0B7xXE3GttkDMuQsJd7qV2yJkXn9rh\nHvdOW+QMjJBwj3unLXIGxtPfuMe90xY5AyMk3OFeaY+cifGpHe5wr7RHzsQICXe4V9ojZ2LZ\nIb3s3j+tU//mBkLycm+0S8bIckPafn6FJPuVdn+PQDj3QrtkjCwzpH3XXz8YHfvuOeOdGDwC\n4dz7bJMxs8yQ+u7jd0OedL/T7vcRCOfeZ5uMmal+Qpanvxvi3mebjJllf2r39RFpn/FODB6B\naO51NkofWu6TDbvb10iv/VP6uzB2BGK5l9kpfWq8+hu/uJfZKX1qhIRf3MvslD41XtmAn9y7\nbJU+NkLCT+5dtkofmyqkV+lLGwjJxr3KXulzyw1pz9dIbXGvslf63LK/j/RF+rJVQrJxr7JX\n+tyyXyL08rbtzudtxy/Rb4J7k92SByd4idDz+0ejU7dNfhdGjkAg9x7bJU9OENKxO/Bau0a4\n99gueXKZIe3eP7U7d5u3V0JqgnuP7ZInlxnS8RrQ7Yf7pC+2IyQP9xr7JY8u9+nv5+v/euq0\nL/4mJBP3Gvslj45XNuA/7jX2Sx4dIeEf9xZXIHl2/BYhfHEvcQ2Sh8dvEcIX9xJXIXV4/BYh\nfHHvcBVSh8dvEcIn9wrXIXV6/BYhfHKvcB1Sp8dvEcIH9wbXInF8/BYhfHAvcC0Sx8cvP8GN\ne3+rkTg/QsKNe3+rkTg/XtmAK/f61iNxgISEK/f61iNxgPkvEbq+tmH3knj8pCNQnnt965E4\nQNVLhKQ/aU5I0dzbW5HECWaGdPj3EqFD4jswegQCuLe3IokTzAxp8+8bspvEd2D0CARwb29F\nEifIS4RARz+kjVD2EYkXrS6Xe3XrkjZDvkYCIf2QNkOetQMh/ZA2Q9GPmvN9pAVzb25tkobI\nKxtWz7231UmaYu5vWtX+Prt7R6Aw995WJ2mKqqe/tQgpkHtvq5M0xeynvy9Jx844AmW517Y+\nSWPMDOmy20r/YqQ7R6Ao99ZWKGmO2Z/a8YN9y+be2golzZGQVs69tRVKmiNPf6+be2lrlDRI\nQlo399LWKGmQmlc2PPFL9JfJvbN1Spkkv0R/1dwrW6eUScp+iT6v/l4g98ZWKmWUsl+iz0/I\nLpB7YyuVMkp+QnbN3BtbqZRRyn6JvvSLJEIK4V7YaiXMMvtvNf/8JfqTfrDv9Xn38czEfuR1\nRYQUwr2v1UqYpe6VDaOvbrhsvv3/HA6PkCK417VeCcMMDGnf9S8fnwiej/3w36dESBHc61qv\nhGEGvrLh6xm+q5HfOkRIAdzbWrGEaQaG1N19uk96BGZwb2vFEqbJR6S1ci9rzRLGGRjS9VUQ\n59tbfI1UAfey1ixhnJGv/t5+e1piM/gj6oRUnHtXq5Ywz9Afo3jd376P1O+e+T6Sm3tX6zZ/\nnvw80jq5N7Vy8wdKSOvk3tTKzR8oIa2Se1FrN3+irpD4PpKVe1FrN3+i9YQ0/bVGyOZe1OrN\nniif2q2Re03rN3ukhLRG7jWt3+yREtIKubd0AWbPlJBWyL2lCzB7poS0Pu4lXYLZQyWk9XEv\n6RLMHmrozyNNfoabkApy7+gizJ5qYEgHQqqCe0cXYfZUIz+1O037XUM5R2CUe0WXYfZYQ79G\nOg3/OJ/iCIxxr+hCzB1r7JMNh28/bV7oCAxzL+hSzJ0rz9qtjHtBl2LuXAlpXdz7uRhzB0tI\nq+Jez+WYO1lCWhX3ei7H3MkS0qq413M55k6WkNbEvZ0LMne0hLQm7u1ckLmjJaQ1cW/ngswd\nLSGtiHs5l2TubAlpPdy7uSwzh0tI6+FezWWZOVxCWg/3ai7LzOES0mq4N3NhZk6XkNbCvZhL\nM3O8hLQW7sVcnHnjJaSVcK/l8sybLyGthHstl2fefAlpHdxbuUDzBkxI6+DeygWaN2BCWgX3\nUi7RvAkT0iq4l3KJ5k2YkFbBvZRLNG/ChLQG7p1cpHkjJqQVcK/kMs2bMSGtgHsll2nejAmp\nfe6NXKpZQyak9rkXcqlmDZmQmufex8WaNWVCap57Hxdr1pQJqXXudVyuWWMmpMa5t3HJ5syZ\nkBrnXsYlmzNnQmqcexmXbM6cCalt7l1ctDmDJqSmuVdx2eZMmpCa5l7FZZszaUJqmnsVl23O\npAmpZe5NXLg5oyaklrk3ceHmjJqQGuZexKWbM2tCapd7DxdvzrAJqV3uPVy+GcMmpHa513D5\nZgybkJrl3sIGzJg2IbXKvYQtmDFuQmqVewlbMGPchNQo9w62Yfq8CalN7g1sxPSBE1Kb3BvY\niOkDJ6Q2uTewEdMHTkhNci9gK6ZPnJCa5F7AVkyfOCG1yL1/zZg+ckJqkHv92jF95oTUIPf6\ntWP6zAmpPe7ta8j0oRNSe9zb15DpQyek5riXrymTp05IzXHvXlMmT52QWuNevbZMHjshNca9\neY2ZPHdCaox78xozee6E1Bb34rVm8uAJqSnuvWvO5MkTUlPce9ecyZMnpKa49645kydPSC1x\nr117Jo+ekFriXrv2TB49ITXEvXUNmjx7QmqHe+maNHX4hNQO9841aerwCakZ7pVr09TpE1Ir\n3BvXqKnjJ6RWuDeuUVPHT0iNcC9cq6bOn5Da4N63dk28AEJqg3vd2jXxAgipCe5ta9jEGyCk\nJri3rWETb4CQmuDetoZNvAFCaoF72Vo28QoIqQHuXWvaxDsgpAa4d61pE++AkJbPvWptm3gJ\nhLR47k1r3MRbIKTFc29a4ybeAiEtnXvRmjftGghp4dxr1r5p90BIC+des/ZNuwdCWjb3lq3A\ntIsgpEVzL9kaTLsJQlo095KtwbSbiA/psOm63bHoEavh3rF1mHQVgSF1tz+47W72RY5YG/eK\nrcOkq4gOad/tL29v5313KHHEyrg3bCUm3UV0SH13ub596TYljlgX94KtxaTLiA6p6779D/kR\n6+JesLWYdBnRIT19hdSXOGJV3Pu1GpNuIzSk3fPh2L28v3nZDz/bQEjj3Ou1HpOuIzSkD7c3\n+0uJI9bEvV7rMek6Ir+PdDodDrvd7SmH/WBHhDTOvV0rMuk+eGXDMrmXa00mXQghLZN7udZk\n0oUQ0iK5d2tVJt2IKyS+j5TDvVorM+VK6gmp+05xRLvci7U2U+6ET+0WyL1YazPlTghpedx7\ntTpTLoWQFse9Vusz5VZCQ3p93t2+AtrtX0sdsQLutVqfKbcSGNJl8+3ZhG2RI9bAvVVrNOFa\nAkPad/3L6fbW+djzotVE7p1apQn3EhhS353+vX3ixygSuXdqlSbcS/jvbLj3P2RHtM+9Uus0\n4WL4iLQo7o1aqQk3E/s10vF8e4uvkRK5F2qtJlxN5NPf22/P2m34wb4E7oVaqwlXE/t9pP3t\n+0j97pnvI6Vw79NqTbgbXtmwHO51Wq8Jl0NIy+Fep/WacDmEtBjubVqxCbdDSEvhXqY1m3A9\nhLQQ7l1at/H7IaRlcG/Syo1fECEtg3uTVm78gghpEdyLtHbjN0RIi+BepLUbvyFCWgD3GmH8\njghpAdxrhPE7IqT6ubcIE0oipOq5dwj/I6QGuFcIV6PXREiVc28QbkbviZAq594g3IzeEyHV\nzb1A+DB6UYRUNff+4NPoTRFSzdzrgy+jV0VINXOvD76MXhUhVcy9Pfhn9K4IqV7u5cF/Ri+L\nkKrl3h18M3pbhFQr9+rgh7HrIqRKuRcHP43dFyHVyb03+GXswgipSu61wW9jN0ZIVXKvDX4b\nuzFCqpF7a/DH2JURUoXcS4O/xu6MkOrj3hncMXZphFQd98rgnrFbI6TauDcGd41dGyFVxr0w\nuG/s3gipLu59wQNjF0dIVXGvCx4ZuzlCqol7W/DYyNURUk3cy4LHRq6OkOrhXhUMGbk8QqqH\ne1UwZOTyCKkW7kXBsJHrI6RKuPcEI0buj5Dq4F4TjBq+QEKqgntJMG74BgmpCu4lwbjhGySk\nCrhXBFMM3yEh+bk3BJMMXyIh2bkXBNMM3yIhubn3AxMNXyMhmbnXA1MN3yMhWbmXA9MN3yQh\nWbmXA9MN3yQhObl3AzMMXyUh+bg3A7MMXyYh2bgXA/MM3yYhubj3AjMNXychebi3AvMNXigh\nebiXAvMNXighWbh3AgkGb5SQDNwbgSSDd0pI8dwLgTSDl0pI4dz7gESDt0pIwdzbgHRD90pI\nodyrgBxDN0tIkdybgCxDV0tIgdyLgDxDd0tIcdx7gExDl0tIUdxbgGxD10tIMdw7AIGhCyak\nEO4VgMLQDRNSBPcGQGLoigmpPPf9Q2TokgmpOPf1Q2XolgmpMPflQ2fongmpKPfVQ2nopgmp\nJPfNQ2roqgmpHPe9Q23gsgmpFPelQ2/gugmpDPeVo4SBCyekEtwXjjIGrpyQ9NzXjVIGLp2Q\n1NyXjXIGrp2QtNxXjaIeXzwhKbnvGYU9vnpC0nHfMop7fPmEpOK+YwR4fP2EJOG+YMR4vAGE\nJOC+XkR5vAOElMt9twj0eA0IKYv7YhHr8SaEhvT6vOuudvvXUkdEct8qwj1ehsCQLpvuP9si\nRwRyXykcHu9DYEj7rn853d46H/tuX+KIKO77hMnjlQgMqe9O/94+dX2JI0K4LxM+j7ciMKSu\ne/Q/ZEcU575IeD3eDD4iTeW+Q1Tg8XrEfo10PN/eWtzXSO77Qy0erkjk09/bb8/abS5FjtBz\n3xyq8nBPYr+PtL99H6nfPS/i+0juS0N9Hi4Lr2y4w31bqNbDnSGkb9y3hOo93B1Ccl8NlkS5\nspItt30fyX0TWDTlyhYKqfvu0R9zzxErp6yisU/tAAFCAgQICRDgB/sAAX6wDxDgB/sAAX6M\nAhDgB/sAAT4iAQL8YB8gwA/2AQL8YB8gwCsbAAFCAgQICRAgJECAkAABQgIECAkQICRAgJAA\nAUICBCoNCViYhC3Xh7OIszmf86XnExLnc35t/7IFnc35nE9InM/5tZ1PSJzP+bX9yxZ0Nudz\nPiFxPufXdj4hcT7n1/YvW9DZnM/5hMT5nF/b+YTE+Zxf279sQWdzPuc3ExLQDEICBAgJECAk\nQICQAAFCAgQICRAgJECAkAABQgIECAkQICRAgJAAAUICBAgJECAkQCA8pH3f9fvL0D8IPv+w\n8Z7/7jXwFv6cf3rquqez7fxL8P2/X/jPaYvOjw5pe/tl/5uBfxB8/v72D/qom7z3cC993C38\nOf/offzn/uP8uJJPP/+uCdX+BYf02vWnt1PfvT78B8Hnn7qny/U/Uk+m8692KX+NiOr8/v0f\nXHbd3nT+0+3kfdT8366Hf5+2bP+CQ9p3x/f/+9I9P/wHwefvPgYQtcr3Hu5L0t/HIzr/5bbI\nl643nd/Fzv/9P5nbH2fJ9i84pF13/Rh+6nYP/0Hw+Z+iLvLO+edfVxt7/lN3ijr77vmfn9VG\nhfz2/t+NH9OW7V9wSH/+AxT8X6QHx126re38bXeOC+nP+Zvu7bm/fXrrOf/581O7oM9I3k6/\nLl+2f4R0dbh9gLec/9y9xH1ic2/+u9sX+67z3w7XZxv6Q9D5vw4nJNn5N+c+6DPLv+ffPqmw\nhnR9suEp6iPCvf+QXEV9QPp1OCHJzr+69EGf2N371Or6xLM1pOvXSOeo7z/8Of9w/dTuPeTA\nD0lNhNT/fr///IPg86+2Yd/F+nP+0+1zyriQ/jz+4P+Q/Tl/012/PLvEfSPx12OV7Z/lWbvz\n72ftzrHP2v047rzZxn038Pf5OX8hveL86Kf//5wf/fT377Nk+xcc0vPtv8DH/77/9+cfBJ//\n/nbY53V3zo8O6cH8z1FD+HP+x0eEsO9jXf2YtWz/1v7KhrAVenD+jfGVDe9fHV2uX6O8mM7f\nd9fXue2j/kN61cQrG94/J766Le/HA/r2DxznP8V+RPj7+H++FX/+s3f+n691i/yv2de0tfsX\nHdLHi30/ju5+/QPH+cGfWv19/D/fMpx/3Drn//nq67Dz336HpNq/6JCAJhESIEBIgAAhAQKE\nBAgQEiBASIAAIQEChAQIEBIgQEiAACEBAoQECBASIEBIgAAhAQKEBAgQEiBASIAAIQEChAQI\nEBIgQEiAACEBAoQECBASIEBIgAAhAQKEBAgQEiBASIAAIQEChAQIEBIgQEiAACEBAoQECBDS\nQh23Xbc9ut8LfCGkZTp0Nwf3+4FPhLRMfXd6e3vpNu73A58IaZm6jk/rqkJIy7Tvut3p5H4v\n8A8hLdRz//4lUn92vxv4REiLddxv+BqpGoS0ZB3XVwtuYpk23QvP2tWEkJbp5eP7SK/u9wOf\nCGmhbq9soKNqEBIgQEiAACEBAoQECBASIEBIgAAhAQKEBAgQEiBASIAAIQEChAQIEBIgQEiA\nACEBAoQECBASIEBIgAAhAQKEBAgQEiBASIAAIQEChAQIEBIgQEiAACEBAoQECBASIEBIgAAh\nAQL/B1gEp1PWHl3AAAAAAElFTkSuQmCC",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "om <- 0.8                      # mode\n",
    "k <- 7                         # concentration (or spread)\n",
    "a <- om * (k - 2) + 1          # can be seen as the number of heads in a coin toss\n",
    "b <- (1 - om) * (k - 2) + 1    # can be seen as the number of tails in a coin toss\n",
    "s <- seq(from = 0, to = 1, by = 0.001)\n",
    "plot(s, s^(a - 1) * (1 - s)^(b - 1) / beta(a, b), type = \"h\", col = \"skyblue\", ylab = \"probability\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gamma distribution\n",
    "The gamma continuous distribution is defined in the interval $[0, + \\infty)$ and characterized by two parameters, the scale s and the rate r."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$Gamma(k| r, s) = \\frac{r^sk^{(s - 1)}e^{-rk}}{\\Gamma(s)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where\n",
    "\n",
    "$$\\Gamma(s) = \\int_0^\\infty t^{(s - 1)}e^{-t}dt$$\n",
    "\n",
    "is the gamma function that can be used to compute the factorial of complex numbers\n",
    "\n",
    "$$\\Gamma(s) = (s - 1)!$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other parameters of the gamma distribution can be computed from the shape and the rate\n",
    "\n",
    "$$\\mu = \\frac{s}{r}$$     \n",
    "$$\\omega = \\frac{s - 1}{r}$$  \n",
    "$$\\sigma = \\frac{\\sqrt{s}}{r}$$\n",
    "\n",
    "where $\\mu$ is the mean of the gamma distribution, $\\omega$ the mode and $\\sigma$ the standard deviation. We can invert the relationships to get shape and rate\n",
    "\n",
    "$$s = 1 + \\omega r$$  \n",
    "$$r = \\frac{\\mu}{\\sigma^2}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAM1BMVEUAAABNTU1oaGh8fHyH\nzuuMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD///8LNvnPAAAACXBIWXMAABJ0\nAAASdAHeZh94AAAeFUlEQVR4nO3dbVfqSpRF4QCXg6iI///XXgHRgIS8rdorVTWfD92e0eMa\nsqtmoyFg8wlgtsb9AIASEBIgQEiAACEBAoQECBASIEBIgAAhAQKEBAgQEiBASIAAIQEChAQI\nEBIgQEiAACEBAoQECBASIEBIgAAhAQKEBAgQEiBASIAAIQEChAQIEBIgQEiAACEBAoQECBAS\nIEBIgAAhAQKEBAgQEiBASIAAIQEChAQIEBIgQEiAACEBAoQECBASIEBIgAAhAQKEBAgQEiBA\nSIAAIQEChAQIEBIgQEiAACEBAoQECBASIEBIgAAhAQKEBAgQEiBASIAAIQEChAQIEBIgQEiA\nACEBAoQECBASIEBIgAAhAQKEBAgQEiBASIAAIQEChAQIEBIgQEiAACEBAoQECBASIEBIgAAh\nAQKEBAgQEiBASIAAIQEChAQIEBIgQEiAACEBAoQECBASIEBIgAAhAQKEBAgQEiBASIAAIQEC\nhAQIBITUAJmZsMv14RgOASgREiBASIAAIQEChAQIEBIgQEiAACEBAoQECBASIEBIgAAhAQKE\nBAgQEiBASIAAIQEChAQIEBIgQEiAACEBAoQECBASIEBISf33zf04kBohJfRfm/vBIClCSua/\nO+7Hg5QIKZX7jkipaISUyIOOKKlghJTGw44oqVyElERHR5RULEJKobMjQioVISXQ3REllYqQ\nEngWEiWViZD0nnZESWUiJLmejiipSIQkR0g1IiS13o4oqUShIb2/bM9/t3a7e091CL8BIVFS\neQJDOq5bfwN6k+QQCzCkI0IqT2BIu2b1ejh/9fG2anYpDrEAg0KipOIEhrRqDj9fH5pVikP4\nDeuIkooTGFLTdP1Ddgi7oR0RUml4RpIaHBIlFSb2d6S3j/NX5f6OREi1irz8vWldtVsfkxzC\nbHhHlFSY2NeRdufXkVbbl0JfRxoTEiUVhTsbhEZ1REhFISShcSFRUkkiQzr+a5rN2/c3KfDy\n98iOCKkkkbcIrS432l2+CSFRUklCL3/vv2rar8632RUY0uiOCKkgoS/Inv/Xx2r9QUiUVBjD\nLULHzYaQCKkwgSGtm+uLsOtNgSFN6IiSyhEY0r759/3VR7MhJEIqSuTl791PPW8NIRFSUUJf\nkD1sr199/CstpEkdUVIxuLNBZGJIlFQIQtKY2hEhFYJPEdKYHBIllYFPEZKY3hEhlWE5nyLU\ntE08hA0h1Y7PbJCYERIlFYFPEZIgpNrxjKQwpyNKKgKfIqRASNXjU4QU5oVESQXgU4QEZnZE\nSAXgzgYBQgIhzTe3I0oqACHNR0iwhVTS60jzQ6Kk7BHSfIQEfrQTICQQ0nyCjigpe4Q0GyGB\nN/bNJ+mIkHLHG/vm0oRESZlbzhv7FIdwICR88jaK2UQdUVLmeGPfTISEE56RZpKFRElZ4419\nMxESTnhj30yEhBPe2DePriNCyhp3NswjDImSckZI8xASzghpFmVHlJQzQpqFkHBBSLMQEi4I\naQ5tR5SUMUKag5DwjZDmICR8I6Q51CFRUrYIaQ5CwjdCmkHeESFli5Bm0IdESbkipBkICVeE\nNF2CjggpV4Q0XYqQKClThDQdIeEHIU1HSPhBSJMl6YiSMkVIkxESfhHSZISEX4Q0WaKQKClL\nhDQZIeEXIU2VqiNCyhIhTZUsJErKESFNRUhoIaSpCAkthDRRuo4oKUeENBEhoY2QJiIktBHS\nNCk7oqQMEdI0hIQbhDQNIeEGIU1DSLhBSNOkDYmSskNI0xASbhDSJIk7oqTsENIkhIRbhDQJ\nIeEWIU2SPCRKygwhTZG+I0LKDCFNQUi4Q0hTEBLuENIUASFRUl4IaYKIjggpL4Q0ASHhHiFN\nEBISJWWFkCYgJNwjpAliQqKknBDSBISEe4Q0XlBHhJQTQhqPkPAHIY1HSPiDkMaLComSMkJI\no4V1REgZIaTR4kKipHwQ0miEhL8IaTRCwl+ENFZgR5SUD0Iai5DwACGNRUh4gJDGIiQ8QEhj\nhYZESbkgpLEICQ8Q0kixHRFSLghppOCQKCkThDQSIeERQhqJkPAIIY0UHRIl5YGQxgnviJDy\nQEjjEBIeIqRxCAkPEdI48SFRUhYIaRRDR4SUBUIahZDwGCGN4giJknJASKMQEh4jpFEICY8R\n0iiWkCgpA4Q0CiHhMUIaw9MRIWWAkMYgJHQgpDFMIVHS8hHSGISEDoQ0gqsjQlo+QhrBFhIl\nLR4hjUBI6BIf0n7dNNu3pIdIhZDQJTCk5vwfbpqzXZJDJOYLiZKWLjqkXbM7fn5+7Jp9ikMk\nRkjoEh3Sqjmevj426xSHSIyQ0CU6pKZp/UN+iLSMHRHS0kWH9O8a0irFIdIiJHQKDWn7sn9r\nXr++PO6eX20gJErKTGhIF+cvV8cUh0iLkNAp8nWkw2G/327Plxx2TzsiJELKDXc2DGbtiJIW\njpAGIyR0I6TBCAndCGkwc0iUtGiENJS7I0JaNMPl79+r4PJDpOTuiJAWLTCkPSERUrFCX0da\nbVIfIiF3R5S0aKG/Ix163oYkOEQ67owIadFiLzbsm0PqQyTjzoiQFo2rdkO5M/qPkpaMkAZy\nR3TingG6EdJA7ohO3DNAt9CQ3l+25yvf2917qkMk447ozD0EdAoM6bhuvYr0/EI4IT3mHgI6\nBYa0a1avl4t2H2+rvxfCB79a6+Fu6Mw9BHQKDGnVuvZ9yO0zG9wJXbingE7hHxD56B+yQ6Tj\nTuibewzowjPSMO6CvrnHgC6xvyO9fZy/evg7kuIQ6bgL+uYeA7pEXv7etK4mrDP7FCF3QVfu\nOaBD7OtIu/PrSKvtS26vI7n7+eEeBDpwZ8Mg7n5+uAeBDoQ0iLufH+5BoAMhDeLu54d7EOjg\nCimz15Hc/fxyTwKPEdIg7nx+uSeBx/jRbhB3Pr/ck8BjhDSIO58W9yjwECEN4Y6nzT0LPMQb\n+4Zwx3PDPQw8whv7hnC3c8M9DDyynDf2KQ6RirudG+5h4BHeRjGAO51b7mngEd7YN4A7nVvu\naeARnpEGcKdzxz0OPMAb+wZwl3PHPQ48wBv7BnCXc8c9DjzAG/v6ucP5wz0Q/MWdDf3c3fzh\nHgj+IqR+7m7+cA8EfxFSP3c3f7kngj8IqZ87m7/cE8EfhNTPnc1f7ongD0Lq587mL/dE8Ach\n9XNn84B7JLhHSL3c0TzingnuEVIvdzSPuGeCe4TUyx3NQ+6h4A4h9XI385B7KLhDSL3czTzk\nHgruEFIfdzId3GPBLULq4y6mg3ssuEVIfdzFdHCPBbcIqY+7mA7useAWIfVxF9PFPRfcIKQe\n7l46uQeDG4TUw91LJ/dgcIOQerh76eaeDNoIqYc7l27uyaCNkHq4c+nmngzaCKmHO5cn3KNB\nCyH1cNfyhHs0aCGkHu5annCPBi2E9Jw7lmfcs0ELIT3njuUZ92zQQkjPuWN5yj0c/CKk59yt\nPOUeDn4R0nPuVp5yDwe/COkpdyo93OPBD0J6yl1KD/d48IOQnnKX0sc9H1wR0lPuUPq454Mr\nQnrKHUof93xwRUhPuUPp454PrgjpGXcnvdwDwhUhPePupJ97QvhGSM+4M+nnnhC+EdIz7kz6\nuSeEb4T0jDuTAdwjwgUhPeOuZAD3iHBBSM+4KxnAPSJcENIT7kgGcQ8JZ4T0hLuRQdxDwtnM\nkNYvH7KH0nEII3cjg7iHhLOZITVNk6IlQhrOPSSczQzp+PovRUuENIJ7SjgR/I70/rJWt0RI\nI7inhBPNxYbD6ut5aT//0Tw5hIG7kIHcY8KJJKS3TXOyETyerkM4uAsZyj0nfCpCOr58PR2t\n345fNW01j4mQxnHPCZ/zQ3o/XWzYHS7/B9n+J6Qx3HPC5/zXkb6ejPbH6/9hpXhE94fwcQcy\nmHtQmP860vZN9lA6DmHjzmM496Qw/3Uk2QPpPISNO4/h3JOC4M6G7y9Wsh/r7g9h485jOPek\nIAvpQ3eh4f4QNu48RnCPCnNCemva1uZHpeeuYwT3qDDrGWnd7ujd/Kj03HWM4B4VZL8jaRHS\nSO5ZgTf2dXLHMYZ7VpgR0unZqPXDnflRybnbGMU9LBBSF3cb47inVT1+tOviTmMc97SqR0hd\n3GmM455W9Wb9aHfD/Kjk3GmM455W9Qipg7uMsdzzqh0/2nVwhzGWe161I6QO7jDGcs+rdlz+\n7uAOYzT3wCpHSB3cXYzmHljl+NGug7uL8dwTqxshPeauYgL3yOo2O6TX7dePdepPbiCkCdwj\nq9vckDbfvyHJPtLu7yEs3FVM4B5Z3WaGtGtWpyejt1XzonpE94fwcFcxhXtmVZsZ0qq5fDbk\nQfeZdveH8HBHMYV7ZlVTvUO2tMvf7iimcM+sarN/tLs+I+00j+fvITzcUUziHlrN5l5s2J5/\nR3pf/RM9ngeHsHA3MYl7aDXj7u+H3ElM455azQjpIXcSE7nHVjHubHjIXcRE7rFVjJAechcx\nkXtsFVOF9C69tYGQpnGPrWJzQ9qN+R3p/WV7uZ9o1/MBx+6Q3EFMZp5bxWa/jnTVf9vqsf1Z\n4c//cDMhTWSeW8Vm3yL0+rlpPj42Tf+H6O+a1evl5duPt9XzF3AJaSLz3ComuEXo5evZ6NDz\nFHNyvS/vpOfePEKayjy4eglCemv2g+61a+7/Q+mjknLnMJ15cPWaGdL260e7j2b9+T4gJJ6R\nApgHV6+ZIb2dAjq/ua//ZrvTe5c+zl8t/Xckdw0zeAdXsbmXv19O//rXDLr5e9O6ard++vfQ\nCWky7+TqFXpnw/vu/DrSavuy7NeR3DHM4Z1cvbhF6AF3DHN4J1cvPkXoAXcMs3hHVy0+RegB\ndwuzeEdXLT5F6AF3C/N4Z1crPkXoAXcK83hnV6vATxEa8Y5ab0juEmayzq5agZ8itCekENbZ\nVSvyU4QOq/47Wx8cIpy7hLmsw6tV6IefDP70O0Kawzq8WsV+itC+dd+q+FEJuUOYyzq8WnFn\nwx/uDuZzTq9WhPSHO4P5nNOr1fxbhE73NmxfRQ/n4SGCuTMQcI6vUqpbhAZdj8vjU4TcFQg4\nx1epmSHtf24R2vf+d7l8ipC7AgHn+Co1M6T1zwuy697/rudThJJ9kPhY7goEnOOrVOAtQpl8\nZoM7Agnj/Cole0bqv2k1k08RcjcgYZxfpQJ/R+IZKY5xfpUKvGqXyacIuRvQMA6wTqK3mg96\nHSmPTxFyJ6BhHGCd+BShe+4ENIwDrNPcT1qV/jHzh4eI5k5AxDjBKqkuf2sR0mzGCVZp9uXv\np7/rTGUMyR2Aim+CdZoZ0nG76f/DSA+/yVJfR3IHIOMbYZVm/2g38a4eQkrNN8IquUIafohg\n7v0v4xthlXhj3x33/tfxzbBGhHTLvfuFbDOskubOhn/DPkQ/gzf2uXe/km2INQr8EP0s3tjn\n3vxKtiHWSPYh+v13f/e8sW/moxJxb34l2xBrJPsQ/f53yGbxNgr35leyDbFGoR+i3/UPxaPS\ncO99LdcUayT7EP3+X5JyeEZyb30t1xRrNPdiw8v3h+gX8sY+99bXck2xRro7G/rvbsjgjX3u\nrS/mGmOFIkPK4I197p0v5hpjhbiz4YZ756u55lgfQrrh3vhqrjnWh5BuuDe+mmuO9SGkG+6N\nr+aaY30Iqc297/VMg6wPIbW5t72eaZD1IaQ297ZPwDTJ6hBSm3vXJ2CaZHUIqc296xMwTbI6\nhNTi3vRJeEZZHUJqce/5JDyjrA4htbj3fBKeUVaHkFrcez4NzyxrQ0gt7i2fhmeWtSGkFveW\nT8Mzy9oQUot7yyfiGWZlCKnFveMT8QyzMoT0y73hk7FMszKE9Mu935OxTLMyhPTLvd+TsUyz\nMoT0y73fk7FMszKE9Mu939OxjLMuhPTDvdsTcoyzMoT0w73bU3LMsy6E9MO92VNyzLMuhPTD\nvdmTcgy0KoT0w73Xk3IMtCqEdOXe6mkZBloXQrpyb/W0DAOtCyFdubd6YoaJVoWQrtw7PTHD\nRKtCSFfunZ6aYaQ1IaQr90ZPzTDSmhDSlXujJ2eYaUUI6cq9z5MzzLQihPTNvc3Ti59pTQjp\nm3ubpxc/05oQ0jf3Ng8QP9SKENI39y4PED/UihDSN/cujxA/1XoQ0oV7j4cIn2pFCOnCvcdj\nhI+1HoR04d7iMcLHWg9CunBv8RjhY60HIV24t3iQ8LlWg5Au3Ds8SPhcq0FIZ+4NHiV6rvUg\npDP3Bg8TPdhqENKZe3+HiR5sNQjpzL2/w0QPthqEdObe33GiJ1sLQjpzb+840ZOtBSGdubd3\noOjRVoKQTtybO1LwaGtBSCfuzR0qeLaVIKQT994OFTzbShDSiXtvhwqebSUI6cS9t2MFD7cO\nhPRZW0eElAIhfVYXEiUlQEifhIT5COmzvpAoSY+QPgkJ8xFShR0Rkh4h1RgSJckREiFBgJCq\nDImS1AiJkCBASHWGRElihERIECAkQoIAIdXZESWJERIhQYCQag2JkqQIiZAgQEjVhkRJSoTk\n3s4+gUMuHyG5t7NP4JDLR0ju7WwUOOXiEZJ7NxsFTrl41Yfk3sxWcWMuHiHVLG7MxSOkqsXN\nuXSEVLW4OZeOkKoWN+fSxYe0XzfN9i3pIcZwb2WzuEEXLjCk5vwfbpqzXZJDTODeyW5xky5b\ndEi7Znf8/PzYNfsUh5jAvZHd4iZdtuiQVs3x9PWxWac4xHjufewXNuqyRYfUNK1/yA8xnnsb\n+4WNumzRIf27hrRKcYjx3NvYL2zUZQsNafuyf2tev7487p5fbSCkQGGzLlpoSBfnL1fHFIcY\nzb2JlyBq1mWLfB3pcNjvt9vzJYfd044IKVTUsItW+Z0N7j28DFHTLhkhgZAEQkN6f9mef0va\n7t5THWIk9xZehqhplywwpOO6+bVJcoix3Dt4KYLGXbLAkHbN6vVw/urjbbWMy9/uDbwUQeMu\nWWBIq+bw8/VhGS/IujfwYgTNu2Dhd38/+ofsEGO59+9yBA28XDwj4SRo4OWK/R3p7eP8Fb8j\nLU7QwMsVefl707pqt17ELULu7bsgQRMvVuzrSLvz60ir7csyXkdyb94liZl4uaq+s8G9eRcl\nZuTFIiR8i5l5qaq+Rci9dZclZualqvkWIffOXZiQmRdrObcINW0TDzGOe+cuTcjQS1XzC7Lu\njbs4IVMvVM23CLn37eKETL1QPCPhV8jYy1TxLULuXbtAEWMvVMW3CLl37RJFzL1MFd8i5N60\nSxQx9zJVfGeDe9MuUsTgi0RIaIsYfJEICTciJl8iV0j+15HcO3ahAiZfJELCrYDRl6jeH+3c\nG3axAmZfIELCnYDZF4iQcCdg9gWq9o197u26YOmHX6Bq39jn3q1Lln765VnOG/sUhxjBvVmX\nLP30y1Pt2yjcm3XR0o+/OLW+sc+9VZct+fjLU+szknurLlzy+Ren1jf2uXfq0iVfgNLU+sY+\n90ZduuQLUJpa39jn3qiLl3wFClPrnQ3ufbp4yVegMJWG5N6mGUi9BIUhJHRIvQZlISR0Sb0I\nRSEkdEm9CEWpMyT3Fs1E4lUoCiGhU+JVKAohoVviZSgJIeGJxOtQkCpDcm/PfKRdh5IQEp5J\nuxAFISQ8lXYlykFIeC7tUhSDkPBc2qUoBiGhR9q1KEWNIbl3ZmaSrkUxCAl9ki5GKQgJvZKu\nRiEICb2SrkYhKgzJvS0zlHI5CkFIGCDlepSBkDBEygUpAiFhiJQLUoT6QnJvyUwlXJEiEBKG\nSbgkJSAkDJRwTQpASBgo4ZoUgJAwVMJFyV91Ibl3Y87SrUr+CAnDpVuW7BESRki3LrmrLST3\nTsxdsoXJHSFhjGQLkztCwijJViZzhIRxki1N3ioLyb0LS5BqbfJGSBgp1drkjZAwVqrFyRoh\nYbxUy5MxQsJ4qZYnY3WF5N6AxUi0PhkjJEyRaIHyRUiYJNEKZauqkNybryhplihbhISJ0qxR\nrggJU6VZpEwREiZLs0p5IiRMl2aZslRTSO5tV6Ak65QlQsIcSRYqR4SEWZKsVIYqCsm95QqV\nYqkyREiYKcVa5YeQMFeKxcpOPSG5t1vBEqxWdggJ8yVYrtwQEhQSLFheCAkKCRYsL9WE5N5p\nxdMvWVYICSL6NcsJIUFFv2gZqSUk9yarg3zZ8kFI0JEvWz4qCcm9w6qhXrhsEBKk1CuXC0KC\nmHrt8kBIUFMvXhbqCMm9tSojXr0sEBL0xMuXgypCcu+rCmkXMAOEhCS0K7h8hIREtGu4dDWE\n5N5RtZIu4tIREtKRLuOyERJSki7kklUQknsvVU65lAtGSEhMuZbLVX5I7n2EKlIiJAQQrudC\nFR+SewvhTLegC0VICKJb0yUiJETRLeoClR6Se/OgTbasy0NIiCRb2KUpPCT3vsFfqrVdFkJC\nONXqLkloSO8v2+Zku3tPdYhb7h2DDqL1XZDAkI7r5tcmySHuufcLOolWeDECQ9o1q9fD+auP\nt1WzS3GIO+7Ngqc0i7wQgSGtmsPP14dmleIQd9w7BX0067wEgSE1Tdc/ZIe45d4lGEKy1H4F\nPyO5dwiGUqy2W+zvSG8f569ifkdybw+MoFhwq8jL35vWVbv1MckhWtxbAyMJ1two9nWk3fl1\npNX2JeB1JPe+wHiCZXcp9s4G957ARPOX3qLUkNzbAXMItlC0Um8Rcm8FzDR7BwQr9BYh9zaA\nxNxtEKjMW4TcGwBCM/dCkCJfkHUvPeTmbYgAJd4i5F50pDJrW6RV4DOSe7WR2JzNkUx5twi5\nlxlRZmwSveJuEXKvLsJN3yxChd0i5F5TWE3dNgJF3dngXkcsiHZH9iooJPfKYdHUe/ROMbcI\nudcJeZq2kyVbdnG3CLnXArWQbdmp/8lZzy1CTVvX93BPEpVTVrH0F2SBeCXeIgSE4xkJECjv\nFiHAoLhbhACHwm4RAjwKurMB8CEkQICQAAFXSLyOhKIQEiDAj3aAACEBAoQECBTzxj7AqYQ3\n9gF2y3ljn+IQgAlvowAEeGMfIMAzEiDAG/sAAd7YBwgs9I19QGbG5+B8uijzqarIs+KkQr9Z\nRsdOp8iz4qRCv1lGx06nyLPipEK/WUbHTqfIs+KkQr9ZRsdOp8iz4qRCv1lGx06nyLPipEK/\nWUbHTqfIs+KkQr9ZRsdOp8iz4qRCv1lGx06nyLPipEK/WUbHTqfIs+KkQr9ZRsdOp8iz4qRC\nv1lGx06nyLPipEK/GVArQgIECAkQICRAgJAAAUICBAgJECAkQICQAAFCAgQICRAgJECAkAAB\nQgIECAkQICRAwBbSbtWsdk//ikVe2p++Xsa57a97o3U62Z/Z9aTky+UK6fJXYdamo+sdWitT\nxrkdrn+ToXU62Z/Z9aT0y2UK6b1ZHT4Pq6bnL8Lk49Bsr1+WcW5fJ3DZG63Tyf7Mfk5Kv1ym\nkHbN29f/fG1ePIfX2/+eShHntm821x98fk8n9zP7PSn9cplC2janv5vZ+v8Luds3++uXRZxb\ns7v+ie3W6eR+Zr8npV8uU0jf5zPlT6Mt07Z5+/f1K+vpyyLO7XB/Hqf/lfuZ/Z6UfrkISWN7\n+eV181nOuZUX0mcrJPVyEZJG07x+fh53p58YSjm3kkPSLxchKR1PV1FLObeSQ7pQLpdpJqv8\nl+Sh0wmVcm7fJ9A6nQLO7PaxC0/KetXuI9vrP11Oy1HKud1ctfv4vWqX9Zn9DUl0UqaQXs4X\n79+anefweqvmdI/JeTlKObfvPdc6nQLO7OdpVr1c3NmgsTstxPH84l4p51bgnQ0/J6VfLteP\nu+ufy49lOK7OJ3T+f2uFnNv1p6DW6eR/Zt8npV8uV0jH8y23poOncDqh9f7nywLO7RpS63Ty\nP7P2SUmXK+MLMMByEBIgQEiAACEBAoQECBASIEBIgAAhAQKEBAgQEiBASIAAIQEChAQIEBIg\nQEiAACEBAoQECBASIEBIgAAhAQKEBAgQEiBASIAAIQEChAQIEBIgQEiAACEBAoQECBASIEBI\ngAAhAQKEBAgQUp7m/jV7iLEeeSKkhWE98kRIC8N65ImQFob1yNMlpF3z4n4guCCkPJ1D2jUb\n9+PAN0LK0ykkOloQQsrTV0h0tCSElKem2TTNu/tR4Ach5alpmlWzdj8K/CCkPH09HR2a5tX9\nMHBFSHk6XWx4aVbuh4ErQsrT+fL3mpeRFoOQ8nQO6euHu6P7geCCkPJ0ubPhpdm6HwguCAkQ\nICRAgJAAAUICBAgJECAkQICQAAFCAgQICRAgJECAkAABQgIECAkQICRAgJAAAUICBAgJECAk\nQICQAAFCAgQICRAgJECAkAABQgIECAkQICRAgJAAAUICBAgJECAkQICQAIH/AemBAE38xwZt\nAAAAAElFTkSuQmCC",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mu <- 50\n",
    "sd <- 20\n",
    "mode <- 42\n",
    "r <- mu / sd^2\n",
    "s <- 1 + mode * r\n",
    "k <- seq(from = 0, to = 150, by = 0.1)\n",
    "plot(k, r^s * k^(s - 1) * exp(-r * k) / gamma(s), type = \"h\", col = \"skyblue\", ylab = \"probability\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Poisson distribution\n",
    "The Poisson discrete distribution is defined over integer values from 0 to $+ \\infty$ and can represent the number of occurrences of an event in a time interval.\n",
    "\n",
    "$$f(k, \\lambda) = \\frac{\\lambda^k e^{-\\lambda}}{k!}$$\n",
    "\n",
    "$\\lambda$ is a real number that represents the mean and the variance of the distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAM1BMVEUAAABNTU1oaGh8fHyH\nzuuMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD///8LNvnPAAAACXBIWXMAABJ0\nAAASdAHeZh94AAAatElEQVR4nO3d61oqybKF4SxYWCAC3v/VLkCdbR9SSSoyRkbU9/7Y2372\nHhbtzK9VTrO8A1isqG8AkAEhAQYICTBASIABQgIMEBJggJAAA4QEGCAkwAAhAQYICTBASIAB\nQgIMEBJggJAAA4QEGCAkwAAhAQYICTBASIABQgIMEBJggJAAA4QEGCAkwAAhAQYICTBASIAB\nQgIMEBJggJAAA4QEGCAkwAAhAQYICTBASIABQgIMEBJggJAAA4QEGCAkwAAhAQYICTBASIAB\nQgIMEBJggJAAA4QEGCAkwAAhAQYICTBASIABQgIMEBJggJAAA4QEGCAkwAAhAQYICTBASIAB\nQgIMEBJggJAAA4QEGCAkwAAhAQYICTBASIABQgIMEBJggJAAA4QEGCAkwAAhAQYICTBASIAB\nQgIMEBJggJAAA4QEGCAkwAAhAQYICTBASIABQgIMEBJggJAAA4QEGHg+pLf9rtzs5jfD2wOE\n9GxIl035y9b0JgHxPBvSXKbX0/2j83Eqs90NAiJ6NqSpnP58fCqTzY0Bono2pFJq/wCsEN+R\nAAMLfkc6nu8f8TsS8Pzd39tv99ptLpY3CYhnweNI8/1xpGm353EkrJ7D3QQFCOaJU24fjuAS\ngCXPkC4vpWyPn5/kx89CSAjGMaTLdP8WuPv4JISETBxDmsvhWtNhuj/NjpCQimNI08fwPG3O\nhIRkHEP6auey3RISknEMaVO+HoTdbAkJuTiGdCgvnx+dy5aQkIrn3d/zn3qO//EA1sJHtwAl\n1wdkT7uvj84vfEdCJjyzATBASIABQgIMqELiXjukQkiAAX60AwwQEmCAkAADriE9/N7fhIRg\nPF/Y9/h7fxMSgnF9Yd/D7/1NSAjG9YV9D7/TKiEhGMEL+/79D2aXGN7/xHv0wnckV4SUle/v\nSI++9zchddqjF8+7vx9/729C6rRHL76PIz363t+E1GmPXnhmgytCyoqQXBFSVoTkipCyIiRX\nhJQVIbkipKwIyRUhZUVIrggpK0JyRUhZEZIrQsqKkFwRUlaE5IqQsiIkV4SUFSG5IqSsCMkV\nIWVFSK4IKStCckVIWRGSK0LKipBcEVJWhOSKkLIiJFeElBUhuSKkrAjJFSFlRUiuCCkrQnJF\nSFkRkitCyoqQXBFSVoTkipCyIiRXhJQVIbkipKwIyRUhZUVIrggpK0JyRUhZEZIrQsqKkFwR\nUlaE5IqQsiIkV4SUFSG5IqSsCMkVIWVFSK4IKStCckVIWRGSK0LKipBcEVJWhOSKkLIiJFeE\nlBUhuSKkrAjJFSFlRUiuCCkrQnJFSFkRkitCyoqQXBFSVoTkipCyIiRXhJQVIbkipKwIyRUh\nZUVIrggpK0JyRUhZEZIrQsqKkFwRUlaE5IqQsiIkV4SUFSG5IqSsCMkVIWVFSK4IKStCckVI\nWRGSK0LKipBcEVJWhOSKkLIiJFeElBUhuSKkrAjJFSFlRUiuCCkrQnJFSFkRkitCyoqQXBFS\nVoTkipCyIiRXhJQVIbkipKwIyRUhZUVIrggpK0JyRUhZEZIrQsqKkFwRUlaE5IqQsiIkV4SU\nFSG5IqSsCMkVIWVFSK4IKStCckVIWRGSK0LKipBcEVJWhOSKkLIiJFeElBUhuSKkrAjJFSFl\nRUiuCCkrQnJFSFkRkitCyoqQXBFSVoTkipCyIiRXhJQVIbkipKwIyRUhZUVIrggpK0JyRUhZ\nEZIrQsqKkFwRUlaE5IqQsiIkV4SUFSG5IqSsCMkVIWVFSK4IKStCckVIWRGSK0LKipBcqUMi\nxF78QzpsStkdu15iXOoQCKkXx5DKfbgtd3OXSwxPHQIh9eId0lzmy/v7eS6HHpcYnjoEQurF\nO6SpXG4fX8qmxyWGpw6BkHrxDqmUb/9gfonhqUMgpF68Q3r5CmnqcYnhqUMgpF5cQ9rtD8fy\nev3wMv98bwMhDbpHjWtIH+4fTpcelxieOgRC6sXzcaTT6XDY7e53Ocw/dkRIo+5RwzMbXKlD\nIKReCMmVOgRC6oWQXKlDIKReCMmVOgRC6oWQXKlDIKReBHd//3UvuPklhqcOgZB6cQzpQEjy\nEAipF9fHkaZt70uMTh0CIfXi+jvS6ZeXIRlcYnDqEAipF987Gw7lVP+0j/7cF5k6BELqhXvt\nXKlDIKReCMmVOgRC6oWQXKlDIKReCMmVOgRC6kUVEo8jhdyjhpBcqUMgpF740c6VOgRC6oWQ\nXKlDIKReCMmVOgRC6sU1pLf97v68hd381usSnakPsnqPGseQLptvzwH6+emrhDToHjWOIc1l\nev14qt35OAV9Xzv1QVbvUeMY0vTtGaunoO+0qj7I6j1q3P9al//6B7NLdKc+yOo9aviO1ER9\nkNV71Pj+jnQ83z/id6Soe9R43v29/Xav3Sbme3+rD7J6jxrfx5Hm++NI027P40gx96jhmQ1N\n1AdZvUcNITVRH2T1HjWE1ER9kNV71BBSE/VBVu9RQ0hN1AdZvUcNITVRH2T1HjWE1ER9kNV7\n1BBSE/VBVu9RQ0hN1AdZvUcNITVRH2T1HjWE1ER9kNV71BBSE/VBVu9RQ0hN1AdZvUcNITVR\nH2T1HjWE1ER9kNV71BBSE/VBVu9RQ0hN1AdZvUcNITVRH2T1HjWE1ER9kNV71BBSE/VBVu9R\nQ0hN1AdZvUcNITVRH2T1HjWE1ER9kNV71BBSE/VBVu9RQ0hN1AdZvUcNITVRH2T1HjWE1ER9\nkNV71BBSE/VBVu9RQ0hN1AdZvUcNITVRH2T1HjWE1ER9kNV71BBSE/VBVu9RQ0hN1AdZvUcN\nITVRH2T1HjWE1ER9kNV71BBSE/VBVu9RQ0hN1AdZvUcNITVRH2T1HjWE1ER9kNV71BBSE/VB\nVu9RQ0hN1AdZvUcNITVRH2T1HjWE1ER9kNV71BBSE/VBVu9RQ0hN1AdZvUcNITVRH2T1HjWE\n1ER9kNV71BBSE/VBVu9RQ0hN1AdZvUcNITVRH2T1HjWE1ER9kNV71BBSE/VBVu9RQ0hN1AdZ\nvUcNITVRH2T1HjWE1ER9kNV71BBSE/VBVu9RQ0hN1AdZvUcNITVRH2T1HjWE1ER9kNV71BBS\nE/VBVu9RQ0hN1AdZvUcNITVRH2T1HjWE1ER9kNV71BBSE/VBVu9RQ0hN1AdZvUcNITVRH2T1\nHjWE1ER9kNV71BBSE/VBVu9RQ0hN1AdZvUcNITVRH2T1HjWE1ER9kNV71BBSE/VBVu9RQ0hN\n1AdZvUcNITVRH2T1HjWE1ER9kNV71BBSE/VBVu9RQ0hN1AdZvUcNITVRH2T1HjWE1ER9kNV7\n1BBSE/VBVu9RQ0hN1AdZvUcNITVRH2T1HjWE1ER9kNV71BBSE/VBVu9RQ0hN1AdZvUfNwpA2\n+7PZTalcYijqg6zeo2ZhSKWUHi0R0qB71CwM6fL60qMlQhp0jxqD35He9hvrlghp0D1qbO5s\nOE3X70uH5bfmh0uMQX2Q1XvUmIR03JabrcHtqV1iEOqDrN6jZnlIl/3129HmeLnWtLO5TYQ0\n7B41S0N6u93ZMJ8+/g9m55+QBt2jZunjSNdvRofL1/9hsrhF/7zEUNQHWb1HzdLHkXZHs5tS\nucRQ1AdZvUfN0seRzG5I9RJDUR9k9R41i5/Z8PnBZPZj3T8vMRT1QVbvUWMU0tnujoZ/XmIo\n6oOs3qNmQUjH8t1GfKt8qA+yeo+aJd+RNt87ehPfKh/qg6zeo8bqdyRbhDToHjW8sK+J+iCr\n96hZENLtu9G3H+7Et8qH+iCr96ghpCbqg6zeo4Yf7ZqoD7J6jxpCaqI+yOo9ahb9aPc34lvl\nQ32Q1XvUEFIT9UFW71HDj3ZN1AdZvUcNITVRH2T1HjXc/d1EfZDVe9QQUhP1QVbvUcOPdk3U\nB1m9R41rSG/73f2b127+5bnihDToHjWLQ3q9tfHQOzdcvr/s4uf3wCOkQfeoWRrS9jOMB97S\nbi7T68f7dp2PU5mNb5UP9UFW71GzMKRrHLdvRtcw9r/upnL68/Hp57fuIqRB96hZGNJXHL+E\n8bF7/BWBhDToHjVWr5B94O5vviPF36Nm8Y92X9+Rfvyd5/P/dzp+/N0v/I601n1eS+9s2N1/\nR3qbXh4Ybr/da7f58a0lCSnpPi/XZ3+/zffHkabdnseR1rnPi5dRNFEfxOj7vHiKUBP1QYy+\nz4uQmqgPYvR9XlYhvZn9bX3VS4xAfRCj7/NaGtLM70jsYfA40pffn7bacOcEISXd57X4KUKv\n79tyPm/L72+ifyCk1e/zMniK0P763ej0y+si7k7TA/9Pz94qH+qDGH2fl0FIx3J48K+leOSJ\nRM/eKh/qgxh9n9fCkHbXH+3OZfP+9tidDYdvz1v916ft9eiuJfVBjL7Pa2FIx9uZvz+H7pEn\n2z11iaGoD2L0fV5L7/7e3/7ppTz6M9szlxiJ+iBG3+fFMxuaqA9i9H1ehNREfRCj7/NyfBeh\nv38SHkda4z4vx3cR+vsnIaQ17vNyfBehJy8xFPVBjL7Py/FdhJ68xFDUBzH6Pi/HdxF68hJD\nUR/E6Pu8HN9FiPf+Zp+X47sI8d7f7PNyfPMT3vubfV6OIfFOq+zzcnxmA+/9zT4vx5D4jsQ+\nr+VPEbo9t2H3+sCO9/5mn5fVU4QeeRE57/29+n1eC0M6/HmK0OGBJe/9vfZ9XgtD2vx5QHZj\nc3v+fYmhqA9i9H1ePEWoifogRt/nZfYdiSetsl8z19+RnrrEUNQHMfo+L8977Z68xEjUBzH6\nPi+jl5o/8jjS05cYiPogRt/nxZufNFEfxOj7vJa+06rt+9n91yWGoj6I0fd5Wd39bYuQku7z\nWnz3949P9XkWISXd57UwpMtu+/tfjLTsEkNRH8To+7wW/2jHX33JHoTUSH0Qo+/z4u7vJuqD\nGH2fFyE1UR/E6Pu8bJ7Z8NL8JvotlxiI+iBG3+elehP9hkuMRH0Qo+/zMnsTfZ79zX7NzN5E\nn1fIsl8zXiHbRH0Qo+/zMnsTfdNfkggp6T6vpXc27D/fRJ8X9rFfNbtnNlg+u4GQku7zIqQm\n6oMYfZ8Xz2xooj6I0fd5EVIT9UGMvs+LkJqoD2L0fV6E1ER9EKPv8yKkJuqDGH2fFyE1UR/E\n6Pu8CKmJ+iBG3+dFSE3UBzH6Pi9CaqI+iNH3eRFSE/VBjL7Pi5CaqA9i9H1ehNREfRCj7/Mi\npCbqgxh9nxchNVEfxOj7vAipifogRt/nRUhN1Acx+j4vQmqiPojR93kRUhP1QYy+z4uQmqgP\nYvR9XoTURH0Qo+/zIqQm6oMYfZ8XITVRH8To+7wIqYn6IEbf50VITdQHMfo+L0Jqoj6I0fd5\nEVIT9UGMvs+LkJqoD2L0fV6E1ER9EKPv8yKkJuqDGH2fFyE1UR/E6Pu8CKmJ+iBG3+dFSE3U\nBzH6Pi9CaqI+iNH3eRFSE/VBjL7Pi5CaqA9i9H1ehNREfRCj7/MipCbqgxh9nxchNVEfxOj7\nvFYWkvogrX2fFyGxd9znRUjsHfd5ERJ7x31ehMTecZ8XIbF33OdFSOwd93kREnvHfV6ExN5x\nnxchsXfc50VI7B33eRESe8d9XoTE3nGfFyGxd9znRUjsHfd5ERJ7x31ehMTecZ8XIbF33OdF\nSOwd93kREnvHfV6ExN5xnxchsXfc50VI7B33eRESe8d9XoTE3nGfFyGxd9znRUjsHfd5ERJ7\nx31ehMTecZ8XIbF33OdFSOwd93kREnvHfV6ExN5xnxchsXfc50VI7B33eRESe8d9XoTE3nGf\nFyGxd9znRUjsHfd5ERJ7x31ehMTecZ8XIbF33OdFSOwD7cdFSOwD7cdFSOwD7cdFSOwD7cdF\nSOwD7cdFSOwD7cdFSOwD7cdFSOwD7cdFSOwD7cdFSOwD7cflH9JhU8ru2PUSdeqDwD4rx5DK\nfbgtd3OXS/xKfRDYZ+Ud0lzmy/v7eS6HHpf4lfogsM/KO6SpXG4fX8qmxyV+pT4I7LPyDqmU\nb/9gfolfqQ8C+6y8Q3r5CmnqcYlfqQ8C+6xcQ9rtD8fyev3wMv98bwMhsQ/GNaQP9w+nS49L\n/Ep9ENhn5fk40ul0OOx297sc5h87IiT20fDMBvaB9uMiJPaB9uNyDeltv7v/lrSb33pd4hfq\ng8A+K8eQLpvyl22XS/xKfRDYZ+UY0lym19P9o/Nx4u5v9qk4hjSV05+PTzwgyz4V92d//9c/\nmF3iV+qDwD4rviOxD7Qfl+/vSMfz/SN+R2Kfjefd39tv99pt/vXUhvLds5f4jfogsM/K93Gk\n+f440rTb8zgS+1x4ZgP7QPtxERL7QPtxERL7QPtxqULicST2qRAS+0D7cfGjHftA+3EREvtA\n+3EREvtA+3Hxwj72gfbj4oV97APtx8UL+9gH2o+Ll1GwD7QfFy/sYx9oPy6+I7EPtB8XL+xj\nH2g/rnFe2Gdyid+oDwL7rHhhH/tA+3HxzAb2gfbjIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbj\nIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbj\nIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbj\nIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbj\nIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbjIiT2gfbj\nIiT2gfbjIiT2K9r3Q0jsV7Tvh5DYr2jfDyGxX9G+H0Jiv6J9P4TEfkX7fgiJ/Yr2/RAS+xXt\n+yEk9iva90NI7Fe074eQ2K9o3w8hsV/Rvh9CYr+ifT+ExH5F+34Iif2K9v0QEvsV7fshJPYr\n2vdDSOxXtO+HkNivaN8PIbFf0b4fQmK/on0/hMR+Rft+CIn9ivb9EBL7Fe37IST2K9r3Q0js\nV7Tvh5DYr2jfDyGxX9G+H0Jiv6J9P4TEfkX7fgiJ/Yr2/RAS+xXt+yEk9iva90NI7Fe074eQ\n2K9o30+wkNR/EOxj7/shJPYr2vdDSOxXtO+HkNivaN8PIbFf0b4fQmK/on0/hMR+Rft+CIn9\nivb9EBL7Fe37IST2K9r3Q0jsV7Tvh5DYr2jfDyGxX9G+H0Jiv6J9P4TEfkX7fgiJPXsDhMSe\nvQFCYs/eACGxZ2+AkNizN0BI7NkbICT27A0QEnv2BgiJPXsDhMSevQFCYs/eACGxZ2+AkNiz\nN0BI7NkbICT27A0QEnv2BgiJPXsDhMSevQFCYs/eACGxZ2+AkNizN0BI7NkbcA3pbb8rN7v5\n7clLqL+Q7Ne9r3MM6bIpf9k+dwn1F5L9uvd1jiHNZXo93T86H6cyP3UJ9ReSPfv/5hjSVE5/\nPj6V6alLjPuFZL/uvWNIpdT+4fFLjPuFZL/uPd+R2LM32Pv+jnQ83z/idyT22faed39vv91r\nt7k8dYlxv5Ds1733fRxpvj+ONO32PI7EPteeZzawZ2+wJyT27A32PEWIPXuDPU8RYs/eYM9T\nhNizN9jzgCx79gZ7niLEnr3Bnu9I7Nkb7HmKEHv2BvtxniJUvqt9inG/kOzXvQ/2FCFgTMGe\n2QCMiZAAA4QEGFCF9OTjSMCYCAkwMOiPdkAwT5xy+3BCXJvrc33T6xMS1+f60k/28Av7Olzb\nBtfn+vJP1vDCPvNrW+H6XF/+yRpe2Gd+bStcn+vLP1nDyyjMr22F63N9+SdreGGf+bWtcH2u\nL/9kfEfi+lzf4JM1vLDP/NpWuD7X13+yx9/72/7aRrg+1x/gkz38wr4O17bB9bn+mJ8s0LW5\nPtdPExKQBiEBBggJMEBIgAFCAgwQEmCAkAADhAQYICTAACEBBggJMEBIgAFCAgwQEmCAkAAD\nhAQYkIU0T2Wan3uJuonDRnv9qzflf8VOL6W8nGWXv0j//A9fX3m7W6H6s/x4y4eN6Oq3N2+5\nv05eWdJlEoZ01P77n6eP62tKPn29f5zhKRT9Wb6V6fR+msqzb/ew1Km8XG7/YXoRXf9m9+Tb\nAZqYrl//y+7Jt39a7uV+5Vnz9b+eu4+vvOUpFP1ZzuV4/Z+vZa+5/PUQ3/+X8ii/PvXX8Jhd\n/HaQL0++IeFyRfj1P5Tt53UtT6Hoz3JXbt/UT2WnufwX4VE+//njVHj59v6eCp8/1UpCvv43\n5PMrb3kKRX+Wyv8i/eXy7N+jYWBbzsJ/+01530/3H2819p8/2il+Ijn98/iZ/DmsOqTD/Vu7\nxL68Kv/tS/l4U0LZDTjc7m2YDqKrE5Kp8yT7yfL+44Q0pNudDS+y31Gv/yG5UV2ekCxdJt0P\ndpvbHc/SkG6/I51ljz8cbj/aXUMWfUvKE9I0QEhb3aNYL/efKaUhff9f/jbl9uvZRRXy57+3\n5SmU3mt3Ft5rd95sdY/rL/l76E2o7/4Xh/y3e+1sTqHo32R//0/yUfaA4PXSup/rBgjp4+t/\nln0RPr4XyB7H+vy6W57ClT6zQXeEvpE+irW53H5HeRVdfy63Z7jNqv+Q5nlmw/tmyV+HvtyL\n+DvCnfLqe+3X//NZbqrrf33lDU+h6s/y49m/oot/+9lKdgvexXe1HLfKr//n865VF//6yhue\nQl6PBBggJMAAIQEGCAkwQEiAAUICDBASYICQAAOEBBggJMAAIQEGCAkwQEiAAUICDBASYICQ\nAAOEBBggJMAAIQEGCAkwQEiAAUICDBASYICQAAOEBBggJMAAIQEGCAkwQEiAAUICDBASYICQ\nAAOEBBggpJi0f9cg/oU/j5gIaTD8ecRESIPhzyMmQhoMfx4xfYQ0l736huADIcV0D2kuW/Xt\nwCdCiukWEh0NhJBiuoZERyMhpJhK2Zbypr4V+IOQYiqlTGWjvhX4g5Biun47OpXyqr4Z+EJI\nMd3ubNiXSX0z8IWQYrrf/b3hYaRhEFJM95CuP9xd1DcEHwgppo9nNuzLTn1D8IGQAAOEBBgg\nJMAAIQEGCAkwQEiAAUICDBASYICQAAOEBBggJMAAIQEGCAkwQEiAAUICDBASYICQAAOEBBgg\nJMAAIQEGCAkwQEiAAUICDBASYICQAAOEBBggJMAAIQEGCAkwQEiAAUICDPwf9Y2eGBIU8kIA\nAAAASUVORK5CYII=",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "k <- seq(from = 0, to = 10, by = 1)\n",
    "lambda <- 4\n",
    "plot(k, lambda^k * exp( - lambda) / factorial(k), type = \"h\", col = \"skyblue\", ylab = \"probability\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gaussian distribution\n",
    "$$f(y|\\mu, \\sigma) = \\frac{1}{\\sqrt{2 \\pi} \\sigma} \\exp^{ - \\frac{1}{2} \\frac{(y - \\mu)^2}{\\sigma^2} }$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exponential distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Student's t distribution\n",
    "The t distribution is used to compute the true mean and the its standard deviation of a population that is normally distributed. When we compute the mean of a normal distribution we still have to compute its standard deviation that is different from the standard deviation of the normal distribution itself. We can collect more samples and compute the mean for each of them and then we compute the distribution of the mean. Usually we use a small number of samples of the population so that the mean doesn't follow exactly the normal distribution and we better use the Student distribution that works better with sample of small size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}