{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Monte Carlo determination of $\\pi$\n", "\n", "Imagine that we have a circle of radius 1 ($r = 1$), centered at the point $(0, 0)$ in Cartesian coordinates.\n", "\n", "We can further imagine that circle being enclosed by a square with the same center and with sides of length $2r$, in this case 2.\n", "\n", "\"Drawing\"\n", "\n", "\n", "Given a square with sides of length $2r$, the square's area is\n", "\n", "$$A_{square} = (2r)^2 = 4r^2$$\n", "\n", "whereas the area of a circle with radius $r$ is\n", "$$A_{circle} = \\pi r^2$$\n", "\n", "Therefore the ratio of the area of the circle to that of the enclosing square is\n", "\n", "$$\\frac{A_{circle}}{A_{square}} = \\frac{\\pi r^2}{4r^2} = \\frac{\\pi}{4}$$.\n", "\n", "By rearranging this expression, we see we can define $\\pi$ as\n", "\n", "$$\\pi = 4\\frac{A_{circle}}{A_{square}}$$\n", "\n", "**And *this* suggests a way to calculate $\\pi$!**\n", "\n", "If we're able to figure out the ratio of the areas of a circle and the smallest square that can hold that circle, we can figure out $\\pi$. And fortunately, this ratio is equal to the probability that a randomly chosen point inside the square also falls inside the circle. \n", "\n", "This means that we can calculate this ratio (and therefore $\\pi$) using a Monte Carlo simulation by choosing a large number of random points inside a square and tracking how many fall inside the circle.\n", "\n", "#### Pseudo-code\n", "\n", "Given the above, our algorithm for determining $\\pi$ looks like this:\n", "\n", "1. For each of $N$ iterations,\n", " 1. Select a random point inside the square depicted above.\n", " 1. Determine if the point also falls inside the largest circle that can be enclosed within this square.\n", " 1. Keep track of whether or not this point fell inside the circle. At the end of $N$ iterations, you want to know $M$ -- the number of the $N$ random points that fell inside the circle!\n", "1. Calculate $\\pi$ as $4\\frac{M}{N}$\n", "\n", "Note: As you increase the value of $N$, you'll be able to estimate pi with greater accuracy!\n", "\n", "### Exercise\n", "\n", "Write a function that calculates $\\pi$ using Julia.\n", "\n", "Feel free to get started immediately, or to make use of the practice exercises below which have you solve pieces of the problem one at a time." ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "dist (generic function with 1 method)" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function dist(x, y)\n", " return sqrt(x^2 + y^2)\n", "end" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "myrand (generic function with 1 method)" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function myrand()\n", " return rand([-1, 1]) * rand(), rand([-1, 1]) * rand() \n", "end" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "inside_circle (generic function with 1 method)" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function inside_circle(x, y)\n", " if dist(x, y) < 1.0\n", " return 1\n", " else\n", " return 0\n", " end\n", "end" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "calc_pi" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"\"\"\n", " calc_pi(N)\n", " \n", "Calculate pi using a Monte Carlo simulation in which `N` sampling points are chosen within a square.\n", "\"\"\"\n", "function calc_pi(N)\n", " points_inside = 0\n", " for i in 1:N\n", " x, y = myrand()\n", " points_inside += inside_circle(x, y)\n", " end\n", " return 4 * (points_inside / N)\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Function practice\n", "\n", "#### Exercise\n", "\n", "Given variables `x` and `y` defining a point $(x, y)$, write a function `dist(x, y)` that returns the distance of $(x, y)$ from $(0, 0)$." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "dist (generic function with 1 method)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function dist(x, y)\n", " return sqrt(x^2 + y^2)\n", "end" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0.09177130543526113, 0.880104215610646)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x, y = rand(), rand()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.8848759251086687" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dist(x, y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Exercise\n", "\n", "Write a function `myrand` that returns `(x, y)`, where `x` and `y` are pseudo-random numbers uniformly distributed between -1 and +1.\n", "\n", "*Hint:* The `rand()` function returns a random number between 0 and 1. Passing a collection as the first argument of `rand` will cause `rand` to return a random value from amongst the elements of the collection.\n", "\n", "Ex. rand(animals) will return one of the strings contained in `animals`:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"cat\"" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "animals = [\"cat\", \"dog\", \"bird\"]\n", "rand(animals)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "myrand (generic function with 1 method)" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function myrand()\n", " return rand([-1, 1]) * rand(), rand([-1, 1]) * rand() \n", "end" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(-0.9230968519773852, -0.38550538642272825)" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "myrand()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Conditional practice\n", "\n", "#### Exercise\n", "\n", "Given variables `x` and `y` defining a point $(x, y)$, write a conditional statement that prints \"Less than 1!\" if the distance of $(x, y)$ from $(0, 0)$ is less than 1." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Array{Float64,1}:\n", " 0.2753495990422792 \n", " 0.011453969623016924" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x, y = rand(2)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Less than 1!\n" ] } ], "source": [ "if sqrt(x^2 + y^2) < 1.0\n", " println(\"Less than 1!\")\n", "end" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Less than 1!\n" ] } ], "source": [ "if dist(x, y) < 1.0\n", " println(\"Less than 1!\")\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Exercise\n", "\n", "Write a function `inside_circle(x, y)` that takes coordinates/variables `x` and `y` and determines whether the point $(x, y)$ falls inside a circle of radius 1 that is centered at $(0, 0)$. Return `1` if the point falls within the circle and `0` if not.\n", "\n", "*Hint:* To fall within the circle, the distance of the point to the circle's center should be less than the radius." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "inside_circle (generic function with 1 method)" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function inside_circle(x, y)\n", " if sqrt(x^2 + y^2) < 1.0\n", " return 1\n", " else\n", " return 0\n", " end\n", "end" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "inside_circle (generic function with 1 method)" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "inside_circle(x, y) = (sqrt(x^2 + y^2) < 1.0) ? 1 : 0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Loop practice" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Exercise\n", "\n", "Use a loop to generate 100 random points $(x_i, y_i)$ where $-1.0 <= x <= 1.0$ and $-1.0 <= y <= 1.0$. Print the number of times the points fell within a circle of radius 1, centered at $(0, 0)$." ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "83 points fell inside the circle!\n" ] } ], "source": [ "N = 0\n", "for i in 1:100\n", " x, y = myrand()\n", " N += inside_circle(x, y)\n", "end\n", "println(\"$N points fell inside the circle!\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Without work that we've already done in previous exercises:" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "78 points fell inside the circle!\n" ] } ], "source": [ "N = 0\n", "for i in 1:100\n", " x = rand([-1, 1]) * rand()\n", " y = rand([-1, 1]) * rand()\n", " if sqrt(x^2 + y^2) < 1.0\n", " N += 1\n", " end\n", "end\n", "println(\"$N points fell inside the circle!\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Array practice" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Exercise\n", "\n", "Use an array comprehension to create an array of 100 points in 2D Cartesian space." ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "100-element Array{Tuple{Float64,Float64},1}:\n", " (0.8235803117901142, -0.6736936601474774) \n", " (0.3558640960890065, -0.7080181200055458) \n", " (0.6161108313295864, 0.10128327704482909) \n", " (0.06767258090743411, 0.48567762731410147) \n", " (-0.01869616622547987, -0.4229478969051781)\n", " (0.8599279235169288, -0.20686834775240293) \n", " (0.25769489684403424, -0.9587876421300074) \n", " (-0.7108573976816956, -0.7344179861513818) \n", " (0.7646953306782522, 0.33019531478167496) \n", " (0.20513540415981146, -0.9106607051488109) \n", " (-0.15357265925894148, -0.4578539987702732)\n", " (0.76617200964695, -0.7296079522154135) \n", " (-0.5640083749041798, -0.8156093136312135) \n", " ⋮ \n", " (0.4825477312267852, -0.6023767911672253) \n", " (0.20483110741326693, 0.16210952105073861) \n", " (-0.3645209544985353, -0.46411948118371793)\n", " (0.1762587416098942, -0.876202029020815) \n", " (-0.4087320814383706, 0.3464881102351882) \n", " (0.803515030181229, 0.227084405652221) \n", " (-0.14559086789114306, -0.543519582178885) \n", " (0.770377554719148, -0.5780172129728407) \n", " (0.9773512106749624, 0.7160989375188618) \n", " (0.5643735149375608, -0.9564993822734364) \n", " (0.31253625667570106, 0.9430302711198242) \n", " (-0.7897705751424646, -0.8060195845846594) " ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[myrand() for i in 1:100]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Exercise\n", "\n", "Use an array comprehension to create an array of the distances of 100 points in 2D Cartesian space from a point at $(0, 0)$.\n", "\n", "*Hint:* You may find \"splat\", i.e. `...`, useful here. Adding `...` after a collection passed as an input argument to a function \"opens up\" that collection, passing each of the elements of the collection as inputs to the function. For example, if\n", "\n", "```\n", "a = [1, 2, 3]\n", "add_three_things(x, y, z) = x + y + z\n", "```\n", "then `add_three_things(a)` will give an error, but `add_three_things(a...)` will yield `6`, just as if you had run `add_three_things(a[1], a[2], a[3])`" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = [1, 2, 3]\n", "add_three_things(x, y, z) = x + y + z\n", "add_three_things(a...)" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "100-element Array{Float64,1}:\n", " 0.9383872456685842 \n", " 0.5839277168087241 \n", " 0.7355502495040855 \n", " 0.8441859785114507 \n", " 0.6807536016582622 \n", " 0.7632412130157962 \n", " 0.7317234224262192 \n", " 0.8824052199190965 \n", " 0.7755182140496781 \n", " 0.9744510181703258 \n", " 0.25005959593047916\n", " 1.0612231389451416 \n", " 1.0606971366415012 \n", " ⋮ \n", " 0.9484156186801598 \n", " 0.9270152199659667 \n", " 0.6614786312111012 \n", " 0.6478811880980543 \n", " 0.23342156646680173\n", " 1.1057076816888292 \n", " 0.5958913950161429 \n", " 1.2529131820213943 \n", " 0.7256403533033896 \n", " 0.910877586092377 \n", " 0.779338060314894 \n", " 0.16368015145789885" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[dist(myrand()...) for i in 1:100]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.0.1", "language": "julia", "name": "julia-1.0" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.0.1" } }, "nbformat": 4, "nbformat_minor": 2 }