diff --git a/ipynb/Dice Baseball.ipynb b/ipynb/Dice Baseball.ipynb new file mode 100644 index 0000000..17e829e --- /dev/null +++ b/ipynb/Dice Baseball.ipynb @@ -0,0 +1,253 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
Peter Norvig
March 2019
\n", + "\n", + "# Dice Baseball\n", + "\n", + "The [March 22, 2019 Riddler](https://fivethirtyeight.com/features/can-you-turn-americas-pastime-into-a-game-of-yahtzee/) asks us to simulate baseball using probabilities from a 19th century dice game. The simulation is pretty straightforward. (Note that I chose not to implement a strike as an event, but rather to count the probability of getting 3 strikes in a row, and calling that an event.)" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import matplotlib.pyplot as plt\n", + "from collections import Counter\n", + "from statistics import mean\n", + "import random" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "def inning(events='2111111EEBBOOooooooofffffD334', strikes=7*'s', verbose=False):\n", + " \"Play a random inning and return number of runs scored.\"\n", + " assert len(events) + len(strikes) == 6 * 6\n", + " K = (len(strikes) / len(strikes+events)) ** 3 # The probability of a strikeout\n", + " events = events.replace('E', '1') # An error is the same as a single\n", + " outs = runs = 0\n", + " runners = [] # A list of runners; [1, 2, 3] means bases loaded\n", + " while True:\n", + " x = 'K' if random.random() <= K else random.choice(events)\n", + " if verbose: \n", + " print('{} outs; {} runs; runners on {}; result is {}'\n", + " .format(outs, runs, runners, x))\n", + " if x in '1234': # single, double, triple, homer\n", + " runners = [r + int(x) + (r == 2) for r in runners + [0]]\n", + " elif x == 'B': # base on balls\n", + " runners = [(r+1 if r==1 or r-1 in runners else r) for r in runners] + [1]\n", + " elif x in 'KOofD': # srikeout, foul out, out at first, fly out, double play\n", + " outs += 1\n", + " if x == 'D' and runners: # double play\n", + " runners.remove(min(runners))\n", + " outs += 1\n", + " elif x == 'o': # out at first (other runners advance)\n", + " runners = [r+1 for r in runners]\n", + " elif x == 'f' and 3 in runners: # fly out; runner on 3rd scores\n", + " runs += 1\n", + " runners.remove(3)\n", + " else:\n", + " raise ValueError('unknown events: ' + x)\n", + " if outs >= 3:\n", + " return runs\n", + " else: # See if anybody scored\n", + " runs += sum(r >= 4 for r in runners)\n", + " runners = [r for r in runners if r < 4]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's peek at some random innings:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0 outs; 0 runs; runners on []; result is o\n", + "1 outs; 0 runs; runners on []; result is o\n", + "2 outs; 0 runs; runners on []; result is 1\n", + "2 outs; 0 runs; runners on [1]; result is 1\n", + "2 outs; 0 runs; runners on [2, 1]; result is O\n" + ] + }, + { + "data": { + "text/plain": [ + "0" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "inning(verbose=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0 outs; 0 runs; runners on []; result is 3\n", + "0 outs; 0 runs; runners on [3]; result is f\n", + "1 outs; 1 runs; runners on []; result is 1\n", + "1 outs; 1 runs; runners on [1]; result is o\n", + "2 outs; 1 runs; runners on [2]; result is 1\n", + "2 outs; 2 runs; runners on [1]; result is 1\n", + "2 outs; 2 runs; runners on [2, 1]; result is K\n" + ] + }, + { + "data": { + "text/plain": [ + "2" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "inning(verbose=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0 outs; 0 runs; runners on []; result is o\n", + "1 outs; 0 runs; runners on []; result is 3\n", + "1 outs; 0 runs; runners on [3]; result is D\n" + ] + }, + { + "data": { + "text/plain": [ + "0" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "inning(verbose=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "That looks good.\n", + "\n", + "Now, simulate a million innings, use them to simulate a million games, and show a histogram of runs scored per team per game:" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 20.7 s, sys: 73.7 ms, total: 20.8 s\n", + "Wall time: 20.9 s\n" + ] + } + ], + "source": [ + "%%time\n", + "N = 1000000\n", + "innings = [inning() for _ in range(N)]\n", + "games = [sum(random.choice(innings) for i in range(9)) for g in range(N)]" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(15.106058, 15.10542)" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD8CAYAAACcjGjIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAFURJREFUeJzt3X+MXeV95/H3pwYalJYagiFej11T1duGRhuSWGCJ1SqF1hga1VQKLVFbHJaVpYhUVOqqIVUlVBKkRKqaH1KWlQVeTJUWUNosVmTqug6ou1JCPAQaAhThUoIHG+yugdKNAgv97h/3cbnxufZcz9hz5868X9LVved7n3vmOWLwZ55znvucVBWSJPX7sVF3QJI0/xgOkqQOw0GS1GE4SJI6DAdJUofhIEnqMBwkSR2GgySpw3CQJHWcNuoOzNS5555bq1evHnU3JGlsPPLII/9UVcuGaTu24bB69WomJydH3Q1JGhtJvj9sW08rSZI6DAdJUofhIEnqMBwkSR2GgySpw3CQJHUYDpKkDsNBktRhOEiSOgwHAbB8YhVJOo/lE6tG3TVJIzC2y2fo5HrxhX389Ce/3ql//49/jSQ/Unv3ipUcmHp+rromaQQMBx3fW/+vExrf/9yHR9QZSXPF00qL0KBTSJLUz5HDIjToFJKjAUn9HDlIkjoMB524Jac7q0la4DytpBPnRWppwXPkIEnqMBwkSR2GgySpw3CQJHUYDpKkDsNBktRhOEiSOgyHBc51lCTNxFBfgkuyFLgDeC9QwH8GngbuBVYDzwG/XlUvp/evzxeBq4AfAB+rqu+0/WwC/rDt9jNVta3VPwjcBZwJ7ABuqqqa/eHJdZQkzcSwI4cvAn9VVT8PvA94CrgZ2F1Va4DdbRvgSmBNe2wGbgdIcg5wC3AJcDFwS5Kz22dub22PfG7D7A5Lc27AkhouqyGNr2lHDknOAv4T8DGAqnoDeCPJRuBDrdk24CHgk8BG4O72l/+3kixNsry13VVVh9t+dwEbkjwEnFVV32z1u4GrgQdOyhFqbgxYUgMcpUjjapiRw88Ah4D/keTRJHckeSdwflUdAGjP57X2K4B9fZ+farXj1acG1DuSbE4ymWTy0KFDQ3RdkjQTw4TDacAHgNur6v3A/+XtU0iDDLriWTOod4tVW6pqbVWtXbZs2fF7LUmasWHCYQqYqqqH2/ZX6YXFS+10Ee35YF/7lX2fnwD2T1OfGFCXJI3ItOFQVS8C+5L8XCtdDjwJbAc2tdom4P72ejtwXXrWAa+20047gfVJzm4XotcDO9t7ryVZ12Y6Xde3L0nSCAx7P4ffAb6S5AzgWeB6esFyX5IbgOeBa1rbHfSmse6lN5X1eoCqOpzk08Ce1u7WIxengY/z9lTWB/BitCSN1FDhUFWPAWsHvHX5gLYF3HiM/WwFtg6oT9L7DoUkaR7wG9KSpA7DQZLUYThIkjoMB0lSh+EgSeowHCRJHYaDTq0Bq7W6Uqs0/w37JThpZgas1upKrdL858hBktRhOCwQg24H6i1BJc2Up5UWiEG3AwVP4UiaGUcOkqQOw0GS1GE4SJI6DAdJUofhIEnqMBwkSR2GgySpw3CQJHUYDpKkDsNBktRhOGjuuYy3NO8NtbZSkueA14C3gDeram2Sc4B7gdXAc8CvV9XL6a329kXgKuAHwMeq6jttP5uAP2y7/UxVbWv1DwJ3AWcCO4CbqqpOwvFpPnIZb2neO5GRwy9W1UVVtbZt3wzsrqo1wO62DXAlsKY9NgO3A7QwuQW4BLgYuCXJ2e0zt7e2Rz63YcZHJEmatdmcVtoIbGuvtwFX99Xvrp5vAUuTLAeuAHZV1eGqehnYBWxo751VVd9so4W7+/YlSRqBYcOhgL9O8kiSza12flUdAGjP57X6CmBf32enWu149akBdUnSiAx7P4dLq2p/kvOAXUn+/jhtB91hpmZQ7+64F0ybAVat8gKmJJ0qQ40cqmp/ez4IfI3eNYOX2ikh2vPB1nwKWNn38Qlg/zT1iQH1Qf3YUlVrq2rtsmXLhum6JGkGpg2HJO9M8pNHXgPrge8B24FNrdkm4P72ejtwXXrWAa+20047gfVJzm4XotcDO9t7ryVZ12Y6Xde3L0nSCAxzWul84GvtfsSnAX9WVX+VZA9wX5IbgOeBa1r7HfSmse6lN5X1eoCqOpzk08Ce1u7WqjrcXn+ct6eyPtAeOoblE6t48YV90zeUpBmaNhyq6lngfQPq/we4fEC9gBuPsa+twNYB9UngvUP0Vwy+X7TfE5B0MvkNaUlSh+EgSeowHCRJHYaDJKnDcJAkdRgOkqQOw0GS1GE4SJI6DAfNDwPuDucd4qTRGXZVVunUGnB3OPCb39KoOHKQJHUYDpKkDsNBktRhOEiSOgwHSVKH4SBJ6jAcJEkdhoMkqcNwkCR1GA6SpA7DQZLUYThIkjqGDockS5I8muTrbfuCJA8neSbJvUnOaPUfb9t72/ur+/bxqVZ/OskVffUNrbY3yc0n7/DG3/KJVZ2VSiXpVDuRVVlvAp4CzmrbnwM+X1X3JPnvwA3A7e355ar62STXtna/keRC4FrgF4B/B/xNkn/f9vVl4JeBKWBPku1V9eQsj21BePGFfZ3VSl2pVNKpNtTIIckE8CvAHW07wGXAV1uTbcDV7fXGtk17//LWfiNwT1W9XlX/COwFLm6PvVX1bFW9AdzT2kqSRmTY00pfAH4f+Ne2/S7glap6s21PASva6xXAPoD2/qut/b/Vj/rMseqSpBGZNhySfBg4WFWP9JcHNK1p3jvR+qC+bE4ymWTy0KFDx+m1FowBd4jz7nDSqTfMNYdLgV9NchXwDnrXHL4ALE1yWhsdTAD7W/spYCUwleQ04KeAw331I/o/c6z6j6iqLcAWgLVr1w4MEC0wA+4Q5zUX6dSbduRQVZ+qqomqWk3vgvI3quo3gQeBj7Rmm4D72+vtbZv2/jeqqlr92jab6QJgDfBtYA+wps1+OqP9jO0n5egkSTMym3tIfxK4J8lngEeBO1v9TuBPk+ylN2K4FqCqnkhyH/Ak8CZwY1W9BZDkE8BOYAmwtaqemEW/JEmzdELhUFUPAQ+118/Sm2l0dJsfAtcc4/O3AbcNqO8AdpxIXyRJp47fkJYkdRgOkqQOw0GS1GE4SJI6DAdJUofhIEnqMBwkSR2GgySpw3CQJHUYDpKkDsNBktRhOEiSOgwHjR9vACSdcrNZslsaDW8AJJ1yjhwkSR2GgySpw3CQJHUYDpKkDsNhnlg+saozAyfJqLslaZFyttI88eIL+zozcMBZOJJGw5GDJKnDcJAkdUwbDknekeTbSf4uyRNJ/qjVL0jycJJnktyb5IxW//G2vbe9v7pvX59q9aeTXNFX39Bqe5PcfPIPU5J0IoYZObwOXFZV7wMuAjYkWQd8Dvh8Va0BXgZuaO1vAF6uqp8FPt/akeRC4FrgF4ANwH9LsiTJEuDLwJXAhcBHW1tJ0ohMGw7V8y9t8/T2KOAy4Kutvg24ur3e2LZp71+e3rSbjcA9VfV6Vf0jsBe4uD32VtWzVfUGcE9rK0kakaGuObS/8B8DDgK7gH8AXqmqN1uTKWBFe70C2AfQ3n8VeFd//ajPHKsuSRqRocKhqt6qqouACXp/6b9nULP2PGhyfs2g3pFkc5LJJJOHDh2avuOSpBk5odlKVfUK8BCwDlia5Mj3JCaA/e31FLASoL3/U8Dh/vpRnzlWfdDP31JVa6tq7bJly06k65KkEzDMbKVlSZa212cCvwQ8BTwIfKQ12wTc315vb9u0979RVdXq17bZTBcAa4BvA3uANW320xn0LlpvPxkHJ0mamWG+Ib0c2NZmFf0YcF9VfT3Jk8A9ST4DPArc2drfCfxpkr30RgzXAlTVE0nuA54E3gRurKq3AJJ8AtgJLAG2VtUTJ+0ItTi0GwAd7d0rVnJg6vkRdEgab9OGQ1V9F3j/gPqz9K4/HF3/IXDNMfZ1G3DbgPoOYMcQ/ZUGG3ADIHD5EWmm/Ia0JKnDcJAkdRgOkqQOw0GS1GE4SJI6DAdJUofhIEnqMBwkSR2GgySpw3CQJHUYDpKkDsNBktRhOGhha6u19j+WT6wada+keW+YJbt1ki2fWMWLL+ybvqFmb8Bqra7UKk3PcBiBF1/Y5z9YkuY1TytJkjoMB0lSh+EgSeowHCRJHYaDJKnDcJAkdRgOkqQOw0GS1DFtOCRZmeTBJE8leSLJTa1+TpJdSZ5pz2e3epJ8KcneJN9N8oG+fW1q7Z9Jsqmv/sEkj7fPfClJTsXBSpKGM8zI4U3g96rqPcA64MYkFwI3A7urag2wu20DXAmsaY/NwO3QCxPgFuAS4GLgliOB0tps7vvchtkfmiRppqYNh6o6UFXfaa9fA54CVgAbgW2t2Tbg6vZ6I3B39XwLWJpkOXAFsKuqDlfVy8AuYEN776yq+mZVFXB3374kSSNwQtcckqwG3g88DJxfVQegFyDAea3ZCqB/VbmpVjtefWpAfdDP35xkMsnkoUOHTqTr0ttcqVWa1tAL7yX5CeAvgN+tqn8+zmWBQW/UDOrdYtUWYAvA2rVrB7aRpuVKrdK0hho5JDmdXjB8par+spVfaqeEaM8HW30KWNn38Qlg/zT1iQF1SdKIDDNbKcCdwFNV9Sd9b20Hjsw42gTc31e/rs1aWge82k477QTWJzm7XYheD+xs772WZF37Wdf17UuSNALDnFa6FPht4PEkj7XaHwCfBe5LcgPwPHBNe28HcBWwF/gBcD1AVR1O8mlgT2t3a1Udbq8/DtwFnAk80B6SpBGZNhyq6n8z+LoAwOUD2hdw4zH2tRXYOqA+Cbx3ur5IkuaG35CWJHUYDpKkDsNBktRhOEiSOgwHSVKH4XAKLZ9Y1VmmwQVn56kBS2q4rIYWs6GXz9CJe/GFfZ1lGsClGualAUtqgP+ttHg5cpAkdRgOkqQOw0GS1GE4SJI6DAdJUofhIEnqMBwkSR2GgySpw3CQJHUYDtLxDFhWwyU1tBi4fIZ0PAOW1XBJDS0GjhwkSR2GgySpw3CQJHUYDpKkjmnDIcnWJAeTfK+vdk6SXUmeac9nt3qSfCnJ3iTfTfKBvs9sau2fSbKpr/7BJI+3z3wp3g1HkkZumJHDXcCGo2o3A7urag2wu20DXAmsaY/NwO3QCxPgFuAS4GLgliOB0tps7vvc0T9LkjTHpg2Hqvpb4PBR5Y3AtvZ6G3B1X/3u6vkWsDTJcuAKYFdVHa6ql4FdwIb23llV9c2qKuDuvn2NlUG3BNUC5XcftAjM9HsO51fVAYCqOpDkvFZfAezrazfVaserTw2oj51BtwR1PvwC5XcftAic7AvSg/5crhnUB+882ZxkMsnkoUOHZthFSdJ0ZhoOL7VTQrTng60+BazsazcB7J+mPjGgPlBVbamqtVW1dtmyZTPsuiRpOjMNh+3AkRlHm4D7++rXtVlL64BX2+mnncD6JGe3C9HrgZ3tvdeSrGuzlK7r25ckaUSmveaQ5M+BDwHnJpmiN+vos8B9SW4Angeuac13AFcBe4EfANcDVNXhJJ8G9rR2t1bVkYvcH6c3I+pM4IH2kCSN0LThUFUfPcZblw9oW8CNx9jPVmDrgPok8N7p+iFJmjt+Q1o6GQZMb3WKq8aZS3ZLJ8OA6a3gFFeNL0cOkqQOw0GS1GE4SJI6DAdJUofhIJ1KLtKnMeVsJelUcpE+jSlHDjPg8tySFjpHDjPg8tySFjpHDtJc8zqExoAjB2mueR1CY8CRgySpw3CQJHUYDtJ84Kqumme85iDNB67qqnnGkYMkqcNwOI5BX3bzC2+aU0571Yh4Wuk4Bn3ZDRzqaw457VUj4shBGjeOJjQHHDlI48bRhOaAIwdpIXAqrE6yeTNySLIB+CKwBLijqj474i5J4+NYU2H/+Nc6kyjevWIlB6aen6ueaUzNi3BIsgT4MvDLwBSwJ8n2qnpyrvqwfGIVL76wb65+nDQ3Bp2CMjA0hHkRDsDFwN6qehYgyT3ARmDOwsFluLVoDBkYAEvOeAdvvfHDH6kZJIvDfAmHFUD/n+1TwCWn6oc5SpCOcpxvaA8TJINC5Fh1A2c8pKpG3QeSXANcUVX/pW3/NnBxVf3OUe02A5vb5s8BT8/wR54L/NMMPzufeVzjZ6Ee20I9LhjvY/vpqlo2TMP5MnKYAlb2bU8A+49uVFVbgC2z/WFJJqtq7Wz3M994XONnoR7bQj0uWNjH1m++TGXdA6xJckGSM4Brge0j7pMkLVrzYuRQVW8m+QSwk95U1q1V9cSIuyVJi9a8CAeAqtoB7JijHzfrU1PzlMc1fhbqsS3U44KFfWz/Zl5ckJYkzS/z5ZqDJGkeWVThkGRDkqeT7E1y86j7MxtJtiY5mOR7fbVzkuxK8kx7PnuUfZyJJCuTPJjkqSRPJLmp1cf62JK8I8m3k/xdO64/avULkjzcjuveNiFj7CRZkuTRJF9v2wvluJ5L8niSx5JMttpY/y4Oa9GEQ98SHVcCFwIfTXLhaHs1K3cBG46q3Qzsrqo1wO62PW7eBH6vqt4DrANubP+dxv3YXgcuq6r3ARcBG5KsAz4HfL4d18vADSPs42zcBDzVt71QjgvgF6vqor7pq+P+uziURRMO9C3RUVVvAEeW6BhLVfW3wOGjyhuBbe31NuDqOe3USVBVB6rqO+31a/T+wVnBmB9b9fxL2zy9PQq4DPhqq4/dcQEkmQB+BbijbYcFcFzHMda/i8NaTOEwaImOFSPqy6lyflUdgN4/ssB5I+7PrCRZDbwfeJgFcGzt1MtjwEFgF/APwCtV9WZrMq6/k18Afh/417b9LhbGcUEvwP86ySNthQZYAL+Lw5g3U1nnwKCbPztVa55K8hPAXwC/W1X/vBDu3V1VbwEXJVkKfA14z6Bmc9ur2UnyYeBgVT2S5ENHygOajtVx9bm0qvYnOQ/YleTvR92hubKYRg5DLdEx5l5KshygPR8ccX9mJMnp9ILhK1X1l628II4NoKpeAR6id01laZIjf6SN4+/kpcCvJnmO3qnay+iNJMb9uACoqv3t+SC9QL+YBfS7eDyLKRwWwxId24FN7fUm4P4R9mVG2vnqO4GnqupP+t4a62NLsqyNGEhyJvBL9K6nPAh8pDUbu+Oqqk9V1URVrab3/9Q3quo3GfPjAkjyziQ/eeQ1sB74HmP+uzisRfUluCRX0fur5sgSHbeNuEszluTPgQ/RWyHyJeAW4H8C9wGrgOeBa6rq6IvW81qS/wj8L+Bx3j6H/Qf0rjuM7bEl+Q/0Ll4uofdH2X1VdWuSn6H3F/c5wKPAb1XV66Pr6cy100r/tao+vBCOqx3D19rmacCfVdVtSd7FGP8uDmtRhYMkaTiL6bSSJGlIhoMkqcNwkCR1GA6SpA7DQZLUYThIkjoMB0lSh+EgSer4/4rD2ben7SopAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.hist(games, ec='black', bins=max(games)-min(games)+1)\n", + "mean(games), mean(innings) * 9" + ] + } + ], + "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.5.3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}