module 2 exercise 1

parent ded3bd3c
#+TITLE: Your title
#+AUTHOR: Your name
#+DATE: Today's date
#+TITLE: On the computation of pi
#+AUTHOR: Miguel Felipe Silva Vasconcelos
#+DATE: 25-02-2021
#+LANGUAGE: en
# #+PROPERTY: header-args :eval never-export
......@@ -11,84 +11,74 @@
#+HTML_HEAD: <script type="text/javascript" src="http://www.pirilampo.org/styles/lib/js/jquery.stickytableheaders.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="http://www.pirilampo.org/styles/readtheorg/js/readtheorg.js"></script>
* Some explanations
This is an org-mode document with code examples in R. Once opened in
Emacs, this document can easily be exported to HTML, PDF, and Office
formats. For more information on org-mode, see
https://orgmode.org/guide/.
When you type the shortcut =C-c C-e h o=, this document will be
exported as HTML. All the code in it will be re-executed, and the
results will be retrieved and included into the exported document. If
you do not want to re-execute all code each time, you can delete the #
and the space before ~#+PROPERTY:~ in the header of this document.
Like we showed in the video, Python code is included as follows (and
is exxecuted by typing ~C-c C-c~):
* Asking the math library
My computer tells me that \pi is approximatively
#+begin_src python :results value :session *python* :exports both
from math import *
pi
#+end_src
#+begin_src python :results output :exports both
print("Hello world!")
#+RESULTS:
: 3.141592653589793
* Buffon's needle
Applying the method of [[https://en.wikipedia.org/wiki/Buffon%27s_needle_problem][Buffon's]] needle, we get the approximation
#+begin_src python :results value :session *python* :exports both
import numpy as np
np.random.seed(seed=42)
N = 10000
x = np.random.uniform(size=N, low=0, high=1)
theta = np.random.uniform(size=N, low=0, high=pi/2)
2/(sum((x+np.sin(theta))>1)/N)
#+end_src
#+RESULTS:
: Hello world!
: 3.128911138923655
And now the same but in an Python session. With a session, Python's
state, i.e. the values of all the variables, remains persistent from
one code block to the next. The code is still executed using ~C-c
C-c~.
* Using a surface fraction argument
A method that is easier to understand and does not make use of the sin
function is based on the fact that if $X \sim U(0,1)$ and $Y \sim U(0,1)$,
then $P[X^2+Y^2 \le 1]=\pi/4$ (see
[[https://en.wikipedia.org/wiki/Monte_Carlo_method]["Monte Carlo
method" on Wikipedia]]). The following code uses this approach:
#+begin_src python :results output :session :exports both
import numpy
x=numpy.linspace(-15,15)
print(x)
#+end_src
#+RESULTS:
#+begin_example
[-15. -14.3877551 -13.7755102 -13.16326531 -12.55102041
-11.93877551 -11.32653061 -10.71428571 -10.10204082 -9.48979592
-8.87755102 -8.26530612 -7.65306122 -7.04081633 -6.42857143
-5.81632653 -5.20408163 -4.59183673 -3.97959184 -3.36734694
-2.75510204 -2.14285714 -1.53061224 -0.91836735 -0.30612245
0.30612245 0.91836735 1.53061224 2.14285714 2.75510204
3.36734694 3.97959184 4.59183673 5.20408163 5.81632653
6.42857143 7.04081633 7.65306122 8.26530612 8.87755102
9.48979592 10.10204082 10.71428571 11.32653061 11.93877551
12.55102041 13.16326531 13.7755102 14.3877551 15. ]
#+end_example
Finally, an example for graphical output:
#+begin_src python :results output file :session :var matplot_lib_filename="./cosxsx.png" :exports results
#+begin_src python :results output :session *python* :var matplot_lib_filename="figure.png" :exports code
import matplotlib.pyplot as plt
plt.figure(figsize=(10,5))
plt.plot(x,numpy.cos(x)/x)
plt.tight_layout()
np.random.seed(seed=42)
N = 1000
x = np.random.uniform(size=N, low=0, high=1)
y = np.random.uniform(size=N, low=0, high=1)
accept = (x*x+y*y) <= 1
reject = np.logical_not(accept)
fig, ax = plt.subplots(1)
ax.scatter(x[accept], y[accept], c='b', alpha=0.2, edgecolor=None)
ax.scatter(x[reject], y[reject], c='r', alpha=0.2, edgecolor=None)
ax.set_aspect('equal')
plt.savefig(matplot_lib_filename)
print(matplot_lib_filename)
#+end_src
#+RESULTS:
[[file:./cosxsx.png]]
Note the parameter ~:exports results~, which indicates that the code
will not appear in the exported document. We recommend that in the
context of this MOOC, you always leave this parameter setting as
~:exports both~, because we want your analyses to be perfectly
transparent and reproducible.
Watch out: the figure generated by the code block is /not/ stored in
the org document. It's a plain file, here named ~cosxsx.png~. You have
to commit it explicitly if you want your analysis to be legible and
understandable on GitLab.
Finally, don't forget that we provide in the resource section of this
MOOC a configuration with a few keyboard shortcuts that allow you to
quickly create code blocks in Python by typing ~<p~, ~<P~ or ~<PP~
followed by ~Tab~.
Now it's your turn! You can delete all this information and replace it
by your computational document.
: figure.png
[[./figure.png]]
It is then straightforward to obtain a (not really good) approximation
to \pi by counting how many times, on average, $X^2+Y^2$ is smaller than 1:
#+begin_src python :results value :session *python* :exports both
4*np.mean(accept)
#+end_src
#+RESULTS:
: 3.112
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment