1DMixture_numpy_random.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1-D mixture of distributions using numpy.random"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This code example shows how to generate samples from a mixture of distributions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define distributions\n",
    "Here the mixture is composed of two normal, $\\mathcal{N}(-3, 2^2), \\mathcal{N}(5, 2.5^2)$, and one uniform, $\\mathcal{U}[1,3]$ distribution. They are stored in a list, with each distribution stored as a dictionary. The entries in the dictionary give the function name of the distirbution and the parameters. The parameters are themselves stored as a dictionary. https://docs.python.org/2/library/stdtypes.html#mapping-types-dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "distributions = [\n",
    "    dict(type= np.random.normal, kwargs=dict(loc = -3, scale = 2)),\n",
    "    dict(type= np.random.uniform, kwargs=dict(low = 1, high = 3)),\n",
    "    dict(type= np.random.normal, kwargs=dict(loc = 5, scale = 2.5)),\n",
    "]\n",
    "num_distr = len(distributions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define mixing coefficients\n",
    "If necessary, the coefficients are scaled to add to one. An exception is raised if the number of coefficients is not equal to the number of distributuions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "coefficients = np.array([0.5, 0.2, 0.3])\n",
    "coefficients /= coefficients.sum()\n",
    "num_coeff = len(coefficients)\n",
    "if num_coeff != num_distr:\n",
    "    raise Exception('number of mixing coefficients has to equal number of distributions')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate samples\n",
    "\n",
    "### Define sample size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_size = 10000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Allocate storage space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.zeros((sample_size, num_distr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sample the distribution which generates each sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "mixing_idx = np.random.choice(len(coefficients), size=(sample_size,), replace=True, p=coefficients)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate enough samples from each distribution\n",
    "The code loops over the distributions. The function name is accessed via the \"type\" entry and the parameters via the \"kwargs\" entry. The function is called using keyword arguments (https://docs.python.org/3/tutorial/controlflow.html#keyword-arguments). The size is one additional argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "for idx, distr in enumerate(distributions):\n",
    "    data[:, idx] = distr[\"type\"](size=(sample_size,), **distr[\"kwargs\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Choose the samples generated by the specified distributions\n",
    "np.arange gives the indices 0 to sample_size-1, equivalent to Matlab 0:samplesize-1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample = data[np.arange(sample_size), mixing_idx]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot histogram of samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAASkUlEQVR4nO3db4xc133e8e8TqlJap5EdiXYdUiwZiHlBp4FRr6kWbdw0rmUqQcQ4JSvKAUKhApgiJfIibRAabWWVzgvLaKoWNVuYrdQoNhxaEWCUqJkwjgUkQKA4pBxXNqWw2TCquKIR0aGqQmkVhdavL/bSGo1nuXe5szu7Z74fYMH759y7Z2aXz/zm3DtnU1VIktr1bZPugCRpZRn0ktQ4g16SGmfQS1LjDHpJatx1k+7AsJtvvrm2bt066W5I0rry5JNPfr2qNo7at+aCfuvWrZw+fXrS3ZCkdSXJ/1pon0M3ktS4XkGfZFeSs0lmkxwasf89Sb6U5HKSPSP2f2eS55N8fBydliT1t2jQJ9kAHAHuAHYAdyfZMdTsOeAe4NMLnOYjwG9dezclSdeqT0W/E5itqnNV9SpwDNg92KCqnq2qp4DXhg9O8i7gbcBvjKG/kqQl6hP0m4DzA+tz3bZFJfk24BeBn1t61yRJ49An6DNiW9+Z0H4aOFFV56/WKMmBJKeTnL548WLPU0uS+uhze+UccMvA+mbgQs/z/23gB5L8NPAdwPVJXq6qN1zQraqjwFGAmZkZp9OUpDHqE/SngO1JtgHPA/uAD/Y5eVX9xJXlJPcAM8MhL0laWYsO3VTVZeAgcBJ4Bni0qs4kOZzkToAk704yB+wFPpHkzEp2WpLUX9baHx6ZmZkpPxmr9Wbroc99c/nZj/7IBHuiaZXkyaqaGbVvzU2BIK0mA1rTwCkQJKlxBr0kNc6hG2kEh3TUEit6SWqcFb10jQarfmkts6KXpMYZ9JLUOINekhpn0EtS47wYK3W8uKpWWdFLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxTIEhL4DQJWo96VfRJdiU5m2Q2yaER+9+T5EtJLifZM7D9nUmeSHImyVNJ7hpn5yVJi1s06JNsAI4AdwA7gLuT7Bhq9hxwD/Dpoe3/F/jJqnoHsAv4d0nevNxOS5L66zN0sxOYrapzAEmOAbuBp680qKpnu32vDR5YVf9zYPlCkheAjcD/XnbPJUm99Bm62QScH1if67YtSZKdwPXAH43YdyDJ6SSnL168uNRTS5Kuok9FnxHbainfJMnbgU8C+6vqteH9VXUUOAowMzOzpHNLK80LsFrv+lT0c8AtA+ubgQt9v0GS7wQ+B/zLqvrdpXVPkrRcfYL+FLA9ybYk1wP7gON9Tt61/yzwy1X1q9feTUnStVo06KvqMnAQOAk8AzxaVWeSHE5yJ0CSdyeZA/YCn0hypjv8HwHvAe5J8uXu650r8kgkSSP1+sBUVZ0ATgxtu29g+RTzQzrDx30K+NQy+yhJWgY/GSuN2eDF22c/+iMT7Ik0z7luJKlxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxjnXjb7JOVqkNlnRS1LjrOg1dfzTgJo2VvSS1DiDXpIa59CNRvLCrNQOg16LMvSl9c2hG0lqnBX9lPMOFKl9vSr6JLuSnE0ym+TQiP3vSfKlJJeT7Bnatz/JH3Zf+8fVcUlSP4sGfZINwBHgDmAHcHeSHUPNngPuAT49dOx3AR8GbgN2Ah9O8pbld1uS1Fefin4nMFtV56rqVeAYsHuwQVU9W1VPAa8NHft+4PNVdamqXgQ+D+waQ78lST31CfpNwPmB9bluWx+9jk1yIMnpJKcvXrzY89SSpD76BH1GbKue5+91bFUdraqZqprZuHFjz1NLkvroE/RzwC0D65uBCz3Pv5xjJUlj0CfoTwHbk2xLcj2wDzje8/wngduTvKW7CHt7t02StEoWDfqqugwcZD6gnwEeraozSQ4nuRMgybuTzAF7gU8kOdMdewn4CPMvFqeAw902SdIq6fWBqao6AZwY2nbfwPIp5odlRh37MPDwMvooSVoGp0CQpMY5BULDVmIysoXO6cRn0tpl0E8h57eRpotDN5LUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc776Buz0D3y3jsvTS8reklqnBW9rpnvEqT1wYpekhpn0EtS4wx6SWqcQS9JjfNirMbOuemltcWKXpIaZ9BLUuMMeklqXK+gT7Irydkks0kOjdh/Q5LPdPu/mGRrt/0vJXkkyVeSPJPkQ+PtviRpMYsGfZINwBHgDmAHcHeSHUPN7gVerKpbgQeBB7rte4EbqupvAO8CfurKi4AkaXX0qeh3ArNVda6qXgWOAbuH2uwGHumWHwPemyRAAW9Kch3wl4FXgf8zlp5LknrpE/SbgPMD63PdtpFtquoy8BJwE/Oh/2fA14DngH9TVZeGv0GSA0lOJzl98eLFJT8ISdLC+gR9Rmyrnm12At8AvhvYBvyzJN/zLQ2rjlbVTFXNbNy4sUeXJEl99Qn6OeCWgfXNwIWF2nTDNDcCl4APAr9eVX9RVS8AvwPMLLfTkqT++gT9KWB7km1Jrgf2AceH2hwH9nfLe4DHq6qYH675ocx7E/C3gD8YT9clSX0sGvTdmPtB4CTwDPBoVZ1JcjjJnV2zh4CbkswCPwtcuQXzCPAdwFeZf8H4r1X11JgfgyTpKnrNdVNVJ4ATQ9vuG1h+hflbKYePe3nUdknS6nFSM60aJzuTJsMpECSpcQa9JDXOoJekxjlG34DBsW9JGmZFL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhrn7ZVaUd76KU2eFb0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcX5gShPh34+VVk+vij7JriRnk8wmOTRi/w1JPtPt/2KSrQP7vj/JE0nOJPlKkm8fX/clSYtZNOiTbACOAHcAO4C7k+wYanYv8GJV3Qo8CDzQHXsd8Cngn1TVO4AfBP5ibL2XJC2qT0W/E5itqnNV9SpwDNg91GY38Ei3/Bjw3iQBbgeeqqr/AVBVf1pV3xhP1yVJffQZo98EnB9YnwNuW6hNVV1O8hJwE/C9QCU5CWwEjlXVx4a/QZIDwAGALVu2LPUxqCGO3Uvj16eiz4ht1bPNdcDfBX6i+/cDSd77LQ2rjlbVTFXNbNy4sUeXJEl99ano54BbBtY3AxcWaDPXjcvfCFzqtv9WVX0dIMkJ4G8CX1hmv6eaU/9KWoo+Ff0pYHuSbUmuB/YBx4faHAf2d8t7gMerqoCTwPcn+SvdC8DfA54eT9clSX0sWtF3Y+4HmQ/tDcDDVXUmyWHgdFUdBx4CPplklvlKfl937ItJ/i3zLxYFnKgqy1G9ge9QpJXV6wNTVXUCODG07b6B5VeAvQsc+ynmb7GUJE2An4xdJ6x6JV0rg15aB7ztVMvhpGaS1DiDXpIa59DNGua4vKRxsKKXpMYZ9JLUOINekhpn0EtS4wx6SWqcd91IU8wPYk0HK3pJapwV/RrjvfNaioUq8uHfI6v16WbQSxPgkIlWk0EvrTMr9a7PF592OUYvSY2zopemjNeBpo8VvSQ1zopeWqPGWXlbxU83K3pJapwVvdYs7wJZGqt2LaRXRZ9kV5KzSWaTHBqx/4Ykn+n2fzHJ1qH9W5K8nOSfj6fbkqS+Fq3ok2wAjgDvA+aAU0mOV9XTA83uBV6sqluT7AMeAO4a2P8g8Gvj67ak1eI7q/Wvz9DNTmC2qs4BJDkG7AYGg343cH+3/Bjw8SSpqkryY8A54M/G1mtpHXJoRZPSJ+g3AecH1ueA2xZqU1WXk7wE3JTk/wE/z/y7gQWHbZIcAA4AbNmypXfnpRZYMWul9Qn6jNhWPdv8a+DBqno5GdWka1h1FDgKMDMzM3zu5lnpSVpJfYJ+DrhlYH0zcGGBNnNJrgNuBC4xX/nvSfIx4M3Aa0leqaqPL7vnmirrdZbGpb6I+6KvldAn6E8B25NsA54H9gEfHGpzHNgPPAHsAR6vqgJ+4EqDJPcDLxvykrS6Fg36bsz9IHAS2AA8XFVnkhwGTlfVceAh4JNJZpmv5PetZKclrSzfWbSl1wemquoEcGJo230Dy68Aexc5x/3X0D9J0jI5BYIkNc4pECQtm7eIrm0GvZqy0Niy4aNp5tCNJDXOil7rjneESEtjRS9JjTPoJalxDt1I6s27a9YnK3pJapxBL0mNM+glqXEGvSQ1zouxE+K94JJWixW9JDXOil7SNfFd6fphRS9JjTPoJalxDt2sIt/qSpoEg36FGe6SJs2hG0lqnBW9pLHyr3ytPb0q+iS7kpxNMpvk0Ij9NyT5TLf/i0m2dtvfl+TJJF/p/v2h8XZfkrSYRSv6JBuAI8D7gDngVJLjVfX0QLN7gRer6tYk+4AHgLuArwM/WlUXknwfcBLYNO4HIS3GayWaZn0q+p3AbFWdq6pXgWPA7qE2u4FHuuXHgPcmSVX9flVd6LafAb49yQ3j6LgkqZ8+Y/SbgPMD63PAbQu1qarLSV4CbmK+or/iHwK/X1V/fu3dlbRe+UdLJqdP0GfEtlpKmyTvYH445/aR3yA5ABwA2LJlS48uSZL66jN0MwfcMrC+GbiwUJsk1wE3Ape69c3AZ4GfrKo/GvUNqupoVc1U1czGjRuX9ggkSVfVp6I/BWxPsg14HtgHfHCozXFgP/AEsAd4vKoqyZuBzwEfqqrfGV+3Ja1nfS6OO7wzPotW9FV1GTjI/B0zzwCPVtWZJIeT3Nk1ewi4Kcks8LPAlVswDwK3Av8qyZe7r7eO/VFIkhbU6wNTVXUCODG07b6B5VeAvSOO+wXgF5bZx3XB2/ckrVV+MnaJvHNA0npj0EtakyyqxsdJzSSpcVb0ktY8q/vlsaKXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalx3ke/DM5vI2k9sKKXpMZZ0S/AT+JJaoUVvSQ1zope0rqy0LUx33kvzKCX1ByHXt9o6oO+zy+Ed9dIa5//Txc29UEvqW1W91MU9P6wJU2rpoN+qW/lfOsnTZdpubDbXNAb1pIWci350MJoQK+gT7IL+PfABuC/VNVHh/bfAPwy8C7gT4G7qurZbt+HgHuBbwA/U1Unx9b7a+SLgaRpsmjQJ9kAHAHeB8wBp5Icr6qnB5rdC7xYVbcm2Qc8ANyVZAewD3gH8N3Abyb53qr6xrgfiCSNy0LFYJ8icbDqX+pdfSv1jqFPRb8TmK2qcwBJjgG7gcGg3w3c3y0/Bnw8Sbrtx6rqz4E/TjLbne+J8XRfktaWtThi0CfoNwHnB9bngNsWalNVl5O8BNzUbf/doWM3DX+DJAeAA93qy0nODjW5Gfh6j75OA5+L1/lcvJHPx+vW3HORB8bT5ir++kI7+gR9Rmyrnm36HEtVHQWOLtiB5HRVzVytk9PC5+J1Phdv5PPxOp+LN+ozqdkccMvA+mbgwkJtklwH3Ahc6nmsJGkF9Qn6U8D2JNuSXM/8xdXjQ22OA/u75T3A41VV3fZ9SW5Isg3YDvzeeLouSepj0aGbbsz9IHCS+dsrH66qM0kOA6er6jjwEPDJ7mLrJeZfDOjaPcr8hdvLwD+9xjtuFhzWmUI+F6/zuXgjn4/X+VwMyHzhLUlqlX94RJIaZ9BLUuPWdNAn2ZvkTJLXkswM7ftQktkkZ5O8f1J9nIQk9yd5PsmXu68fnnSfVluSXd3PfjbJoUn3Z5KSPJvkK93vwulJ92e1JXk4yQtJvjqw7buSfD7JH3b/vmWSfZy0NR30wFeBHwd+e3Dj0NQKu4D/2E3VME0erKp3dl8nJt2Z1TQwLccdwA7g7u53Ypr9/e53YRrvHf8l5nNg0CHgC1W1HfhCtz611nTQV9UzVTX8KVkYmFqhqv4YuDK1gqbDN6flqKpXgSvTcmgKVdVvM3+336DdwCPd8iPAj61qp9aYNR30VzFqWoZvmVqhcQeTPNW9bZ22t6X+/N+ogN9I8mQ3nYjgbVX1NYDu37dOuD8TNfH56JP8JvDXRuz6F1X13xY6bMS2pu4TvdrzAvwn4CPMP+aPAL8I/OPV693ENf/zX6K/U1UXkrwV+HySP+iqXAlYA0FfVf/gGg5rfmqFvs9Lkv8M/PcV7s5a0/zPfymq6kL37wtJPsv80Na0B/2fJHl7VX0tyduBFybdoUlar0M3Uz21QveLe8UHmL9oPU36TMsxFZK8KclfvbIM3M70/T6MMjgty35godGBqTDxiv5qknwA+A/ARuBzSb5cVe8f49QK69XHkryT+eGKZ4Gfmmx3VtdC03JMuFuT8jbgs/N//oHrgE9X1a9PtkurK8mvAD8I3JxkDvgw8FHg0ST3As8BeyfXw8lzCgRJatx6HbqRJPVk0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TG/X98sm3sJYmF5QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(sample, bins=100, density=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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": 2
}