notes/ImportanceSampling.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "from scipy.stats import uniform,beta,norm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Monte Carlo Sampling** method is to sample $x$ from the $p(x)$ and take the average of all samples to estimate the expectation.\n",
    "$$\\mathbb{E}[f(x)]=\\int f(x)p(x)dx \\approx \\frac{1}{n}\\sum_{i=1}^n f(x_i) \\qquad x_i\\sim p$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAPE0lEQVR4nO3dYYzkdX3H8fenHAShEkCWy5WTLiYXKmki2AliSUzKicFC4B6ogbTmYkiuD6zR2sSefdDEpA/OpKn6oDG5gLpNFcSr5Ig2VnJKWpMW3QOs4GFQeuLJebdaqKJNEf32wf4vXPdmd2dnZ3b2t/t+JZeZ/29nbj6TYz/89rf/339SVUiS2vMbkw4gSRqOBS5JjbLAJalRFrgkNcoCl6RGbVnLF7vkkktqenp6LV9Skpp3+PDhH1fV1MLxZQs8yZXAZ08beg3wV8Dfd+PTwFHgHVX13FJ/1/T0NLOzs4OnliSR5Pv9xpddQqmq71TV1VV1NfB7wC+A+4G9wKGq2gEc6o4lSWtkpWvgO4HvVdX3gduAmW58Btg1ymCSpKWttMBvB+7p7m+tquMA3e2lowwmSVrawAWe5BzgVuBzK3mBJHuSzCaZnZubW2k+SdIiVjIDfyvwSFWd6I5PJNkG0N2e7PekqtpfVb2q6k1NnfFLVEnSkFZS4Hfw8vIJwAPA7u7+buDgqEJJkpY3UIEnOQ+4Efj8acP7gBuTPNV9bd/o40mSFjPQRp6q+gXwqgVjP2H+rBRJ0gS4lV6SGrWmW+mltTK994t9x4/uu3mNk0jj4wxckhplgUtSoyxwSWqUBS5JjbLAJalRFrgkNcoCl6RGWeCS1Cg38qhpi23YkTYDZ+CS1CgLXJIaZYFLUqMscElqlAUuSY3yLBQ1wbNNpDM5A5ekRlngktQoC1ySGmWBS1KjLHBJatRABZ7kwiQHkjyZ5EiSNya5OMmDSZ7qbi8ad1hJ0ssGnYF/DPhSVf0O8DrgCLAXOFRVO4BD3bEkaY0sW+BJLgDeBNwNUFUvVtXzwG3ATPewGWDXuEJKks40yAz8NcAc8Mkkjya5K8n5wNaqOg7Q3V7a78lJ9iSZTTI7Nzc3suCStNkNUuBbgNcDH6+qa4Cfs4LlkqraX1W9qupNTU0NGVOStNAgBX4MOFZVD3fHB5gv9BNJtgF0tyfHE1GS1M+yBV5VPwJ+kOTKbmgn8G3gAWB3N7YbODiWhJKkvga9mNV7gE8nOQd4GngX8+V/X5I7gWeAt48noiSpn4EKvKoeA3p9vrRztHEkSYNyJ6YkNcoCl6RGWeCS1CgLXJIaZYFLUqMscElqlAUuSY2ywCWpURa4JDVq0K300oYwvfeLfceP7rt5jZNIq+cMXJIaZYFLUqMscElqlAUuSY2ywCWpURa4JDXKApekRlngktQoC1ySGmWBS1KjLHBJatRA10JJchT4GfAr4KWq6iW5GPgsMA0cBd5RVc+NJ6YkaaGVzMD/oKqurqped7wXOFRVO4BD3bEkaY2sZgnlNmCmuz8D7Fp9HEnSoAYt8AK+nORwkj3d2NaqOg7Q3V46joCSpP4GvR749VX1bJJLgQeTPDnoC3SFvwfg8ssvHyKiJKmfgWbgVfVsd3sSuB+4FjiRZBtAd3tykefur6peVfWmpqZGk1qStHyBJzk/yStP3QfeAjwOPADs7h62Gzg4rpCSpDMNsoSyFbg/yanHf6aqvpTkG8B9Se4EngHePr6Y2mj8aDNp9ZYt8Kp6Gnhdn/GfADvHEUqStDx3YkpSoyxwSWqUBS5JjbLAJalRg27kkdbEYmenSDqTM3BJapQFLkmNssAlqVEWuCQ1ygKXpEZ5ForGyrNKpPFxBi5JjbLAJalRFrgkNcoCl6RGWeCS1CjPQtHAljqjxE/SkdaeM3BJapQFLkmNcglFWoIfvqz1zBm4JDXKApekRg28hJLkLGAW+GFV3ZLkCuBe4GLgEeCdVfXieGJK4+U1W9SilczA3wscOe34w8BHqmoH8Bxw5yiDSZKWNlCBJ9kO3Azc1R0HuAE40D1kBtg1joCSpP4GnYF/FPgA8Ovu+FXA81X1Und8DLis3xOT7Ekym2R2bm5uVWElSS9btsCT3AKcrKrDpw/3eWj1e35V7a+qXlX1pqamhowpSVpokF9iXg/cmuQPgXOBC5ifkV+YZEs3C98OPDu+mJKkhZadgVfVB6tqe1VNA7cDX6mqPwK+Cryte9hu4ODYUkqSzrCa88D/Anh/ku8yvyZ+92giSZIGsaKt9FX1EPBQd/9p4NrRR5IkDcKdmJLUKAtckhplgUtSoyxwSWqUBS5JjbLAJalRFrgkNcoCl6RGWeCS1CgLXJIaZYFLUqMscElq1IouZiUtxg8FltaeM3BJapQFLkmNcglFGqFhlpKO7rt5DEm0GTgDl6RGWeCS1CgLXJIaZYFLUqMscElq1LIFnuTcJF9P8s0kTyT5UDd+RZKHkzyV5LNJzhl/XEnSKYPMwP8XuKGqXgdcDdyU5Drgw8BHqmoH8Bxw5/hiSpIWWrbAa94L3eHZ3Z8CbgAOdOMzwK6xJJQk9TXQGniSs5I8BpwEHgS+BzxfVS91DzkGXLbIc/ckmU0yOzc3N4rMkiQGLPCq+lVVXQ1sB64FXtvvYYs8d39V9aqqNzU1NXxSSdL/s6KzUKrqeeAh4DrgwiSntuJvB54dbTRJ0lKWvRZKkingl1X1fJJXAG9m/heYXwXeBtwL7AYOjjOo1o6XhpXaMMjFrLYBM0nOYn7Gfl9VfSHJt4F7k/w18Chw9xhzSpIWWLbAq+o/gGv6jD/N/Hq4JGkC3IkpSY2ywCWpURa4JDXKApekRlngktQoC1ySGmWBS1Kj/FT6Tcwdl1LbnIFLUqMscElqlEsoG8hiSyJH9928xkk2PpeftB44A5ekRlngktQoC1ySGmWBS1KjLHBJalSq+n4W8Vj0er2anZ1ds9fbqDwDYnPw7CGdkuRwVfUWjjsDl6RGWeCS1Cg38kgbnBu8Ni5n4JLUKAtckhq1bIEneXWSryY5kuSJJO/txi9O8mCSp7rbi8YfV5J0yiAz8JeAP6+q1wLXAe9OchWwFzhUVTuAQ92xJGmNLFvgVXW8qh7p7v8MOAJcBtwGzHQPmwF2jSukJOlMKzoLJck0cA3wMLC1qo7DfMknuXSR5+wB9gBcfvnlq8kqCc8q0csG/iVmkt8E/hF4X1X9dNDnVdX+qupVVW9qamqYjJKkPgYq8CRnM1/en66qz3fDJ5Js676+DTg5noiSpH6WXUJJEuBu4EhV/e1pX3oA2A3s624PjiWhtEl5zRstZ5A18OuBdwLfSvJYN/aXzBf3fUnuBJ4B3j6eiJKkfpYt8Kr6GpBFvrxztHEkSYPyWijrgGcVSBqGW+klqVEWuCQ1ygKXpEZZ4JLUKAtckhrlWSjSBuHGn83HGbgkNcoCl6RGWeCS1CgLXJIaZYFLUqM8C0XapLwGT/ucgUtSoyxwSWqUSyjrmBszJC3FGbgkNcoCl6RGuYSyCv4WXxuR/123wxm4JDXKApekRi1b4Ek+keRkksdPG7s4yYNJnupuLxpvTEnSQoPMwD8F3LRgbC9wqKp2AIe6Y0nSGlq2wKvqX4D/WjB8GzDT3Z8Bdo04lyRpGcOugW+tquMA3e2liz0wyZ4ks0lm5+bmhnw5SdJCY/8lZlXtr6peVfWmpqbG/XKStGkMW+AnkmwD6G5Pji6SJGkQwxb4A8Du7v5u4OBo4kiSBjXIaYT3AP8GXJnkWJI7gX3AjUmeAm7sjiVJa2jZrfRVdcciX9o54iySpBXwWiiSBuI1UtYft9JLUqMscElqlEsoY+An6WgzcWllcpyBS1KjLHBJapQFLkmNssAlqVEWuCQ1ygKXpEZZ4JLUKAtckhplgUtSo9yJKWndc7dnf87AJalRFrgkNaqZJZSV/gjlBaWkyRrme3CzL4mslDNwSWqUBS5JjUpVrdmL9Xq9mp2dHeq5LolIGtRGW4pJcriqegvHnYFLUqMscElq1KrOQklyE/Ax4CzgrqraN5JUkrQOjOrst3Et6Qw9A09yFvB3wFuBq4A7klw1qmCSpKWtZgnlWuC7VfV0Vb0I3AvcNppYkqTlrGYJ5TLgB6cdHwPesPBBSfYAe7rDF5J8ZxWvOWmXAD+edIgR8H2sLxvlfcA6eS/58Kr/iiXfx0r//hHk+e1+g6sp8PQZO+OcxKraD+xfxeusG0lm+53K0xrfx/qyUd4HbJz30sr7WM0SyjHg1acdbweeXV0cSdKgVlPg3wB2JLkiyTnA7cADo4klSVrO0EsoVfVSkj8F/pn50wg/UVVPjCzZ+rQhloLwfaw3G+V9wMZ5L028jzXdSi9JGh13YkpSoyxwSWqUBT6AJOcm+XqSbyZ5IsmHJp1pNZKcleTRJF+YdJZhJTma5FtJHksy3CUu14EkFyY5kOTJJEeSvHHSmVYqyZXdv8OpPz9N8r5J5xpGkj/rvscfT3JPknMnnWkproEPIEmA86vqhSRnA18D3ltV/z7haENJ8n6gB1xQVbdMOs8wkhwFelU18U0jq5FkBvjXqrqrO5vrvKp6ftK5htVdYuOHwBuq6vuTzrMSSS5j/nv7qqr6nyT3Af9UVZ+abLLFOQMfQM17oTs8u/vT5P/5kmwHbgbumnSWzS7JBcCbgLsBqurFlsu7sxP4XmvlfZotwCuSbAHOY53vbbHAB9QtOzwGnAQerKqHJ51pSB8FPgD8etJBVqmALyc53F2uoUWvAeaAT3ZLWnclOX/SoVbpduCeSYcYRlX9EPgb4BngOPDfVfXlyaZamgU+oKr6VVVdzfyO02uT/O6kM61UkluAk1V1eNJZRuD6qno981fDfHeSN0060BC2AK8HPl5V1wA/B/ZONtLwuiWgW4HPTTrLMJJcxPwF+a4Afgs4P8kfTzbV0izwFep+xH0IuGnCUYZxPXBrt358L3BDkn+YbKThVNWz3e1J4H7mr47ZmmPAsdN+mjvAfKG36q3AI1V1YtJBhvRm4D+raq6qfgl8Hvj9CWdakgU+gCRTSS7s7r+C+X/oJyebauWq6oNVtb2qppn/UfcrVbWuZxj9JDk/yStP3QfeAjw+2VQrV1U/An6Q5MpuaCfw7QlGWq07aHT5pPMMcF2S87oTF3YCRyacaUmr+kSeTWQbMNP9hv03gPuqqtlT8DaArcD9899jbAE+U1Vfmmykob0H+HS3/PA08K4J5xlKkvOAG4E/mXSWYVXVw0kOAI8ALwGPss631HsaoSQ1yiUUSWqUBS5JjbLAJalRFrgkNcoCl6RGWeCS1CgLXJIa9X92ENotdMHiqgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mu,sig,N = 5,1,1000\n",
    "X = sig*np.random.randn(N)+mu\n",
    "_ = plt.hist(X,bins=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.045865815217911 0.9028513930783005\n"
     ]
    }
   ],
   "source": [
    "mu_ = sum(X)/N\n",
    "sig_ = sum((X-mu_)**2)/N\n",
    "print(mu_,sig_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 0.9049550660839032\n"
     ]
    }
   ],
   "source": [
    "f = lambda x: (x-mu)**2\n",
    "target = sig\n",
    "ans = sum(f(X))/N\n",
    "print(target,ans)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Importance Sampling**\n",
    "If $p(x)$ is hard to sample from, we can sample from easy distribution $q(x)$ and calculate the sampling ratio as correction weight to offset the probability sampling from a different distribution.\n",
    "$$\\mathbb{E}[f(x)]=\\int f(x)p(x)dx = \\int f(x)\\frac{p(x)}{q(x)}q(x)dx \\approx \\frac{1}{n}\\sum^n_{i=1}f(x_i)\\frac{p(x_i)}{q(x_i)} \\qquad x_i\\sim q$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.591056097376715\n"
     ]
    }
   ],
   "source": [
    "xmin,xmax = 0,9\n",
    "X = np.random.uniform(xmin,xmax,N)\n",
    "mu_is = sum(X*norm.pdf(X,mu,sig)/uniform.pdf(X,xmin,xmax))/N\n",
    "print(mu_is)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.150207478521478\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAPOklEQVR4nO3dbaxlVX3H8e9PHoqPQcuRTGe4HbTEakwczA2lncRY0AbFCCaaSFqdNCTXJtpiS2rRN2rTJpio2BfGdBRkklKUIAZCrZUgxpg0Y2dg1MHRYJHqyMhohAp9oR3898Xdk453zp1z7rnn4a57vp/k5p69zt5n/w8DP9bstddeqSokSe15xqwLkCSNxgCXpEYZ4JLUKANckhplgEtSo06f5snOOeec2r59+zRPKUnN279//0+rqreyfaoBvn37dvbt2zfNU0pS85L8V792L6FIUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRg0d4ElOS/JAkru77fOT7E3yUJLPJjlzcmVKklZaSw/8GuDQCdsfAm6oqguAx4Grx1mYJOnUhgrwJNuAy4FPddsBLgFu73bZA1w5iQIlSf0NOxPzY8B7gOd2278JPFFVx7rtw8DWfgcmWQKWABYWFkavVNpgtl/3L33bH7n+8ilXonk1sAee5A3A0araf2Jzn137Lu1TVburarGqFnu9k6byS5JGNEwPfCfwxiSvB84Cnsdyj/zsJKd3vfBtwKOTK1OStNLAHnhVvbeqtlXVduCtwJer6o+B+4A3d7vtAu6cWJWSpJOs5z7wvwH+Ksn3WL4mfuN4SpIkDWNNj5Otqq8AX+lePwxcNP6SNM8cGJSG50xMSWqUAS5JjTLAJalRBrgkNWqqa2JKo3JwUzqZPXBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRwyxqfFaSryf5RpIHk3ywa785yfeTHOh+dky+XEnSccM8zOoXwCVV9VSSM4CvJfnX7r2/rqrbJ1eeJGk1AwO8qgp4qts8o/upSRYlSRpsqGvgSU5LcgA4CtxTVXu7t/4+yTeT3JDkNyZWpSTpJEM9D7yqngZ2JDkb+HySlwPvBX4MnAnsZnmV+r9deWySJWAJYGFhYUxlay2m8Sxtn9ctTd+a7kKpqidYXpX+sqo6Ust+AXyaVVaor6rdVbVYVYu9Xm/dBUuSlg1zF0qv63mT5JnAa4DvJNnStQW4Ejg4yUIlSb9umEsoW4A9SU5jOfBvq6q7k3w5SQ8IcAD4swnWKUlaYZi7UL4JXNin/ZKJVCRJGoqLGm8iqw0kzvIcDm5Kk+NUeklqlAEuSY0ywCWpUQa4JDXKQUzNFQdVtZnYA5ekRhngktQoA1ySGmWAS1KjHMSUcHBTbbIHLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckho1zJqYZyX5epJvJHkwyQe79vOT7E3yUJLPJjlz8uVKko4bpgf+C+CSqnoFsAO4LMnFwIeAG6rqAuBx4OrJlSlJWmlggNeyp7rNM7qfAi4Bbu/a97C8Mr0kaUqGmonZrUi/H/gd4OPAfwJPVNWxbpfDwNZVjl0ClgAWFhbWW6/0ayY9g3Ia64xKoxpqELOqnq6qHcA24CLgpf12W+XY3VW1WFWLvV5v9EolSb9mTXehVNUTwFeAi4GzkxzvwW8DHh1vaZKkUxnmLpRekrO7188EXgMcAu4D3tzttgu4c1JFSpJONsw18C3Anu46+DOA26rq7iTfBj6T5O+AB4AbJ1inJGmFgQFeVd8ELuzT/jDL18MlSTPg88DnmHdYSG1zKr0kNcoAl6RGGeCS1CgDXJIalaq+EygnYnFxsfbt2ze1880bByU3NhdI1qiS7K+qxZXt9sAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGDbMm5nlJ7ktyKMmDSa7p2j+Q5EdJDnQ/r598uZKk44ZZkecYcG1V3Z/kucD+JPd0791QVR+eXHmSpNUMsybmEeBI9/rJJIeArZMuTJJ0amtaEzPJdpYXON4L7ATeleTtwD6We+mP9zlmCVgCWFhYWGe50uaz2mOAffysBhl6EDPJc4DPAe+uqp8DnwBeDOxguYf+kX7HVdXuqlqsqsVerzeGkiVJMGSAJzmD5fC+paruAKiqx6rq6ar6FfBJ4KLJlSlJWmmYu1AC3AgcqqqPntC+5YTd3gQcHH95kqTVDHMNfCfwNuBbSQ50be8DrkqyAyjgEeAdE6lQktTXMHehfA1In7e+MP5yJEnDWtNdKBrOWhcX9m4DSaNwKr0kNcoAl6RGGeCS1CgDXJIa5SCmNCVrHdyWBrEHLkmNMsAlqVEGuCQ1ygCXpEY5iNkgB8Pmg88J1yD2wCWpUQa4JDXKAJekRhngktQoBzHXwcFEzYKDmzrOHrgkNWqYNTHPS3JfkkNJHkxyTdf+giT3JHmo+/38yZcrSTpumB74MeDaqnopcDHwziQvA64D7q2qC4B7u21J0pQMDPCqOlJV93evnwQOAVuBK4A93W57gCsnVaQk6WRrGsRMsh24ENgLnFtVR2A55JO8cJVjloAlgIWFhfXUOnccJJV0KkMPYiZ5DvA54N1V9fNhj6uq3VW1WFWLvV5vlBolSX0MFeBJzmA5vG+pqju65seSbOne3wIcnUyJkqR+hrkLJcCNwKGq+ugJb90F7Ope7wLuHH95kqTVDHMNfCfwNuBbSQ50be8DrgduS3I18APgLZMpUZLUz8AAr6qvAVnl7UvHW87G5GCipI3ImZiS1CgDXJIaZYBLUqMMcElqlAEuSY3yeeAbgHe5SBqFPXBJapQBLkmNMsAlqVEGuCQ1ykFMaZNY62D4aosgj+tzNHn2wCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjhlkT86YkR5McPKHtA0l+lORA9/P6yZYpSVppmB74zcBlfdpvqKod3c8XxluWJGmQgQFeVV8FfjaFWiRJa7CemZjvSvJ2YB9wbVU93m+nJEvAEsDCwsI6TtffarPGRpkd5mNdJbVk1EHMTwAvBnYAR4CPrLZjVe2uqsWqWuz1eiOeTpK00kgBXlWPVdXTVfUr4JPAReMtS5I0yEgBnmTLCZtvAg6utq8kaTIGXgNPcivwauCcJIeB9wOvTrIDKOAR4B0TrFGS1MfAAK+qq/o03ziBWiRJa+BMTElqlAEuSY0ywCWpUQa4JDXKNTGlOeXM4/bZA5ekRhngktQoA1ySGmWAS1KjDHBJatTc3YXiyLs0Xqf6b2qU5/JrePbAJalRBrgkNcoAl6RGGeCS1Ki5G8SUND3jXHRcJ7MHLkmNGhjgSW5KcjTJwRPaXpDkniQPdb+fP9kyJUkrDdMDvxm4bEXbdcC9VXUBcG+3LUmaooEBXlVfBX62ovkKYE/3eg9w5ZjrkiQNMOog5rlVdQSgqo4keeFqOyZZApYAFhYWRjydpHngoOfaTHwQs6p2V9ViVS32er1Jn06S5saoAf5Yki0A3e+j4ytJkjSMUQP8LmBX93oXcOd4ypEkDWuY2whvBf4deEmSw0muBq4HXpvkIeC13bYkaYoGDmJW1VWrvHXpmGuRpL4c3OzPmZiS1CgDXJIaZYBLUqMMcElq1KZ9nKxrX0ra7OyBS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWrUpp1KL2l+rfX54a0+b9weuCQ1al098CSPAE8CTwPHqmpxHEVJkgYbxyWUP6yqn47hcyRJa+AlFElq1Hp74AV8KUkB/1hVu1fukGQJWAJYWFhY5+kkbQbjel7/vD/3f7098J1V9UrgdcA7k7xq5Q5VtbuqFqtqsdfrrfN0kqTj1hXgVfVo9/so8HngonEUJUkabOQAT/LsJM89/hr4I+DguAqTJJ3aeq6Bnwt8Psnxz/nnqvriWKqSJA00coBX1cPAK8ZYyynN+2CFpPXbbDnibYSS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ar2LGkvSprXa88Mfuf7yNe1/qmPWwx64JDVqXQGe5LIk303yvSTXjasoSdJg61nU+DTg48DrgJcBVyV52bgKkySd2np64BcB36uqh6vql8BngCvGU5YkaZD1DGJuBX54wvZh4PdW7pRkCVjqNp9K8t01nucc4KcjVdg+v/t88rtvcPnQRI451Xf/7X6N6wnw9GmrkxqqdgO7Rz5Jsq+qFkc9vmV+d7/7vPG7r+27r+cSymHgvBO2twGPruPzJElrsJ4A/w/ggiTnJzkTeCtw13jKkiQNMvIllKo6luRdwL8BpwE3VdWDY6vs/418+WUT8LvPJ7/7fFrzd0/VSZetJUkNcCamJDXKAJekRm3oAJ/XqfpJzktyX5JDSR5Mcs2sa5q2JKcleSDJ3bOuZZqSnJ3k9iTf6f78f3/WNU1Lkr/s/n0/mOTWJGfNuqZJSXJTkqNJDp7Q9oIk9yR5qPv9/EGfs2EDfM6n6h8Drq2qlwIXA++co+9+3DXAoVkXMQP/AHyxqn4XeAVz8s8gyVbgL4DFqno5yzdGvHW2VU3UzcBlK9quA+6tqguAe7vtU9qwAc4cT9WvqiNVdX/3+kmW/yPeOtuqpifJNuBy4FOzrmWakjwPeBVwI0BV/bKqnphtVVN1OvDMJKcDz2ITzyupqq8CP1vRfAWwp3u9B7hy0Ods5ADvN1V/bkLsuCTbgQuBvbOtZKo+BrwH+NWsC5myFwE/AT7dXT76VJJnz7qoaaiqHwEfBn4AHAH+u6q+NNuqpu7cqjoCy5044IWDDtjIAT7UVP3NLMlzgM8B766qn8+6nmlI8gbgaFXtn3UtM3A68ErgE1V1IfA/DPHX6M2gu957BXA+8FvAs5P8yWyr2vg2coDP9VT9JGewHN63VNUds65ninYCb0zyCMuXzS5J8k+zLWlqDgOHq+r437ZuZznQ58FrgO9X1U+q6n+BO4A/mHFN0/ZYki0A3e+jgw7YyAE+t1P1k4Tl66CHquqjs65nmqrqvVW1raq2s/xn/uWqmoueWFX9GPhhkpd0TZcC355hSdP0A+DiJM/q/v2/lDkZwD3BXcCu7vUu4M5BB2zYNTGnOFV/I9oJvA34VpIDXdv7quoLM6xJ0/HnwC1dp+Vh4E9nXM9UVNXeJLcD97N8F9YDbOJp9UluBV4NnJPkMPB+4HrgtiRXs/w/tLcM/Byn0ktSmzbyJRRJ0ikY4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalR/weZxuUc6LN2+wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "a,b = 2,2.7\n",
    "scale = mu*2\n",
    "X = np.random.beta(a,b,N)*scale\n",
    "_ = plt.hist(X,bins=50)\n",
    "mu_is = sum(X*norm.pdf(X,mu,sig)/beta.pdf(X,a,b,scale=scale))/N\n",
    "print(mu_is)"
   ]
  },
  {
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}