HSG-MCS-HS21_Julia/Problemsets/MC_pi.ipynb

312 lines
95 KiB
Plaintext
Raw Normal View History

2021-11-15 20:14:51 +00:00
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"## Load Packages"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2020-09-25T11:59:02.029Z",
"iopub.status.busy": "2020-09-25T11:59:01.320Z",
"iopub.status.idle": "2020-09-25T11:59:08.538Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"printyellow (generic function with 1 method)"
]
},
"execution_count": 1,
"metadata": {
},
"output_type": "execute_result"
}
],
"source": [
"using Printf, Distributions\n",
"include(\"jlFiles/printmat.jl\")"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2020-09-25T11:59:08.640Z",
"iopub.status.busy": "2020-09-25T11:59:08.559Z",
"iopub.status.idle": "2020-09-25T11:59:25.236Z"
}
},
"outputs": [
],
"source": [
"using Plots\n",
"gr(size=(480,480))\n",
"default(fmt = :svg)"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"# Find $\\pi$\n",
"\n",
"Find $\\pi$ by checking how often a pair of random variables are inside a circle. \n",
"\n",
"Clearly, this is not the smartest way to find $pi$, but the steps require some useful skills. Who knows, maybe it's smarter to remember how many letters there are in each of the words in \"How I wish I could calculate pi\" (3.141592).\n",
"\n",
"The next few cells will give you a useful starting point."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2020-09-25T11:59:40.900Z",
"iopub.status.busy": "2020-09-25T11:59:40.890Z",
"iopub.status.idle": "2020-09-25T11:59:42.059Z"
}
},
"outputs": [
{
"data": {
"image/svg+xml": "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"480\" height=\"480\" viewBox=\"0 0 1920 1920\">\n<defs>\n <clipPath id=\"clip240\">\n <rect x=\"0\" y=\"0\" width=\"1920\" height=\"1920\"/>\n </clipPath>\n</defs>\n<path clip-path=\"url(#clip240)\" d=\"\nM0 1920 L1920 1920 L1920 0 L0 0 Z\n \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n<defs>\n <clipPath id=\"clip241\">\n <rect x=\"384\" y=\"192\" width=\"1345\" height=\"1345\"/>\n </clipPath>\n</defs>\n<path clip-path=\"url(#clip240)\" d=\"\nM189.136 1785.97 L1872.76 1785.97 L1872.76 153.712 L189.136 153.712 Z\n \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n<defs>\n <clipPath id=\"clip242\">\n <rect x=\"189\" y=\"153\" width=\"1685\" height=\"1633\"/>\n </clipPath>\n</defs>\n<polyline clip-path=\"url(#clip242)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 236.786,1785.97 236.786,153.712 \n \"/>\n<polyline clip-path=\"url(#clip242)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 633.866,1785.97 633.866,153.712 \n \"/>\n<polyline clip-path=\"url(#clip242)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 1030.95,1785.97 1030.95,153.712 \n \"/>\n<polyline clip-path=\"url(#clip242)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 1428.03,1785.97 1428.03,153.712 \n \"/>\n<polyline clip-path=\"url(#clip242)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 1825.11,1785.97 1825.11,153.712 \n \"/>\n<polyline clip-path=\"url(#clip240)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 189.136,1785.97 1872.76,1785.97 \n \"/>\n<polyline clip-path=\"url(#clip240)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 236.786,1785.97 236.786,1767.07 \n \"/>\n<polyline clip-path=\"url(#clip240)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 633.866,1785.97 633.866,1767.07 \n \"/>\n<polyline clip-path=\"url(#clip240)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 1030.95,1785.97 1030.95,1767.07 \n \"/>\n<polyline clip-path=\"url(#clip240)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 1428.03,1785.97 1428.03,1767.07 \n \"/>\n<polyline clip-path=\"url(#clip240)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 1825.11,1785.97 1825.11,1767.07 \n \"/>\n<path clip-path=\"url(#clip240)\" d=\"M183.314 1825.97 L212.989 1825.97 L212.989 1829.91 L183.314 1829.91 L183.314 1825.97 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip240)\" d=\"M223.892 1838.86 L231.531 1838.86 L231.531 1812.5 L223.221 1814.17 L223.221 1809.91 L231.485 1808.24 L236.161 1808.24 L236.161 1838.86 L243.799 1838.86 L243.799 1842.8 L223.892 1842.8 L223.892 1838.86 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip240)\" d=\"M253.244 1836.92 L258.128 1836.92 L258.128 1842.8 L253.244 1842.8 L253.244 1836.92 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip240)\" d=\"M278.313 1811.32 Q274.702 1811.32 272.873 1814.88 Q271.068 1818.43 271.068 1825.55 Q271.068 1832.66 272.873 1836.23 Q274.702 1839.77 278.313 1839.77 Q281.947 1839.77 283.753 1836.23 Q285.582 1
},
"execution_count": 3,
"metadata": {
},
"output_type": "execute_result"
}
],
"source": [
"θ = range(0,2*pi,length=201) #angles (in radians) from 0 to 2*pi\n",
"\n",
"xCirc = cos.(θ)\n",
"yCirc = sin.(θ)\n",
"\n",
"p1 = plot( xCirc,yCirc,\n",
" legend=nothing,\n",
" title = \"A circle with radius of 1 and center at (0,0)\" )\n",
"display(p1)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" x y x.^2 + y.^2 v\n",
" 0.492 0.961 1.166 0 \n",
" -0.118 -0.275 0.089 1 \n",
" -0.227 0.742 0.602 1 \n",
" -0.558 -0.404 0.475 1 \n",
" 0.183 0.431 0.219 1 \n",
" -0.769 -0.847 1.309 0 \n",
" 0.526 0.731 0.811 1 \n",
" 0.777 0.530 0.885 1 \n",
" 0.427 -0.490 0.423 1 \n",
" -0.940 0.348 1.005 0 \n",
" 0.999 0.006 0.999 1 \n",
" -0.833 0.594 1.046 0 \n",
" -0.885 0.907 1.606 0 \n",
" 0.665 0.225 0.492 1 \n",
" -0.326 -0.178 0.138 1 \n",
" 0.875 -0.778 1.372 0 \n",
" 0.688 -0.440 0.667 1 \n",
" 0.630 0.141 0.416 1 \n",
" 0.092 0.742 0.559 1 \n",
" -0.447 -0.630 0.596 1 \n",
" -0.596 0.589 0.703 1 \n",
" 0.532 0.403 0.446 1 \n",
" -0.805 -0.740 1.195 0 \n",
" -0.017 -0.562 0.316 1 \n",
" -0.899 -0.723 1.332 0 \n",
"\n"
]
}
],
"source": [
"N = 25\n",
"\n",
"dist = Uniform(-1,1)\n",
"x = rand(dist,N) #draw uniformly [-1,1] distributed numbers\n",
"y = rand(dist,N)\n",
"\n",
"v = (x.^2 + y.^2) .<= 1 #true if inside\n",
"\n",
"printmat(x,y,x.^2 + y.^2,v,colNames=[\"x\",\"y\",\"x.^2 + y.^2\",\"v\"],width=14)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2020-09-25T11:59:43.020Z",
"iopub.status.busy": "2020-09-25T11:59:43.009Z",
"iopub.status.idle": "2020-09-25T11:59:43.481Z"
}
},
"outputs": [
{
"data": {
"image/svg+xml": "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"480\" height=\"480\" viewBox=\"0 0 1920 1920\">\n<defs>\n <clipPath id=\"clip280\">\n <rect x=\"0\" y=\"0\" width=\"1920\" height=\"1920\"/>\n </clipPath>\n</defs>\n<path clip-path=\"url(#clip280)\" d=\"\nM0 1920 L1920 1920 L1920 0 L0 0 Z\n \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n<defs>\n <clipPath id=\"clip281\">\n <rect x=\"384\" y=\"192\" width=\"1345\" height=\"1345\"/>\n </clipPath>\n</defs>\n<path clip-path=\"url(#clip280)\" d=\"\nM276.164 1698.94 L1872.76 1698.94 L1872.76 153.712 L276.164 153.712 Z\n \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n<defs>\n <clipPath id=\"clip282\">\n <rect x=\"276\" y=\"153\" width=\"1598\" height=\"1546\"/>\n </clipPath>\n</defs>\n<polyline clip-path=\"url(#clip282)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 321.351,1698.94 321.351,153.712 \n \"/>\n<polyline clip-path=\"url(#clip282)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 697.905,1698.94 697.905,153.712 \n \"/>\n<polyline clip-path=\"url(#clip282)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 1074.46,1698.94 1074.46,153.712 \n \"/>\n<polyline clip-path=\"url(#clip282)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 1451.01,1698.94 1451.01,153.712 \n \"/>\n<polyline clip-path=\"url(#clip282)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n 1827.57,1698.94 1827.57,153.712 \n \"/>\n<polyline clip-path=\"url(#clip280)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 276.164,1698.94 1872.76,1698.94 \n \"/>\n<polyline clip-path=\"url(#clip280)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 321.351,1698.94 321.351,1680.04 \n \"/>\n<polyline clip-path=\"url(#clip280)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 697.905,1698.94 697.905,1680.04 \n \"/>\n<polyline clip-path=\"url(#clip280)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 1074.46,1698.94 1074.46,1680.04 \n \"/>\n<polyline clip-path=\"url(#clip280)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 1451.01,1698.94 1451.01,1680.04 \n \"/>\n<polyline clip-path=\"url(#clip280)\" style=\"stroke:#000000; stroke-linecap:butt; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n 1827.57,1698.94 1827.57,1680.04 \n \"/>\n<path clip-path=\"url(#clip280)\" d=\"M267.879 1738.94 L297.554 1738.94 L297.554 1742.88 L267.879 1742.88 L267.879 1738.94 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip280)\" d=\"M308.457 1751.84 L316.096 1751.84 L316.096 1725.47 L307.786 1727.14 L307.786 1722.88 L316.05 1721.21 L320.726 1721.21 L320.726 1751.84 L328.364 1751.84 L328.364 1755.77 L308.457 1755.77 L308.457 1751.84 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip280)\" d=\"M337.809 1749.89 L342.693 1749.89 L342.693 1755.77 L337.809 1755.77 L337.809 1749.89 Z\" fill=\"#000000\" fill-rule=\"evenodd\" fill-opacity=\"1\" /><path clip-path=\"url(#clip280)\" d=\"M362.878 1724.29 Q359.267 1724.29 357.438 1727.86 Q355.633 1731.4 355.633 1738.53 Q355.633 1745.63 357.438 1749.2 Q359.267 1752.74 362.878 1752.74 Q366.512 1752.74 368.318 1749.2 Q370.147
},
"execution_count": 5,
"metadata": {
},
"output_type": "execute_result"
}
],
"source": [
"p1 = plot( xCirc,yCirc,\n",
" legend = nothing,\n",
" title = \"What fraction is inside the circle?\",\n",
" xlabel = \"x\",\n",
" ylabel = \"y\" )\n",
"scatter!(x[v],y[v],color=:blue)\n",
"scatter!(x[.!v],y[.!v],color=:red)\n",
"display(p1)"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"## Hints\n",
"\n",
"Draw many more random numbers (1_000_000 or more) and calculate $pi$ as the fraction of points inside the circle times the area of the box.\n",
"\n",
"The area of the box is $2x2=4$, and a pair is inside or on the circle if $x^2 + y^2 \\le 1$"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2020-09-25T11:59:43.499Z",
"iopub.status.busy": "2020-09-25T11:59:43.491Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.785 3.141 π\n",
"\n"
]
}
],
"source": [
"N = 25_000_000\n",
"\n",
"x = rand(dist,N)\n",
"y = rand(dist,N)\n",
"\n",
"FreqInCircle = mean((x.^2 + y.^2) .<= 1)\n",
"\n",
"piEstimate1 = FreqInCircle*4\n",
"\n",
"printmat(FreqInCircle,piEstimate1,pi)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"false"
]
},
"execution_count": 8,
"metadata": {
},
"output_type": "execute_result"
}
],
"source": [
"1 > 2"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"collapsed": false
},
"outputs": [
],
"source": [
]
}
],
"metadata": {
"@webio": {
"lastCommId": null,
"lastKernelId": null
},
"anaconda-cloud": {
},
"kernelspec": {
"display_name": "Julia 1.6",
"env": {
"JULIA_DEPOT_PATH": "/home/user/.julia/:/ext/julia/depot/",
"JULIA_PROJECT": "/home/user/.julia/environment/v1.6"
},
"language": "julia",
"metadata": {
"cocalc": {
"description": "The Julia Programming Language",
"priority": 10,
"url": "https://julialang.org/"
}
},
"name": "julia-1.6",
"resource_dir": "/ext/jupyter/kernels/julia-1.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.6.2"
},
"nteract": {
"version": "0.25.0"
}
},
"nbformat": 4,
"nbformat_minor": 4
}