1. An Introductory Example#
1.1. Overview#
ഇനി നമുക്ക് Python language-നെ പറ്റി പഠിക്കാം.
ഈ lecture-ൽ, നമ്മൾ ചെറിയ Python programs എഴുതി, അവയെ വിശദമായി പരിശോധിക്കും.
Basic Python syntax-ഉം data structures-ഉം നിങ്ങൾക്ക് introduce ചെയ്തു തരിക എന്നതാണ് ഈ lecture-ന്റെ ലക്ഷ്യം.
കൂടുതൽ ആഴത്തിലുള്ള concepts പിന്നീടുള്ള lectures-ൽ cover ചെയ്യും.
ഈ lecture-ലേക്കു കടക്കുന്നതിനു മുമ്പ് നിങ്ങൾ getting started with Python എന്ന lecture വായിച്ചിട്ടുണ്ടാകും എന്ന് കരുതുന്നു.
1.2. The Task: Plotting a White Noise Process#
നമുക്ക് ഒരു white noise process (\(\epsilon_0, \epsilon_1, \ldots, \epsilon_T\)) simulate ചെയ്ത്, plot ചെയ്യണം എന്ന് കരുതുക — ഓരോ draw-ഉം (\(\epsilon_t\)) ഒരു independent standard normal ആണ്.
അതായത്, താഴെ കൊടുത്തിരിക്കുന്നത് പോലെയുള്ള ഒരു figure നമുക്ക് generate ചെയ്യണം:
(ഇവിടെ \(t\) horizontal axis-ലും, \(\epsilon_t\) vertical axis-ലും ആണ്.)
ഈ white noise process plotting നമ്മൾ പല വിധത്തിൽ ചെയ്യും — ഓരോ രീതിയിൽ ചെയ്യുമ്പോഴും നമ്മൾ Python-നെക്കുറിച്ച് കൂടുതൽ കാര്യങ്ങൾ പഠിക്കും.
1.3. Version 1#
നമ്മൾ set ചെയ്ത task ചെയ്യുന്ന കുറച്ച് lines of code താഴെ കാണാം.
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng()
ϵ_values = rng.standard_normal(100)
plt.plot(ϵ_values)
plt.show()
നമുക്ക് ഈ program-നെ break down ചെയ്ത്, അത് എങ്ങനെ work ചെയ്യുന്നു എന്ന് നോക്കാം.
1.3.1. Imports#
ഈ program-ന്റെ ആദ്യത്തെ രണ്ട് lines, external code libraries-ൽ നിന്നും functionality import ചെയ്യുന്നു.
ആദ്യത്തെ line NumPy import ചെയ്യുന്നു — NumPy, താഴെ കൊടുത്തിരിക്കുന്ന പോലത്തെ tasks-കൾ ചെയ്യാനുള്ള ഒരു Python package ആണ്:
working with arrays (vectors and matrices)
common mathematical functions like
cosandsqrtgenerating random numbers
linear algebra, etc.
import numpy as np ചെയ്ത് കഴിഞ്ഞാൽ, np.attribute എന്ന syntax ഉപയോഗിച്ച് ഈ attributes നമുക്ക് access ചെയ്യാം.
Here's two more examples
np.sqrt(4)
np.float64(2.0)
np.log(4)
np.float64(1.3862943611198906)
1.3.1.1. Why So Many Imports?#
Python programs-ൽ സാധാരണയായി multiple import statements ആവശ്യമായിവരുന്നു.
കാരണം, core language മനഃപൂർവ്വം ചെറുതായി നിലനിർത്തിയിരിക്കുന്നു, അതുകൊണ്ട് അത് പഠിക്കാനും maintain ചെയ്യാനും improve ചെയ്യാനും easy ആണ്.
Python ഉപയോഗിച്ച് interesting ആയ എന്തെങ്കിലും ചെയ്യണമെങ്കിൽ, മിക്ക സമയത്തും additional functionality import ചെയ്യേണ്ടി വരും.
1.3.1.2. Packages#
മുകളിൽ പറഞ്ഞതുപോലെ, NumPy ഒരു Python package ആണ്.
Share ചെയ്യാൻ ആഗ്രഹിക്കുന്ന code-നെ organize ചെയ്യാനാണ് developers packages ഉപയോഗിക്കുന്നത്.
In fact, ഒരു package എന്നത് താഴെ കൊടുത്തിരിക്കുന്നവ അടങ്ങിയ ഒരു directory ആണ്:
Python code ഉള്ള files — Python-ന്റെ ഭാഷയിൽ ഇവയെ modules എന്ന് വിളിക്കുന്നു
Python-ന് access ചെയ്യാൻ കഴിയുന്ന compiled code (e.g., functions compiled from C or FORTRAN code)
__init__.pyഎന്ന ഒരു file — നമ്മൾimport package_nametype ചെയ്യുമ്പോൾ എന്ത് execute ചെയ്യണം എന്ന് ഇത് specify ചെയ്യുന്നു
NumPy-യുടെ __init__.py-യുടെ location check ചെയ്യാൻ, താഴെ കൊടുത്തിരിക്കുന്ന code Python-ൽ run ചെയ്യുക:
import numpy as np
print(np.__file__)
1.3.1.3. Subpackages#
rng = np.random.default_rng() എന്ന line നോക്കുക.
ഇവിടെ np എന്നത് NumPy package-നെ refer ചെയ്യുന്നു, അതേ സമയം random എന്നത് NumPy-യുടെ ഒരു subpackage ആണ്.
Subpackages എന്നത് മറ്റൊരു package-ന്റെ subdirectory ആയ packages മാത്രമാണ്.
ഉദാഹരണത്തിന്, NumPy-യുടെ directory-യിൽ random എന്ന folder കാണാം.
1.3.2. Importing Names Directly#
മുകളിൽ കണ്ട ഈ code ഓർക്കുക
import numpy as np
np.sqrt(4)
np.float64(2.0)
NumPy-യുടെ square root function access ചെയ്യാനുള്ള മറ്റൊരു രീതി താഴെ കാണാം:
from numpy import sqrt
sqrt(4)
np.float64(2.0)
ഇതും correct ആണ്.
ഇതിന്റെ advantage — നമ്മുടെ code-ൽ sqrt പലപ്പോഴും ഉപയോഗിക്കുകയാണെങ്കിൽ, ഇതുവഴി നമുക്ക് typing കുറക്കാൻ സാധിക്കും.
ഇതിന്റെ disadvantage — ഒരു long program-ൽ, ഈ രണ്ട് lines-ന്റെയും ഇടയിൽ മറ്റ് പല lines-ഉം വന്നേക്കാം.
അങ്ങനെ വരുമ്പോൾ, sqrt എവിടെ നിന്ന് വന്നു എന്ന് readers-ന് മനസ്സിലാക്കാൻ ബുദ്ധിമുട്ടായേക്കാം.
1.3.3. Random Draws#
White noise plot ചെയ്യുന്ന നമ്മുടെ program-ലേക്ക് തിരികെ വരാം. Import statements കഴിഞ്ഞുള്ള ബാക്കി മൂന്ന് lines ഇവയാണ്:
ϵ_values = rng.standard_normal(100)
plt.plot(ϵ_values)
plt.show()
ആദ്യത്തെ line, 100 (quasi) independent standard normals-നെ generate ചെയ്ത്, ϵ_values-ൽ store ചെയ്യുന്നു.
അടുത്ത രണ്ട് lines plot generate ചെയ്യുന്നു.
ഈ plot configure ചെയ്യാനും improve ചെയ്യാനുമുള്ള പല വഴികൾ നമുക്ക് താഴെ കാണാം.
1.4. Alternative Implementations#
Standard normal distribution-ൽ നിന്നും IID draws plot ചെയ്ത നമ്മുടെ ആദ്യത്തെ program, മറ്റു പല രീതികളിൽ എങ്ങനെ എഴുതാം എന്ന് നോക്കാം.
താഴെ കൊടുത്തിരിക്കുന്ന programs, original-ന്റെ അത്ര efficient അല്ല - അവ ഒരു ആശയം വിശദീകരിക്കാനായി മാത്രം നിർമ്മിച്ചവയാണ്.
എന്നാൽ ഇവ ഒരു familiar setting-ൽ ചില പ്രധാന Python syntax-ഉം, semantics-ഉം illustrate ചെയ്യാൻ സഹായിക്കുന്നു.
1.4.1. A Version with a For Loop#
for loops-ഉം, Python lists-ഉം illustrate ചെയ്യുന്ന ഒരു version താഴെ കാണാം:
ts_length = 100
ϵ_values = [] # empty list
for i in range(ts_length):
e = rng.standard_normal()
ϵ_values.append(e)
plt.plot(ϵ_values)
plt.show()
ചുരുക്കത്തിൽ,
ആദ്യത്തെ line, time series-ന് ആവശ്യമായ length set ചെയ്യുന്നു.
അടുത്ത line,
ϵ_valuesഎന്നൊരു empty list create ചെയ്യുന്നു — അതിൽ ആയിരിക്കും നമ്മൾ generate ചെയ്യുന്ന \(\epsilon_t\) values store ചെയ്യുക.# empty listഎന്ന statement ഒരു comment ആണ്, Python-ന്റെ interpreter ഇത് ignore ചെയ്യും.അടുത്ത മൂന്ന് lines ആണ്
forloop — ഇത് repeatedly ഒരു പുതിയ random number \(\epsilon_t\) draw ചെയ്ത്ϵ_valueslist-ന്റെ അവസാനം append ചെയ്യുന്നു.അവസാനത്തെ രണ്ട് lines, plot generate ചെയ്ത്, user-ന് display ചെയ്യുന്നു.
ഈ program-ന്റെ കുറച്ചു ഭാഗങ്ങൾ നമുക്ക് വിശദമായി പഠിക്കാം.
1.4.2. Lists#
ϵ_values = [] എന്ന statement നോക്കുക. ഇത് ഒരു empty list create ചെയ്യുന്നു.
ഒരു കൂട്ടം objects-നെ ഒരുമിച്ച് group ചെയ്യാൻ ഉപയോഗിക്കുന്ന Python-ന്റെ ഒരു native data structure ആണ് Lists.
Lists-ലെ items ordered ആണ്, കൂടാതെ lists-ൽ duplicates അനുവദനീയമാണ്.
ഉദാഹരണത്തിന്, ഇത് try ചെയ്യുക
x = [10, 'foo', False]
type(x)
list
x-ന്റെ ആദ്യത്തെ element ഒരു integer ആണ്, അടുത്തത് ഒരു string ആണ്, മൂന്നാമത്തേത് ഒരു Boolean value ആണ്.
ഒരു list-ലേക്ക് ഒരു value add ചെയ്യാൻ, list_name.append(some_value) എന്ന syntax നമുക്ക് ഉപയോഗിക്കാം
x
[10, 'foo', False]
x.append(2.5)
x
[10, 'foo', False, 2.5]
ഇവിടെ append() എന്നത് ഒരു method ആണ്. ഒരു object-നോട് "attach" ആയിരിക്കുന്ന ഒരു function-നെയാണ് method എന്ന് വിളിക്കുന്നത്. ഇവിടെ ആ object x എന്ന list ആണ്.
Methods-നെ പറ്റി നമ്മൾ പിന്നീട് വിശദമായി പഠിക്കും, പക്ഷേ ഇപ്പോൾ ഒരു idea കിട്ടാൻ നിങ്ങൾ ഇത്രെയും മനസിലാക്കുക:
Lists, strings തുടങ്ങിയ Python objects-ന് എല്ലാം, അവയിൽ അടങ്ങിയിരിക്കുന്ന data manipulate ചെയ്യാൻ ഉപയോഗിക്കുന്ന methods ഉണ്ട്.
String objects-ന് string methods ഉണ്ട്, list objects-ന് list methods ഉണ്ട്, അങ്ങനെ ഓരോ object-ഇനും അതിന് suitable ആയ methods ഉണ്ടായിരിക്കും.
മറ്റൊരു useful list method ആണ് pop()
x
[10, 'foo', False, 2.5]
x.pop()
2.5
x
[10, 'foo', False]
Python-ലെ lists zero-based ആണ് (as in C, Java or Go). അതിനാൽ list-ൽ, ആദ്യത്തെ element-ന്റെ reference x[0] ആയി ഉപയോഗിക്കുന്നു
x[0] # first element of x
10
x[1] # second element of x
'foo'
1.4.3. The For Loop#
ഇനി നമുക്ക് മുകളിലുള്ള program-ലെ for loop നോക്കാം. അവിടെ നമ്മൾ ഉപയോഗിച്ച for loop താഴെ കാണാം:
for i in range(ts_length):
e = rng.standard_normal()
ϵ_values.append(e)
Indent ചെയ്തിരിക്കുന്ന ഈ രണ്ട് lines, ts_length തവണ execute ചെയ്ത ശേഷമേ Python മുന്നോട്ട് പോകൂ.
ഈ രണ്ട് lines-നെ നാം code block എന്ന് വിളിക്കുന്നു - കാരണം നമ്മൾ for loop ഉപയോഗിച്ച് ഈ "block" of code-നെയാണ് loop ചെയ്യിക്കുന്നത്.
ഒരു code block എവിടെ വരെ extent ചെയ്യുന്നു എന്ന് Python മനസ്സിലാക്കുന്നത് അതിന്റെ indentation മാത്രം ഉപയോഗിച്ചാണ്. ഇത്, മറ്റു പല programming languages-ൽ നിന്നും Python-നെ വ്യത്യസ്തമാക്കുന്നു.
നമ്മുടെ program-ൽ, ϵ_values.append(e) എന്ന line-ന് ശേഷം indentation കുറയുന്നു. ഇതിലൂടെ, ആ code block അവിടെ അവസാനിക്കുന്നു എന്ന് Python മനസ്സിലാക്കുന്നു.
Indentation-നെ പറ്റി കൂടുതൽ താഴെ കാണാം — ഇപ്പോൾ for loop-ന്റെ മറ്റൊരു example നോക്കാം.
animals = ['dog', 'cat', 'bird']
for animal in animals:
print("The plural of " + animal + " is " + animal + "s")
The plural of dog is dogs
The plural of cat is cats
The plural of bird is birds
ഈ example, for loop എങ്ങനെ പ്രവർത്തിക്കുന്നു എന്ന് clarify ചെയ്യാൻ സഹായിക്കുന്നു: താഴെ കൊടുത്തിരിക്കുന്ന രീതിയിൽ ഒരു loop execute ചെയ്യുമ്പോൾ,
for variable_name in sequence:
<code block>
Python interpreter ഇവ perform ചെയ്യുന്നു:
sequence-ലെ ഓരോ element-ഇനും, Python ആ element-നെvariable_nameഎന്ന name "bind" ചെയ്യുന്നു. തുടർന്ന് code block execute ചെയ്യുന്നു.
1.4.4. A Comment on Indentation#
for loop discuss ചെയ്തപ്പോൾ, loop ചെയ്യപ്പെടുന്ന code block-ന്റെ delimit, അതിന്റെ indentation ഉപയോഗിച്ചാണ് Python മനസ്സിലാക്കുന്നതെന്ന് നമ്മൾ discuss ചെയ്തിരുന്നു.
In fact, Python-ൽ, എല്ലാ code blocks-ഉം (അതായത്, loops-ന്റെ code block, if clauses-ന്റെ code block, function definitions-ന്റെ code block, etc.) indentation ഉപയോഗിച്ചാണ് delimit ചെയ്യപ്പെടുന്നത്.
അതിനാൽ, മറ്റു മിക്ക programming languages-ൽ നിന്നും വ്യത്യസ്തമായി, Python code-ലെ whitespace, program-ന്റെ output-നെ affect ചെയ്യുന്നു.
ഒരിക്കൽ ഇത് ശീലമായാൽ, ഇത് ഒരു നല്ല കാര്യമാണ്.
clean-ഉം consistent-ഉം ആയ indentation വഴി readability improve ചെയ്യുന്നു
മറ്റ് languages-ൽ ഉപയോഗിക്കുന്ന brackets അല്ലെങ്കിൽ end statements പോലുള്ള clutter remove ചെയ്യുന്നു
On the other hand, ഇത് correct ആയി ഉപയോഗിക്കാൻ ഒരല്പം care ആവശ്യമാണ്. അതിനാൽ താഴെപ്പറയുന്ന കാര്യങ്ങൾ ഓർത്തിരിക്കുക:
ഒരു code block ആരംഭിക്കുന്നതിന് മുമ്പുള്ള line എപ്പോഴും colon-ൽ അവസാനിക്കണം
for i in range(10):if x > y:while x < 100:etc.
ഒരു code block-ലെ എല്ലാ lines-ഇനും ഒരേ amount of indentation ഉണ്ടായിരിക്കണം.
Python-ന്റെ standard 4 spaces ആണ്. അതിനാൽ നിങ്ങളും 4 spaces ഉപയോഗിക്കണം.
1.4.5. While Loops#
Python-ൽ iteration ചെയ്യാൻ ഏറ്റവും common ആയി ഉപയോഗിക്കുന്ന technique ആണ് for loop.
എന്നാൽ, illustration purpose-ന് വേണ്ടി, മുൻപത്തെ program-ൽ for loop-നു പകരം while loop ഉപയോഗിച്ച് എങ്ങനെ ചെയ്യാം എന്ന് നോക്കാം.
ts_length = 100
ϵ_values = []
i = 0
while i < ts_length:
e = rng.standard_normal()
ϵ_values.append(e)
i = i + 1
plt.plot(ϵ_values)
plt.show()
Indentation ഉപയോഗിച്ച് delimit ചെയ്തിരിക്കുന്ന while loop-ന്റെ code block, (i < ts_length) എന്ന condition satisfy ആകുന്നത് വരെ execute ചെയ്ത് കൊണ്ടേയിരിക്കും.
ഈ case-ൽ, i ts_length-ന് equal ആകുന്നത് വരെ program ϵ_values list-ലേക്ക് values add ചെയ്ത് കൊണ്ടേയിരിക്കും:
i == ts_length #the ending condition for the while loop
True
ശ്രദ്ധിക്കുക,
whileloop-ന്റെ code block, indentation മാത്രം ഉപയോഗിച്ചാണ് delimit ചെയ്തിരിക്കുന്നത്.i = i + 1എന്ന statement-ന് പകരംi += 1എന്നും എഴുതാം.
1.5. Another Application#
Exercises-ലേക്ക് കടക്കുന്നതിന് മുമ്പ് ഒരു application കൂടി നോക്കാം.
ഈ application-ൽ, സമയം കടന്നുപോകുന്നതിനനുസരിച്ച് ഒരു bank account-ന്റെ balance എങ്ങനെ മാറുന്നു എന്ന് നാം plot ചെയ്യുന്നു.
ഈ application-നു വേണ്ടി നമ്മൾ consider ചെയ്യുന്ന time period-ൽ withdraws ഒന്നുമില്ല. കൂടാതെ, നമ്മുടെ time period-ന്റെ last date \(T\) എന്ന് denote ചെയ്യുന്നു.
Initial balance \(b_0\) ആണ്, interest rate \(r\) ആണ്.
സമയം \(t\)-ൽ നിന്നും \(t+1\) ആകുമ്പോൾ, balance update ചെയ്യേണ്ട formula: \(b_{t+1} = (1 + r) b_t\)
താഴെയുള്ള code-ൽ, \(b_0, b_1, \ldots, b_T\) എന്ന sequence നാം generate ചെയ്ത് plot ചെയ്യുന്നു.
ഈ sequence store ചെയ്യാൻ ഒരു Python list ഉപയോഗിക്കുന്നതിന് പകരം, നമ്മൾ ഒരു NumPy array ഉപയോഗിക്കും.
r = 0.025 # interest rate
T = 50 # end date
b = np.empty(T+1) # an empty NumPy array, to store all b_t
b[0] = 10 # initial balance
for t in range(T):
b[t+1] = (1 + r) * b[t]
plt.plot(b, label='bank balance')
plt.legend()
plt.show()
b = np.empty(T+1) എന്ന statement, T+1 (floating point) numbers-നുള്ള storage space, memory-യിൽ allocate ചെയ്യുന്നു.
ഈ numbers for loop വഴി fill ചെയ്യപ്പെടുന്നു.
തുടക്കത്തിൽ തന്നെ memory allocate ചെയ്യുന്നത്, Python list-ഉം append-ഉം ഉപയോഗിക്കുന്നതിനേക്കാൾ efficient ആണ് — കാരണം, രണ്ടാമത്തെ രീതിയിൽ (list, append), ഓരോ തവണയും storage space നൽകണമെന്ന് operating system-നോട് ആവശ്യപ്പെടേണ്ടി വരും.
Plot-ൽ നമ്മൾ ഒരു legend add ചെയ്തത് ശ്രദ്ധിക്കുക — exercises-ൽ നിങ്ങളോട് ഇത് ഉപയോഗിക്കാൻ ആവശ്യപ്പെടും.
1.6. Exercises#
ഇനി നമ്മൾ exercises-ലേക്ക് കടക്കുന്നു. ഇവ complete ചെയ്തതിന് ശേഷം മാത്രം മുന്നോട്ട് പോകുക — കാരണം, ഇവിടെ പരിചയപ്പെടുത്തുന്ന concepts നമുക്ക് പിന്നീട് ആവശ്യമായി വരും.
Exercise 1.1
Your first task is to simulate and plot the correlated time series
The sequence of shocks \(\{\epsilon_t\}\) is assumed to be IID and standard normal.
In your solution, restrict your import statements to
import numpy as np
import matplotlib.pyplot as plt
Set \(T=200\) and \(\alpha = 0.9\).
Solution
Here's one solution.
α = 0.9
T = 200
x = np.empty(T+1)
x[0] = 0
rng = np.random.default_rng()
for t in range(T):
x[t+1] = α * x[t] + rng.standard_normal()
plt.plot(x)
plt.show()
Exercise 1.2
Starting with your solution to exercise 1, plot three simulated time series, one for each of the cases \(\alpha=0\), \(\alpha=0.8\) and \(\alpha=0.98\).
Use a for loop to step through the \(\alpha\) values.
If you can, add a legend, to help distinguish between the three time series.
Hint
If you call the
plot()function multiple times before callingshow(), all of the lines you produce will end up on the same figure.For the legend, noted that suppose
var = 42, the expressionf'foo{var}'evaluates to'foo42'.
Solution
α_values = [0.0, 0.8, 0.98]
T = 200
x = np.empty(T+1)
rng = np.random.default_rng()
for α in α_values:
x[0] = 0
for t in range(T):
x[t+1] = α * x[t] + rng.standard_normal()
plt.plot(x, label=f'$\\alpha = {α}$')
plt.legend()
plt.show()
Note
f'$\\alpha = {α}$' in the solution is an application of f-String, which allows you to use {} to contain an expression.
The contained expression will be evaluated, and the result will be placed into the string.
Exercise 1.3
Similar to the previous exercises, plot the time series
Use \(T=200\), \(\alpha = 0.9\) and \(\{\epsilon_t\}\) as before.
Search online for a function that can be used to compute the absolute value \(|x_t|\).
Solution
Here's one solution:
α = 0.9
T = 200
x = np.empty(T+1)
x[0] = 0
rng = np.random.default_rng()
for t in range(T):
x[t+1] = α * np.abs(x[t]) + rng.standard_normal()
plt.plot(x)
plt.show()
Exercise 1.4
One important aspect of essentially all programming languages is branching and conditions.
In Python, conditions are usually implemented with if--else syntax.
Here's an example, that prints -1 for each negative number in an array and 1 for each nonnegative number
numbers = [-9, 2.3, -11, 0]
for x in numbers:
if x < 0:
print(-1)
else:
print(1)
-1
1
-1
1
Now, write a new solution to Exercise 3 that does not use an existing function to compute the absolute value.
Replace this existing function with an if--else condition.
Solution
Here's one way:
α = 0.9
T = 200
x = np.empty(T+1)
x[0] = 0
rng = np.random.default_rng()
for t in range(T):
if x[t] < 0:
abs_x = - x[t]
else:
abs_x = x[t]
x[t+1] = α * abs_x + rng.standard_normal()
plt.plot(x)
plt.show()
Here's a shorter way to write the same thing:
α = 0.9
T = 200
x = np.empty(T+1)
x[0] = 0
rng = np.random.default_rng()
for t in range(T):
abs_x = - x[t] if x[t] < 0 else x[t]
x[t+1] = α * abs_x + rng.standard_normal()
plt.plot(x)
plt.show()
Exercise 1.5
Here's a harder exercise, that takes some thought and planning.
The task is to compute an approximation to \(\pi\) using Monte Carlo.
Use no imports besides
import numpy as np
Hint
Your hints are as follows:
If \(U\) is a bivariate uniform random variable on the unit square \((0, 1)^2\), then the probability that \(U\) lies in a subset \(B\) of \((0,1)^2\) is equal to the area of \(B\).
If \(U_1,\ldots,U_n\) are IID copies of \(U\), then, as \(n\) gets large, the fraction that falls in \(B\), converges to the probability of landing in \(B\).
For a circle, \(area = \pi * radius^2\).
Solution
Consider the circle of diameter 1 embedded in the unit square.
Let \(A\) be its area and let \(r=1/2\) be its radius.
If we know \(\pi\) then we can compute \(A\) via \(A = \pi r^2\).
But here the point is to compute \(\pi\), which we can do by \(\pi = A / r^2\).
Summary: If we can estimate the area of a circle with diameter 1, then dividing by \(r^2 = (1/2)^2 = 1/4\) gives an estimate of \(\pi\).
We estimate the area by sampling bivariate uniforms and looking at the fraction that falls into the circle.
n = 1000000 # sample size for Monte Carlo simulation
rng = np.random.default_rng()
count = 0
for i in range(n):
# drawing random positions on the square
u, v = rng.uniform(), rng.uniform()
# check whether the point falls within the boundary
# of the unit circle centred at (0.5,0.5)
d = np.sqrt((u - 0.5)**2 + (v - 0.5)**2)
# if it falls within the inscribed circle,
# add it to the count
if d < 0.5:
count += 1
area_estimate = count / n
print(area_estimate * 4) # dividing by radius**2
3.14132