From eb26a94e0a20697bde057c51f4659f8a2ba30611 Mon Sep 17 00:00:00 2001 From: Alex Montiel Date: Sun, 12 Jan 2020 14:39:00 -0600 Subject: [PATCH] Lab bayesian statistics --- your-code/main.ipynb | 222 ++++++++++++++++++++++++++++++++++++++----- 1 file changed, 197 insertions(+), 25 deletions(-) diff --git a/your-code/main.ipynb b/your-code/main.ipynb index 95bfcb9..e79b040 100644 --- a/your-code/main.ipynb +++ b/your-code/main.ipynb @@ -11,13 +11,14 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", - "import matplotlib.pyplot as plt" + "import matplotlib.pyplot as plt\n", + "%matplotlib inline" ] }, { @@ -31,10 +32,39 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "import numpy as np\n", + "\n", + "def bayes_rule(priors, likelihoods):\n", + " marg = sum(np.multiply(priors, likelihoods))\n", + " post = np.divide(np.multiply(priors, likelihoods), marg)\n", + " return post" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0.6, 0.4])" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "priors = [1/2, 1/2]\n", + "likelihoods = [30/40, 20/40]\n", + "bayes_rule(priors, likelihoods)" + ] }, { "cell_type": "markdown", @@ -45,10 +75,13 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 9, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "#Para el BOWL 1 = 0.6\n", + "#Para el BOWL 2 = 0.4" + ] }, { "cell_type": "markdown", @@ -59,10 +92,25 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 10, "metadata": {}, - "outputs": [], - "source": [] + "outputs": [ + { + "data": { + "text/plain": [ + "array([0.33333333, 0.66666667])" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "chocoprior = [1/2, 1/2]\n", + "chocolike = [10/40, 20/40]\n", + "bayes_rule(chocoprior, chocolike)" + ] }, { "cell_type": "markdown", @@ -95,10 +143,47 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 17, "metadata": {}, - "outputs": [], - "source": [] + "outputs": [ + { + "data": { + "text/plain": [ + "array([0.58823529, 0.41176471])" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "bags = [1/2, 1/2]\n", + "yellow_candy = [0.20,0.14]\n", + "bayes_rule(bags, yellow_candy)" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0.33333333, 0.66666667])" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "bags = [1/2, 1/2]\n", + "green_candy = [0.10,0.20]\n", + "bayes_rule(bags,green_candy)" + ] }, { "cell_type": "markdown", @@ -112,7 +197,9 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "# Probability of 67%" + ] }, { "cell_type": "markdown", @@ -141,10 +228,27 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 23, "metadata": {}, - "outputs": [], - "source": [] + "outputs": [ + { + "data": { + "text/plain": [ + "array([0.33333333, 0.66666667, 0. ])" + ] + }, + "execution_count": 23, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#La probabilidad de elegir la puerta correcta es 1/3. \n", + "#0.33 por cada \n", + "doors = [1/3, 1/3, 1/3]\n", + "chance = [1/2,1,0]\n", + "bayes_rule(doors, chance)\n" + ] }, { "cell_type": "markdown", @@ -157,10 +261,15 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 11, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "new_users = 14\n", + "def modelo_generativo(parametro):\n", + " resultado = np.random.binomial(100, parametro)\n", + " return resultado" + ] }, { "cell_type": "markdown", @@ -171,10 +280,43 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 19, "metadata": {}, - "outputs": [], - "source": [] + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAP+UlEQVR4nO3df4wc9XnH8fdTSCjFaTByOLnGzRHJrWp6LUmuFBWpPYu2/FJjIoUKSsEGKqeVURPl/jFJpUSNkFBVJ2pUinoRCKOkcVCTCEuQVI7LCvGHG2zkchiXYswVjC2jBBc4iKjOPP3jxmIxd7693Z3b2++9X9JqZ787M/s8zPHx3PzYi8xEklSWX+h1AZKk7jPcJalAhrskFchwl6QCGe6SVKAze10AwIoVK3JwcLDXZfTEm2++yTnnnNPrMnrG/u1/qfbfjd737t3708z8yEzvLYpwHxwcZM+ePb0uoycajQYjIyO9LqNn7N/+l2r/3eg9Iv5ntvc8LCNJBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQVaFHeoqn8Mbnm4q+sbHZpiYwvrnLjrmq5+rlQ699wlqUCGuyQVyHCXpAIZ7pJUoDnDPSJWR8SjEXEgIvZHxOeq8a9ExMsRsa96XN20zB0RcTAino2IK+psQJL0fq1cLTMFjGbmkxHxIWBvROys3vt6Zv5988wRsRa4HrgI+BXgxxHxa5l5opuFS5JmN+eee2Yezcwnq+k3gAPAqtMssh7YnplvZ+YLwEHgkm4UK0lqTWRm6zNHDAKPAb8JfAHYCLwO7GF67/54RPwjsDszv1Utcy/ww8z811PWtQnYBDAwMPDJ7du3d9pLX5qcnGTZsmW9LqNl4y+/1tX1DZwNx34+93xDqz7c1c9dLPpt+3fbUu6/G72vW7dub2YOz/ReyzcxRcQy4HvA5zPz9Yi4B/gqkNXzVuBWIGZY/H3/gmTmGDAGMDw8nP6prf7Qyg1H8zE6NMXW8bl/DCduHOnq5y4W/bb9u20p91937y1dLRMRH2A62L+dmd8HyMxjmXkiM98Bvsm7h14OA6ubFr8AONK9kiVJc2nlapkA7gUOZObXmsZXNs32aeDpanoHcH1EnBURFwJrgJ90r2RJ0lxaOSxzGXATMB4R+6qxLwI3RMTFTB9ymQA+C5CZ+yPiQeAZpq+02eyVMpK0sOYM98x8nJmPoz9ymmXuBO7soC5JUge8Q1WSCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalALf8lJi0eg13+a0iSyuOeuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoH8Vkj1hV5+E+bEXdf07LOldrnnLkkFMtwlqUBzhntErI6IRyPiQETsj4jPVePnRcTOiHiuel5ejUdEfCMiDkbEUxHxibqbkCS9Vyt77lPAaGb+BnApsDki1gJbgF2ZuQbYVb0GuApYUz02Afd0vWpJ0mnNGe6ZeTQzn6ym3wAOAKuA9cC2arZtwLXV9HrggZy2Gzg3IlZ2vXJJ0qzmdcw9IgaBjwP/AQxk5lGY/gcAOL+abRXwUtNih6sxSdICaflSyIhYBnwP+Hxmvh4Rs846w1jOsL5NTB+2YWBggEaj0WopRZmcnJx376NDU/UU0wMDZy/+fur82Wxn+5dkKfdfd+8thXtEfIDpYP92Zn6/Gj4WESsz82h12OWVavwwsLpp8QuAI6euMzPHgDGA4eHhHBkZaa+DPtdoNJhv7xt7eM13t40OTbF1fHHfbjFx40ht625n+5dkKfdfd++tXC0TwL3Agcz8WtNbO4AN1fQG4KGm8Zurq2YuBV47efhGkrQwWtllugy4CRiPiH3V2BeBu4AHI+I24EXguuq9R4CrgYPAW8AtXa1YkjSnOcM9Mx9n5uPoAJfPMH8CmzusS5LUAe9QlaQCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQHOGe0TcFxGvRMTTTWNfiYiXI2Jf9bi66b07IuJgRDwbEVfUVbgkaXat7LnfD1w5w/jXM/Pi6vEIQESsBa4HLqqW+aeIOKNbxUqSWjNnuGfmY8CrLa5vPbA9M9/OzBeAg8AlHdQnSWrDmR0se3tE3AzsAUYz8ziwCtjdNM/haux9ImITsAlgYGCARqPRQSn9a3Jyct69jw5N1VNMDwycvfj7qfNns53tX5Kl3H/dvbcb7vcAXwWyet4K3ArEDPPmTCvIzDFgDGB4eDhHRkbaLKW/NRoN5tv7xi0P11NMD4wOTbF1vJN9jPpN3DhS27rb2f4lWcr91917W1fLZOaxzDyRme8A3+TdQy+HgdVNs14AHOmsREnSfLUV7hGxsunlp4GTV9LsAK6PiLMi4kJgDfCTzkqUJM3XnL8PR8R3gBFgRUQcBr4MjETExUwfcpkAPguQmfsj4kHgGWAK2JyZJ+opXZI0mznDPTNvmGH43tPMfydwZydFSZI64x2qklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXozF4XIC12g1serm3do0NTbJxl/RN3XVPb56p87rlLUoEMd0kqkOEuSQUy3CWpQJ5Q7UA3TrSd7oSaJLXLPXdJKpDhLkkFmjPcI+K+iHglIp5uGjsvInZGxHPV8/JqPCLiGxFxMCKeiohP1Fm8JGlmrey53w9cecrYFmBXZq4BdlWvAa4C1lSPTcA93SlTkjQfc4Z7Zj4GvHrK8HpgWzW9Dbi2afyBnLYbODciVnarWElSa9q9WmYgM48CZObRiDi/Gl8FvNQ03+Fq7OipK4iITUzv3TMwMECj0WizlN4ZHZrqeB0DZ3dnPf3K/mfvvx//n5ivycnJJdHnTOruvduXQsYMYznTjJk5BowBDA8P58jISJdLqV83LmEcHZpi6/jSvSLV/mfvf+LGkYUtpgcajQb9+P9+N9Tde7tXyxw7ebilen6lGj8MrG6a7wLgSPvlSZLa0W647wA2VNMbgIeaxm+urpq5FHjt5OEbSdLCmfP34Yj4DjACrIiIw8CXgbuAByPiNuBF4Lpq9keAq4GDwFvALTXULEmaw5zhnpk3zPLW5TPMm8DmTouSJHXGO1QlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSrQmZ0sHBETwBvACWAqM4cj4jzgu8AgMAH8aWYe76xMSdJ8dGPPfV1mXpyZw9XrLcCuzFwD7KpeS5IWUB2HZdYD26rpbcC1NXyGJOk0IjPbXzjiBeA4kMA/Z+ZYRPxvZp7bNM/xzFw+w7KbgE0AAwMDn9y+fXvbdfTK+MuvdbyOgbPh2M+7UEyfsv/Z+x9a9eGFLaYHJicnWbZsWa/L6Ilu9L5u3bq9TUdN3qOjY+7AZZl5JCLOB3ZGxH+1umBmjgFjAMPDwzkyMtJhKQtv45aHO17H6NAUW8c73Qz9y/5n73/ixpGFLaYHGo0G/fj/fjfU3XtHh2Uy80j1/ArwA+AS4FhErASonl/ptEhJ0vy0He4RcU5EfOjkNPDHwNPADmBDNdsG4KFOi5QkzU8nvw8PAD+IiJPr+ZfM/FFEPAE8GBG3AS8C13VepiRpPtoO98w8BPz2DOM/Ay7vpChJUme8Q1WSCmS4S1KBDHdJKpDhLkkFMtwlqUBL99ZAaZEb7MId0O2YuOuannyuuss9d0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgrU99/n3qvvvJakxcw9d0kqkOEuSQUy3CWpQH1/zF1Sdy3keazRoSk2Nn2ef7+1e9xzl6QCGe6SVCDDXZIKZLhLUoEMd0kqUG3hHhFXRsSzEXEwIrbU9TmSpPer5VLIiDgDuBv4I+Aw8ERE7MjMZ+r4PEllWEpfJ3LyMtC6Lv+sa8/9EuBgZh7KzP8DtgPra/osSdIpIjO7v9KIzwBXZuZfVK9vAn43M29vmmcTsKl6+evAs10vpD+sAH7a6yJ6yP7tf6n2343eP5qZH5npjbruUI0Zxt7zr0hmjgFjNX1+34iIPZk53Os6esX+7X+p9l9373UdljkMrG56fQFwpKbPkiSdoq5wfwJYExEXRsQHgeuBHTV9liTpFLUclsnMqYi4Hfg34AzgvszcX8dnFWCpH5qy/6VtKfdfa++1nFCVJPWWd6hKUoEMd0kqkOFeo7m+giEifj8inoyIqeregOb3NkTEc9Vjw8JV3T0d9n8iIvZVj747Gd9C71+IiGci4qmI2BURH216byls+9P139fbHlrq/y8jYrzq8fGIWNv03h3Vcs9GxBVtF5GZPmp4MH0i+XngY8AHgf8E1p4yzyDwW8ADwGeaxs8DDlXPy6vp5b3uaaH6r96b7HUPNfe+DvilavqvgO8usW0/Y//9vu3n0f8vN01/CvhRNb22mv8s4MJqPWe0U4d77vWZ8ysYMnMiM58C3jll2SuAnZn5amYeB3YCVy5E0V3USf/9rpXeH83Mt6qXu5m+FwSWzrafrf8StNL/600vz+HdmzzXA9sz8+3MfAE4WK1v3gz3+qwCXmp6fbgaq3vZxaLTHn4xIvZExO6IuLa7pdVuvr3fBvywzWUXo076h/7e9tBi/xGxOSKeB/4O+Ov5LNsK/0B2feb8Coaall0sOu3hVzPzSER8DPj3iBjPzOe7VFvdWu49Iv4cGAb+YL7LLmKd9A/9ve2hxf4z827g7oj4M+BvgA2tLtsK99zr08lXMJTw9Q0d9ZCZR6rnQ0AD+Hg3i6tZS71HxB8CXwI+lZlvz2fZRa6T/vt928P8t+F24ORvKN3b/r0++VDqg+nfig4xfVLk5EmVi2aZ937ef0L1BaZPqC2vps/rdU8L2P9y4KxqegXwHKeckFrMj1Z6ZzqwngfWnDK+JLb9afrv620/j/7XNE3/CbCnmr6I955QPUSbJ1R7/h+i5AdwNfDf1Q/xl6qxv2V6TwXgd5j+l/pN4GfA/qZlb2X6ZMpB4JZe97KQ/QO/B4xXP+TjwG297qWG3n8MHAP2VY8dS2zbz9h/Cdu+xf7/Adhf9f5oc/gz/dvM80x/DfpV7dbg1w9IUoE85i5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoH+H13iiTiNp/QQAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "datos = list()\n", + "n_draws = 100000\n", + "prior = pd.Series(np.random.uniform(0, 1, size=n_draws))\n", + "\n", + "for p in prior:\n", + " datos.append(modelo_generativo(p))\n", + " \n", + "posterior = prior[list(map(lambda x: x == new_users, datos))]\n", + "posterior.hist()" + ] }, { "cell_type": "markdown", @@ -185,10 +327,20 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 22, "metadata": {}, - "outputs": [], - "source": [] + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.08928295277336795 | 0.19081587446568182\n" + ] + } + ], + "source": [ + "print(posterior.quantile(.025),'|',posterior.quantile(.90))" + ] }, { "cell_type": "markdown", @@ -197,6 +349,26 @@ "What is the Maximum Likelihood Estimate?" ] }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Maximum Likelihood Estimate: 0.14 | 0.1342828077314344\n" + ] + } + ], + "source": [ + "rounded = posterior.round(2)\n", + "mode = rounded.mode()[0]\n", + "probability = list(rounded).count(mode)/len(rounded)\n", + "print('Maximum Likelihood Estimate: ', mode, '|',probability)" + ] + }, { "cell_type": "code", "execution_count": null, @@ -221,7 +393,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.2" + "version": "3.7.3" } }, "nbformat": 4,