Commit d4b18a45 authored by Dorinel Bastide's avatar Dorinel Bastide

Proceeded to completing the exercise, first attempt

parent 8e98e7ad
#+TITLE: Some Title
#+TITLE: On the computation of pi
#+AUTHOR: Dorinel Bastide
#+DATE: Today's date
#+LANGUAGE: en
......@@ -11,84 +11,77 @@
#+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
* Table of Contents
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/.
1. Asking the math libraryo
2. * Buffon's needle
3. Using a surface fraction argument
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~):
* 1 Asking the math library
My computer tells me that $\pi$ is /approximatively/
#+begin_src python :results output :exports both
print("Hello world!")
import math
print(math.pi)
#+end_src
#+RESULTS:
: Hello world!
: 3.141592653589793
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
qC-c~.
* 2 * 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 output :session :exports both
import numpy
x=numpy.linspace(-15,15)
print(x)
#+begin_src python :results output :exports both
import math
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=math.pi/2)
print(2/(sum((x+np.sin(theta))>1)/N))
#+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
: 3.128911138923655
* 3. 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 \leq 1]=\pi/4$ (see [[https://en.wikipedia.org/wiki/Monte_Carlo_method]["Monte Carlo method"
on Wikipedia]]). The following ocde uses this approach:
#+begin_src python :results output file :session :var matplot_lib_filename="C:/Users/Utilisateur/mooc-rr/module2/exo1/PictureRes.png" :exports results
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
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)
plt.figure(figsize=(10,5))
plt.plot(x,numpy.cos(x)/x)
plt.tight_layout()
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.
[[file:Python 3.7.4 (tags/v3.7.4:e09359112e, Jul 8 2019, 19:29:22) [MSC v.1916 32 bit (Intel)] on win32
Type "help", "copyright", "credits" or "license" for more information.
C:/Users/Utilisateur/mooc-rr/module2/exo1/PictureRes.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 output :exports both
import numpy as np
4*np.mean(accept)
#+end_src
#+RESULTS:
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