{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Playing Darts and Estimating Pi\n", "You and some friends are playing darts. You all are **very** bad at darts. Each throw is guaranteed to hit the square board depicted below, but otherwise each throw will land in a completely random position within the square. To entertain yourself during this pathetic display, you decide to use this as an opportunity to estimate the irrational number $\\pi \\approx 3.14159$.\n", "\n", "![A dart board](attachments/circle_square_small.png)\n", "\n", "Because each throw falls randomly within the square, you realize that the probability of a dart landing within the circle is given by the ratio of the circle's area to the square's area:\n", "\n", "\\begin{equation}\n", "P_{circle} = \\frac{Area_{circle}}{Area_{square}} = \\frac{\\pi r^2}{(2r)^2}\n", "\\end{equation}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Furthermore, we can interpret $P_{circle}$ as being approximated by the fraction of darts thrown that land in the circle. Thus, we find:\n", "\n", "\\begin{equation}\n", "\\frac{N_{circle}}{N_{total}} \\approx \\frac{\\pi r^2}{(2r)^2} = \\frac{\\pi}{4}\n", "\\end{equation}\n", "\n", "where $N_{total}$ is the total number of darts thrown, and $N_{circle}$ is the number of darts that land within the circle. Thus simply by keeping tally of where the darts land, you can begin to estimate the value of $\\pi$!\n", "\n", "## Problem 1\n", "\n", "Write code that simulates the dart throwing and tallying process, and keep a running estimate of $\\pi$ as \"darts are being thrown\". For simplicity, you can assume that the board is centered at $(0, 0)$, and that $r = 1$ (the radius of the circle). Use `numpy.random.rand` ([link to docs](https://numpy.org/doc/stable/reference/generated/numpy.random.rand.html)) to randomly generate the positions on the board where the darts land. Do this for $N = 10,000$ darts in total. For each dart thrown determine whether or not it landed within the circle, and update your estimate of $\\pi$ according to the formula: $N_{circle} / N_{total} \\approx \\pi / 4$\n", "\n", "Keep in mind that each dart can land in $(x \\in [-1, 1], y \\in [-1, 1])$ and that a dart that lands at $(x, y)$ falls within the circle if \n", "\n", "\\begin{equation}\n", "\\sqrt{x^2 + y^2} < 1\n", "\\end{equation}\n", "\n", "You can start this problem by writing a solution that uses explicit for-loops to help make clear the solution. That being said, you should strive to write a fully-vectorized solution (i.e. compute the running estimation of $\\pi$ during $10,000$ dart throws without any explicit for-loops)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Relevant Reading\n", "\n", "You will need to be familiar with NumPy's [vectorized operations](http://www.pythonlikeyoumeanit.com/Module3_IntroducingNumpy/VectorizedOperations.html) and with [summing over axes](http://www.pythonlikeyoumeanit.com/Module3_IntroducingNumpy/VectorizedOperations.html#Specifying-the-axis-Keyword-Argument-in-Sequential-NumPy-Functions). It will also likely be useful to leverage [boolean indexing](http://www.pythonlikeyoumeanit.com/Module3_IntroducingNumpy/AdvancedIndexing.html#Boolean-Array-Indexing). \n", "\n", "It is strongly recommended that you use matplotlib to visualize your running estimate of $\\pi$. You can refer to [this segment of PLYMI](http://www.pythonlikeyoumeanit.com/Module5_OddsAndEnds/Matplotlib.html#Plotting-and-Saving-a-Figure) to see how to create a simple line plot." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tips\n", "\n", "It is helpful to know about NumPy's [cumulative-sum function](https://numpy.org/doc/stable/reference/generated/numpy.cumsum.html), `numpy.cumsum`. This is useful for keeping a running tally - i.e. keeping the history of the number of darts that have fallen within the circle, and not just the current count." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solution (Unvectorized)\n", "\n", "To start, we want to generate the $(x, y)$ coordinates for our $N = 10,000$ darts. Using `numpy.random.rand(N, 2)`, we can generate a 2D array with $N$ rows - each row contains the $(x, y)$ coordinate for a single dart. \n", "\n", "We want the $x$ and $y$ coordinate of each dart to fall within $[-1, 1]$, respectively. `numpy.random.rand` generates numbers on the interval $[0, 1)$. We can multiply the generated numbers by two and then subtract one to instead generate numbers on the interval $[-1, 1)$." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "N = 10000\n", "dart_positions = 2 * np.random.rand(N, 2) - 1 # generates numbers in [-1, 1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can loop over the dart positions, determine if whether a given dart fell within the circle, and update $N_{circle}$ accordingly." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "Ncircle = [0] # start the count with 0 to make our loop-logic easier\n", "\n", "for x,y in dart_positions: \n", " if np.sqrt(x**2 + y**2) < 1:\n", " Ncircle.append(Ncircle[-1] + 1) # another dart has fallen in the circle\n", " else:\n", " Ncircle.append(Ncircle[-1]) # the dart fell outside of the circle - Ncircle is unchanged" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now lets use our list, `Ncircle`, to compute the running approximation for $\\pi$. We will " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "running_estimate = []\n", "\n", "for n_total, n_circle in enumerate(Ncircle[1:]): # skip the inital 0-count\n", " # n_total starts at 0, so we need to add 1\n", " running_estimate.append(4 * n_circle / (n_total + 1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's print our estimate for our first 10 throws and our last 10 throws. We should see that our estimate is extremely crude and noisy for the first 10 throws and is much more tightly bound for the last 10 throws." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[4.0,\n", " 4.0,\n", " 2.6666666666666665,\n", " 2.0,\n", " 2.4,\n", " 2.6666666666666665,\n", " 2.2857142857142856,\n", " 2.5,\n", " 2.6666666666666665,\n", " 2.8]" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "running_estimate[:10]" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[3.1184065659093183,\n", " 3.1184947958366696,\n", " 3.118583008105674,\n", " 3.118671202721633,\n", " 3.1187593796898447,\n", " 3.118847539015606,\n", " 3.1189356807042112,\n", " 3.119023804760952,\n", " 3.119111911191119,\n", " 3.1192]" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "running_estimate[-10:]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solution (Vectorized)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To start, we want to generate the $(x, y)$ coordinates for our $N = 10,000$ darts. Using `numpy.random.rand(N, 2)`, we can generate a 2D array with $N$ rows - each row contains the $(x, y)$ coordinate for a single dart. \n", "\n", "We want the $x$ and $y$ coordinate of each dart to fall within $[-1, 1]$, respectively. `numpy.random.rand` generates numbers on the interval $[0, 1)$. We can multiply the generated numbers by two and then subtract one to instead generate numbers on the interval $[-1, 1)$." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "N = 10000\n", "dart_positions = 2 * np.random.rand(N, 2) - 1 # generates numbers in [-1, 1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we want to compute the distance from the origin, $\\sqrt{x^2 + y^2}$, for every dart in `dart_positions`. We can do this by squaring each element of `dart_positions` and summing across its *columns* (axis-1), and then taking the square root of the result. This will produce a shape-$(N,)$ array that stores the distance of each dart from the origin." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "dist_from_origin = np.sqrt((dart_positions**2).sum(axis=1)) # shape-(N,) array" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "More concisely, you can use the built-in function `np.linalg.norm` to the same effect." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we want to determine which of those darts fall within the circle. That is, find where $\\sqrt{x^2 + y^2} < 1$. We can simply use `<` to perform an elementwise comparison. This will produce a *boolean-valued* array whose value is `True` for each dart that fell within the circle and `False` for those that didn't." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "is_in_circle = dist_from_origin < 1 # shape-(N,) boolean array" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we want to compute the total number of darts that had landed within the circle at each given point in the sequence of throws. [Recall that](https://www.pythonlikeyoumeanit.com/Module2_EssentialsOfPython/Basic_Objects.html#Boolean-Objects-are-Integers) `True` behaves like `1` and `False` behaves like `0`. Thus we can perform a cumulative sum on `in_circle` to perform this tally." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# cumulative sum: num_in_circle[i] = sum(is_in_circle[:i])\n", "num_in_circle = np.cumsum(is_in_circle)\n", "\n", "num_thrown = np.arange(1, N+1) # 1, 2, ..., N" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally we can compute our approximated value of $\\pi$ for each dart thrown via $N_{circle} / N_{total} \\approx \\pi / 4$." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "running_estimate = 4 * num_in_circle / num_thrown" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's inspect our results by plotting them. We will create a simple line plot, and will include the actual value of pi as a dashed horizontal line. Because we are throwing so many darts, it is useful to plot the number of darts thrown on a log-scale. This will let us to see how our approximation improves on the scale of throwing tens of darts versus hundred versus thousands, etc." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support.' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " if (mpl.ratio != 1) {\n", " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", " }\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var backingStore = this.context.backingStorePixelRatio ||\n", "\tthis.context.webkitBackingStorePixelRatio ||\n", "\tthis.context.mozBackingStorePixelRatio ||\n", "\tthis.context.msBackingStorePixelRatio ||\n", "\tthis.context.oBackingStorePixelRatio ||\n", "\tthis.context.backingStorePixelRatio || 1;\n", "\n", " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width * mpl.ratio);\n", " canvas.attr('height', height * mpl.ratio);\n", " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
')\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('');\n", " button.click(method_name, toolbar_event);\n", " button.mouseover(tooltip, toolbar_mouse_event);\n", " nav_element.append(button);\n", " }\n", "\n", " // Add the status bar.\n", " var status_bar = $('');\n", " nav_element.append(status_bar);\n", " this.message = status_bar[0];\n", "\n", " // Add the close button to the window.\n", " var buttongrp = $('
');\n", " var button = $('');\n", " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", " buttongrp.append(button);\n", " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", " titlebar.prepend(buttongrp);\n", "}\n", "\n", "mpl.figure.prototype._root_extra_style = function(el){\n", " var fig = this\n", " el.on(\"remove\", function(){\n", "\tfig.close_ws(fig, {});\n", " });\n", "}\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(el){\n", " // this is important to make the div 'focusable\n", " el.attr('tabindex', 0)\n", " // reach out to IPython and tell the keyboard manager to turn it's self\n", " // off when our div gets focus\n", "\n", " // location in version 3\n", " if (IPython.notebook.keyboard_manager) {\n", " IPython.notebook.keyboard_manager.register_events(el);\n", " }\n", " else {\n", " // location in version 2\n", " IPython.keyboard_manager.register_events(el);\n", " }\n", "\n", "}\n", "\n", "mpl.figure.prototype._key_event_extra = function(event, name) {\n", " var manager = IPython.notebook.keyboard_manager;\n", " if (!manager)\n", " manager = IPython.keyboard_manager;\n", "\n", " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", " event.shiftKey = false;\n", " // Send a \"J\" for go to next cell\n", " event.which = 74;\n", " event.keyCode = 74;\n", " manager.command_mode();\n", " manager.handle_keydown(event);\n", " }\n", "}\n", "\n", "mpl.figure.prototype.handle_save = function(fig, msg) {\n", " fig.ondownload(fig, null);\n", "}\n", "\n", "\n", "mpl.find_output_cell = function(html_output) {\n", " // Return the cell and output element which can be found *uniquely* in the notebook.\n", " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", " // IPython event is triggered only after the cells have been serialised, which for\n", " // our purposes (turning an active figure into a static one), is too late.\n", " var cells = IPython.notebook.get_cells();\n", " var ncells = cells.length;\n", " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", " data = data.data;\n", " }\n", " if (data['text/html'] == html_output) {\n", " return [cell, data, j];\n", " }\n", " }\n", " }\n", " }\n", "}\n", "\n", "// Register the function which deals with the matplotlib target/channel.\n", "// The kernel may be null if the page has been refreshed.\n", "if (IPython.notebook.kernel != null) {\n", " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", "}\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "%matplotlib notebook\n", "\n", "fig, ax = plt.subplots()\n", "\n", "ax.plot(num_thrown, mean_in_circled, label=\"mean\");\n", "ax.fill_between(num_thrown, y1=mean_in_circled-std_in_circle, y2=mean_in_circled+std_in_circle, \n", " alpha=0.2, label=\"standard deviation\")\n", "ax.hlines(y=np.pi, xmin=1, xmax=N+1, linestyles=\"--\")\n", "\n", "ax.set_xscale(\"log\")\n", "ax.grid(True)\n", "ax.set_ylabel(\"Estimated value of pi\")\n", "ax.set_xlabel(\"Number of darts thrown\")\n", "ax.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we deduced earlier, our approximation reliably improves as we throw more and more darts (as is expected). The standard deviation gives us a nice estimate of how many darts should be thrown in order to reach a given degree of accuracy in our approximation.\n", "\n", "Hopefully you enjoyed this fun hypothetical experiment, and are proud to have run fully-vectorized numerical simulations in NumPy!" ] } ], "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.2" } }, "nbformat": 4, "nbformat_minor": 2 }