From 6bc8347997fde6798cb44d09a3bd024635d8b52f Mon Sep 17 00:00:00 2001 From: boulanlo Date: Wed, 14 Oct 2020 23:43:42 +0200 Subject: [PATCH] Module 2 exercise 1 --- module2/exo1/toy_document_orgmode_R_en.org | 94 +++++++++------------- 1 file changed, 39 insertions(+), 55 deletions(-) diff --git a/module2/exo1/toy_document_orgmode_R_en.org b/module2/exo1/toy_document_orgmode_R_en.org index 2b73d64..a232e76 100644 --- a/module2/exo1/toy_document_orgmode_R_en.org +++ b/module2/exo1/toy_document_orgmode_R_en.org @@ -1,6 +1,6 @@ -#+TITLE: Your title -#+AUTHOR: Your name -#+DATE: Today's date +#+TITLE: On the computation of pi +#+AUTHOR: Louis Boulanger +#+DATE: October 13th 2020 #+LANGUAGE: en # #+PROPERTY: header-args :eval never-export @@ -11,71 +11,55 @@ #+HTML_HEAD: #+HTML_HEAD: -* 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, R code is included as follows (and is -exxecuted by typing ~C-c C-c~): +* Asking the maths library +My computer tells me that $\pi$ is /approximately/ #+begin_src R :results output :exports both -print("Hello world!") +pi #+end_src #+RESULTS: -: [1] "Hello world!" +: [1] 3.141593 -And now the same but in an R session. This is the most frequent -situation, because R is really an interactive language. With a -session, R'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~. +* 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 R :results output :session *R* :exports both -summary(cars) +#+begin_src R :results output :exports both +set.seed(42) +N = 100000 +x = runif(N) +theta = pi/2*runif(N) +2/(mean(x+sin(theta)>1)) #+end_src #+RESULTS: -: speed dist -: Min. : 4.0 Min. : 2.00 -: 1st Qu.:12.0 1st Qu.: 26.00 -: Median :15.0 Median : 36.00 -: Mean :15.4 Mean : 42.98 -: 3rd Qu.:19.0 3rd Qu.: 56.00 -: Max. :25.0 Max. :120.00 - -Finally, an example for graphical output: -#+begin_src R :results output graphics :file "./cars.png" :exports results :width 600 :height 400 :session *R* -plot(cars) +: [1] 3.14327 + +* Using a surface fraction argument +A method that is easier to understand and does not make use of the sin +function is base on the fact that if $X \sim U(0,1)$ and $Y \sim +U(0,1)$, then $P\left[X^2 + Y^2 \le 1 \right] = \pi/4$ (see [[https://en.wikipedia.org/wiki/Monte_Carlo_method]["Monte +Carlo method" on Wikipedia]]). The following code uses this approach: + +#+begin_src R :results output graphics file :file "./monte_carlo.png" :exports both :session *R* +set.seed(42) +N = 1000 +df = data.frame(X = runif(N), Y = runif(N)) +df$Accept = (df$X**2 + df$Y**2 <=1) +library(ggplot2) +ggplot(df, aes(x=X,y=Y,color=Accept)) + geom_point(alpha=.2) + coord_fixed() + theme_bw() #+end_src #+RESULTS: -[[file:./cars.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. +[[file:./monte_carlo.png]] -Watch out: the figure generated by the code block is /not/ stored in -the org document. It's a plain file, here named ~cars.png~. You have -to commit it explicitly if you want your analysis to be legible and -understandable on GitLab. +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: -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 R by typing ~