{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Uma Simulação Simples e o Problema de Monty Hall\n", "### Nelson Kuhl - IME/USP" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np, matplotlib.pyplot as plt # o mínimo\n", "import scipy.stats as st # Há várias ferramentas estatísticas no Python. O módulo stats do scipy integra bem com o\n", " # numpy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exemplo de Simulação com Inteiros\n", "O objeto *randint* do módulo scipy.stats é usado para a distribuição uniforme discreta. Aqui usaremos apenas o *método* rvs, que sorteia valores identicamente distribuidos entre um conjunto discreto de pontos: dados naturais $i$ e $j$ com $i < j-1$ e um inteiro positivo $n$, o comando
\n", "$$\n", "scipy.stats.randint.rvs(i, j, size=n)\n", "$$
\n", "gera $n$ números distribuidos uniformemente no conjunto $\\{i, i+1,\\dots,j-1\\}$. Faremos um teste simples para verificar a distribuição uniforme." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0.321 0.336 0.343]\n" ] } ], "source": [ "n = 1000 # Tamanho da amostra\n", "rs = st.randint.rvs(1, 4, size=n) # n sorteios de números uniformemente distribuidos no conjunto {1, 2, 3}\n", "y = np.array([sum(rs==1), sum(rs==2), sum(rs==3)])/n # array contendo a frequência dos resultados de cada valor\n", "print(y)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAXLklEQVR4nO3dfXBU9b3H8fcXCSI+QEuwOgQI9mp5BiWg4pTiExOwIjNSSyZQY9VY8aGO3nuLYkEpaFHb6+BVkT7ZanxAtB1E0BYrY6uixjY+AMIgQom3I0mkKE1Bke/9Y5d0k2yyJ7LZ3fzyec1k2HN+v5zzzY/dT87+ztkTc3dERKTj65LtAkREJD0U6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigUgZ6Gb2SzPbaWbvtNBuZrbYzLaY2Vtmdkr6yxQRkVSiHKE/CBS30j4JODH+VQ7cf+hliYhIW6UMdHd/EfiolS4XAL/xmHVALzM7Pl0FiohINF3TsI2+wI6E5er4ur837Whm5cSO4jnyyCNHDxo0KA27FxHpPN54441ad++TrC0dgR6Zuy8FlgIUFRV5ZWVlJncvItLhmdn2ltrScZXLB0C/hOWC+DoREcmgdAT6CuA78atdTgN2u3uz6RYREWlfKadczOxRYAKQb2bVwDwgD8DdlwCrgMnAFqAeuKS9ihURkZalDHR3L0nR7sBVaatIRHLKZ599RnV1NXv37s12KZ1K9+7dKSgoIC8vL/L3ZPSkqIh0PNXV1Rx99NEUFhZiZtkup1Nwd+rq6qiurmbgwIGRv08f/ReRVu3du5fevXsrzDPIzOjdu3eb3xUp0EUkJYV55n2RMVegi4gEQoEuIpLgk08+4fbbb++QJ4EV6CKSXhUVUFgIXbrE/q2oSMtmFy9ezODBgyktLU3L9lrygx/8gOHDhzNnzpwvvI25c+eyZs2aNFYVjcWuOsw8ffRfpGPYuHEjgwcPjta5ogLKy6G+/t/revSApUvhEIN40KBBrFmzhoKCgoZ1+/fvp2vXcC/WSzb2ZvaGuxcl668jdBFJnzlzGoc5xJYP4WgX4Hvf+x5bt25l0qRJ9OzZk5kzZ3LGGWcwc+ZMampquPDCCxkzZgxjxozhpZdeAqCuro6JEycydOhQLrvsMgYMGEBtbS3btm1j2LBhDdu+6667uOWWWwB47733KC4uZvTo0Xz961/n3XffBaCsrIxrr72WcePGccIJJ7B8+fKG71+0aBHDhw9n5MiRzJ49u6H/wT7z589nzJgxDBs2jPLycg4eRC9evJghQ4YwYsQIpk+ffkjj08Dds/I1evRoF5Hct2HDhuidzdyh+ZfZIdcxYMAAr6mp8Xnz5vkpp5zi9fX17u5eUlLif/rTn9zdffv27T5o0CB3d7/mmmv81ltvdXf3lStXOuA1NTX+/vvv+9ChQxu2e+edd/q8efPc3f2ss87yzZs3u7v7unXr/Mwzz3R394svvtinTZvmn3/+ua9fv96/+tWvurv7qlWr/PTTT/d//vOf7u5eV1fX0P+JJ55otM7dfcaMGb5ixQp3dz/++ON979697u6+a9eupD9zsrEHKr2FXA33vYqIZF7//rA9yc0A+/dP626mTJnCEUccAcCaNWvYsGFDQ9vHH3/Mnj17ePHFF3nqqacAOO+88/jSl77U6jb37NnDyy+/zLe+9a2Gdfv27Wt4PHXqVLp06cKQIUP48MMPG/Z9ySWX0KNHDwC+/OUvN9vuCy+8wB133EF9fT0fffQRQ4cO5fzzz2fEiBGUlpYydepUpk6d+sUGogkFuoikz8KFyefQFy5M626OPPLIhscHDhxg3bp1dO/ePdL3du3alQMHDjQsH7ya5cCBA/Tq1Yuqqqqk33f44Yc3PPaI5x737t3LrFmzqKyspF+/ftxyyy0N+3vmmWd48cUXefrpp1m4cCFvv/32IZ8P0By6iKRPaWnsBOiAAWAW+zcNJ0RbM3HiRO65556G5YOBPH78eB555BEAVq9eza5duwD4yle+ws6dO6mrq2Pfvn2sXLkSgGOOOYaBAwfyxBNPALHQfvPNN1vd97nnnsuvfvUr6uO/wD76qPEfdzsY3vn5+ezZs6dhXv3AgQPs2LGDM888k0WLFrF792727NlzKMMA6AhdRNKttLRdA7ypxYsXc9VVVzFixAj279/P+PHjWbJkCfPmzaOkpIShQ4cybtw4+senffLy8pg7dy5jx46lb9++JP7ltIqKCq688koWLFjAZ599xvTp0xk5cmSL+y4uLqaqqoqioiK6devG5MmTue222xrae/XqxeWXX86wYcM47rjjGDNmDACff/45M2bMYPfu3bg71157Lb169TrksdBliyLSqjZdtpjDCgsLqaysJD8/P9ulRKbLFkVEOilNuYhIp7Bt27Zsl9DudIQuIhIIBbqISCAU6CIigVCgi4gEQoEuIjnvu9/9Lscee2yjm2oBrFu3jssvv5zXXnuNUaNGMWrUKEaOHMlvf/vbLFWaXQp0Ecl5ZWVlPPvss83Wr169muLiYoYNG0ZlZSVVVVU8++yzXHHFFezfvz8LlWaXLlsUkTaZMGFCs3UXXXQRs2bNor6+nsmTJzdrLysro6ysjNraWqZNm9aobe3atSn3OX78+KSXHT7//PNcf/31DTfHgtjH7Tvr30DVEbqIdEi1tbXk5eXRs2dPAF599VWGDh3K8OHDWbJkSdB/+KIlne8nFpFD0toRdY8ePVptz8/Pj3REHsXvf/97Jk6c2LB86qmnsn79ejZu3MjFF1/MpEmTIt+BMRQ6QheRDung/HlTgwcP5qijjuKdd97JQlXZpUAXkQ7H3XnrrbcYNWoUAO+//37DSdDt27fz7rvvUlhYmL0Cs0RTLiKS80pKSli7di21tbUUFBRwzTXXcPLJJzec/Pzzn//Mj3/8Y/Ly8ujSpQv33Xdfh7qrYroo0EUk5z366KONlhcsWNBoumXmzJnMnDkz02XlHAW6iHQ4N998c7ZLyEmaQxcRCYQCXURSytZfNuvMvsiYK9BFpFXdu3enrq5OoZ5B7k5dXV2br6PXHLqItKqgoIDq6mpqamqyXUqn0r17dwoKCtr0PQp0EWlVXl4eAwcOzHYZEkGkKRczKzazTWa2xcxmJ2nvb2YvmNlfzewtM2t+dx4REWlXKQPdzA4D7gUmAUOAEjMb0qTbzcAydz8ZmA7cl+5CRUSkdVGO0McCW9x9q7t/CjwGXNCkjwPHxB/3BP4vfSWKiEgUUQK9L7AjYbk6vi7RLcAMM6sGVgHXJNuQmZWbWaWZVeoEi4hIeqXrssUS4EF3LwAmAw+ZWbNtu/tSdy9y96I+ffqkadciIgLRAv0DoF/CckF8XaJLgWUA7v4K0B3ofHfGERHJoiiB/jpwopkNNLNuxE56rmjS52/A2QBmNphYoGtORUQkg1IGurvvB64GngM2EruaZb2ZzTezKfFuNwCXm9mbwKNAmetjZSIiGRXpg0XuvorYyc7EdXMTHm8AzkhvaSIi0ha6l4uISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigYgU6GZWbGabzGyLmc1uoc9FZrbBzNab2SPpLVNERFLpmqqDmR0G3AucC1QDr5vZCnffkNDnROBG4Ax332Vmx7ZXwSIiklyUI/SxwBZ33+runwKPARc06XM5cK+77wJw953pLVNERFKJEuh9gR0Jy9XxdYlOAk4ys5fMbJ2ZFSfbkJmVm1mlmVXW1NR8sYpFRCSpdJ0U7QqcCEwASoCfmVmvpp3cfam7F7l7UZ8+fdK0axERgWiB/gHQL2G5IL4uUTWwwt0/c/f3gc3EAl5ERDIkSqC/DpxoZgPNrBswHVjRpM/viB2dY2b5xKZgtqavTBERSSVloLv7fuBq4DlgI7DM3deb2XwzmxLv9hxQZ2YbgBeA/3L3uvYqWkREmjN3z8qOi4qKvLKyMiv7FhHpqMzsDXcvStamT4qKiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIICIFupkVm9kmM9tiZrNb6XehmbmZFaWvRBERiSJloJvZYcC9wCRgCFBiZkOS9Dsa+D7warqLFBGR1KIcoY8Ftrj7Vnf/FHgMuCBJvx8Bi4C9aaxPREQiihLofYEdCcvV8XUNzOwUoJ+7P9Pahsys3MwqzayypqamzcWKiEjLDvmkqJl1AX4K3JCqr7svdfcidy/q06fPoe5aREQSRAn0D4B+CcsF8XUHHQ0MA9aa2TbgNGCFToyKiGRWlEB/HTjRzAaaWTdgOrDiYKO773b3fHcvdPdCYB0wxd0r26ViERFJKmWgu/t+4GrgOWAjsMzd15vZfDOb0t4FiohINF2jdHL3VcCqJuvmttB3wqGXJSIibaVPioqIBEKBLiISCAW6iEggFOgi8sVVVEBhIXTpEvu3oiLbFXVqkU6Kiog0U1EB5eVQXx9b3r49tgxQWpq9ujoxHaGLJNIRZ3Rz5vw7zA+qr4+tl6zQEbrIQTribJu//a1t66XddchAnzBhQrN1F110EbNmzaK+vp7Jkyc3ay8rK6OsrIza2lqmTZvWrP3KK6/k29/+Njt27GDmzJnN2m+44QbOP/98Nm3axBVXXNGs/eabb+acc86hqqqK6667rln7bbfdxrhx43j55Ze56aabmrXffffdjBo1ijVr1rBgwYJm7Q888ABf+9rXePrpp/nJT37SrP2hhx6iX79+PP7449x///3N2pcvX05+fj4PPvggDz74YLP2VatW0aNHD+677z6WLVvWrH3t2rUA3HXXXaxcubJR2xFHHMHq1asB+NGPfsTzzz/fqL137948+eSTANx444288sorjdoLCgp4+OGHAbjuuuuoqqpq1H7SSSexdOlSAMrLy9m8eXOj9lGjRnH33XcDMGPGDKqrqxu1n3766dx+++0AXHjhhdTV1TVqP/vss/nhD38Ic+Ywqb6efyU21tfzzauv5j/jga7nXsJzr1s32LePh4jdG+Rx4P6D6xPGSc+95s+9gz9TumnKJXQffgjr1kGfPrEphCZPaEnQ0pHlP/6R0TI6jIEDY1NTibp0ia2XrDB3z8qOi4qKvLJSt3tpV02nEAB69IClSzWFkExhYWyapakBA2DbtkxX0zFUVMCll8K+fbFxWrhQz612ZmZvuHvSmx/qCD1kOmnVNgsXxn7hJerRI7ZekisthdNOg298I/ZLT2GeVQr0kOmkVduUlsbevRx+eGx5wAC9m5EOpUOeFJWI+vdPPoXQv3/ma+koSkvhZz+LPW6nE1ci7UVH6CHTFIJIp6JAD5mmEEQ6FU25hE5TCCKdho7QRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQlEpEA3s2Iz22RmW8xsdpL2681sg5m9ZWbPm9mA9JcqIiKtSRnoZnYYcC8wCRgClJjZkCbd/goUufsIYDlwR7oLFRGR1kU5Qh8LbHH3re7+KfAYcEFiB3d/wd3r44vrgIL0likiIqlECfS+wI6E5er4upZcCqxO1mBm5WZWaWaVNTU10asUEZGU0npS1MxmAEXAncna3X2puxe5e1GfPn3SuWsRkU6va4Q+HwD9EpYL4usaMbNzgDnAN9x9X3rKExGRqKIcob8OnGhmA82sGzAdWJHYwcxOBh4Aprj7zvSXKSIiqaQMdHffD1wNPAdsBJa5+3ozm29mU+Ld7gSOAp4wsyozW9HC5kREpJ1EmXLB3VcBq5qsm5vw+Jw01yUiIm2kT4qKiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIICIFupkVm9kmM9tiZrOTtB9uZo/H2181s8K0VyoiIq1KGehmdhhwLzAJGAKUmNmQJt0uBXa5+38A/wMsSnehIiLSuihH6GOBLe6+1d0/BR4DLmjS5wLg1/HHy4GzzczSV6aIiKTSNUKfvsCOhOVq4NSW+rj7fjPbDfQGahM7mVk5UB5f3GNmm75I0UB+023niFytCyAfs1ysLVfHTOPVNhqvtjmUuga01BAl0NPG3ZcCSw91O2ZW6e5FaSgprXK1Lsjd2lRX26iutulsdUWZcvkA6JewXBBfl7SPmXUFegJ16ShQRESiiRLorwMnmtlAM+sGTAdWNOmzArg4/nga8Ed39/SVKSIiqaScconPiV8NPAccBvzS3deb2Xyg0t1XAL8AHjKzLcBHxEK/PR3ytE07ydW6IHdrU11to7raplPVZTqQFhEJgz4pKiISCAW6iEggcjrQc/WWAxHqKjOzGjOrin9dlqG6fmlmO83snRbazcwWx+t+y8xOyZG6JpjZ7oTxmpuhuvqZ2QtmtsHM1pvZ95P0yfiYRawr42NmZt3N7DUzezNe161J+mT8NRmxrqy8JuP7PszM/mpmK5O0pXe83D0nv4idgH0POAHoBrwJDGnSZxawJP54OvB4jtRVBvxvFsZsPHAK8E4L7ZOB1YABpwGv5khdE4CVWRiv44FT4o+PBjYn+b/M+JhFrCvjYxYfg6Pij/OAV4HTmvTJxmsySl1ZeU3G93098Eiy/690j1cuH6Hn6i0HotSVFe7+IrGrjFpyAfAbj1kH9DKz43Ogrqxw97+7+1/ijz8BNhL71HOijI9ZxLoyLj4Ge+KLefGvpldVZPw1GbGurDCzAuA84OctdEnreOVyoCe75UDTJ3WjWw4AB285kO26AC6Mv0Vfbmb9krRnQ9Tas+H0+Fvm1WY2NNM7j7/VPZnY0V2irI5ZK3VBFsYsPn1QBewE/uDuLY5XBl+TUeqC7Lwm7wb+GzjQQntaxyuXA70jexoodPcRwB/4929gSe4vwAB3HwncA/wukzs3s6OAJ4Hr3P3jTO67NSnqysqYufvn7j6K2CfGx5rZsEzsN5UIdWX8NWlm3wR2uvsb7b2vg3I50HP1lgMp63L3OnffF1/8OTC6nWuKKsqYZpy7f3zwLbO7rwLyzCw/E/s2szxioVnh7k8l6ZKVMUtVVzbHLL7PfwAvAMVNmrJ6G5CW6srSa/IMYIqZbSM2NXuWmT3cpE9axyuXAz1XbzmQsq4mc6xTiM2B5oIVwHfiV26cBux2979nuygzO+7gvKGZjSX2vGz3EIjv8xfARnf/aQvdMj5mUerKxpiZWR8z6xV/fARwLvBuk24Zf01GqSsbr0l3v9HdC9y9kFhO/NHdZzTpltbxyujdFtvCc/OWA1HrutbMpgD743WVtXddAGb2KLGrH/LNrBqYR+wEEe6+BFhF7KqNLUA9cEmO1DUNuNLM9gP/AqZn4BczxI6gZgJvx+dfAW4C+ifUlo0xi1JXNsbseODXFvujN12AZe6+MtuvyYh1ZeU1mUx7jpc++i8iEohcnnIREZE2UKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEoj/B6ZfjZaRrisZAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Gráfico\n", "\n", "x = np.array([1, 2, 3]) # abscissas. As ordenadas estão em y\n", "plt.plot(x, y, 'ro', label=\"frequências\")\n", "plt.vlines(x, 0, y, color='r', linestyles='-')\n", "plt.ylim(0,1)\n", "xx = np.arange(0, 4, 0.01); yy = np.ones(np.size(xx))/3 # Para a ver a linha 1/3 na figura.\n", "plt.plot(xx, yy, '--k', label=\"1/3\") # Deve haver maneira melhor\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## O Problema de Monty Hall\n", "Vamos relembrar o enunciado. Em um programa de auditório, há três portas fechadas, sendo que uma delas contém um prêmio e nas outras duas há duas cabras. O apresentador pede para o participante escolher uma porta. Após a escolha, o apresentador abre uma porta contendo uma cabra e pergunta ao participante se ele quer trocar de porta.\n", "\n", "Simularemos aqui a situação na qual o participante sempre troca de porta para estimar a probabilidade de ganhar o prêmio. A ideia é ignorar solução discutida em aula e resolver o problema usando força bruta. Não é elegante, mas funciona.\n", "\n", "Note que se o participante escolher a porta que tem o prêmio, o apresentador pode abrir qualquer uma das outras portas e, ao mudar de porta, o participante obterá uma cabra. Caso a escolha inicial tenha sido uma porta contendo uma cabra, só há uma escolha para o apresentador abrir e, ao mudar a porta, o participante ganha o prêmio.\n", "\n", "As portas serão numeradas por 1, 2 e 3. Para a simulação, usaremos um vetor $p$ no qual a componente a componente $p_i$ armazena o número da porta contendo o prêmio no $i$-ésimo desafio, e um vetor $x$ cuja componente $x_i$ contém a escolha do participante no $i$-ésimo desafio. O prêmio é ganho de $x_i \\ne p_i$. Os valores de $x$ serão sorteados entre os números 1, 2 e 3 com mesma probabilidade. Os valores de $p$ podem ser atribuídos arbitrariamente. Usaremos também a distribuição uniforme entre os valores 1, 2 e 3, mas você pode testar outras maneiras." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "A probabilidade de ganhar o prêmio é estimada por 0.66539\n" ] } ], "source": [ "n = 100000 # número de sorteios (modifique se quiser)\n", "p = st.randint.rvs(1, 4, size=n) # sorteio da porta que contém o prêmio\n", "x = st.randint.rvs(1, 4, size=n) # escolha inicial do participante\n", "ganhou = x != p # vetor Booleano armazenando o resultado do desafio\n", "prob_ganhar = sum(ganhou)/n # frequência de ganho\n", "print(\"A probabilidade de ganhar o prêmio é estimada por %.5f\" % prob_ganhar)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Observações\n", "\n", "- Se executarmos a célula acima várias vêzes com o mesmo valor de $n$, obteremos resultados diferentes, mas próximos. Há uma **margem de erro** para o resultado, e veremos adiante como estimá-la.\n", "- A simulação acima parece fácil e é. Podemos aproximar a solução do problema sem pensar muito, deixando o trabalho braçal para o computador. Porém as hipóteses ficam mascaradas" ] }, { "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.6.9" } }, "nbformat": 4, "nbformat_minor": 2 }