From bdc12117a791f51ffb45be2b2280b76b0e9f2828 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Gi=C3=B9=20Marcer?= <giuliamarcer@yahoo.it>
Date: Tue, 19 May 2020 16:01:52 +0200
Subject: [PATCH] ex-5: revised and typo-fixed

Some images and an article have also been added to the folder.
The todo has been updated too as well as the readme.
---
 README.md                                 |  45 ++-
 ex-5/main.c                               |  29 +-
 ex-5/plot.py                              |   4 +-
 notes/docs/bibliography.bib               |  32 +++
 notes/images/{MC_some.pdf => MC_MC.pdf}   | Bin 17005 -> 14782 bytes
 notes/images/{MC_all.pdf => MC_MC_MI.pdf} | Bin 17593 -> 16296 bytes
 notes/images/{MC_few.pdf => MC_MI_VE.pdf} | Bin
 notes/sections/4.md                       |   3 +-
 notes/sections/5.md                       | 319 +++++++++++++---------
 notes/todo                                |   2 -
 10 files changed, 276 insertions(+), 158 deletions(-)
 mode change 100644 => 100755 ex-5/plot.py
 rename notes/images/{MC_some.pdf => MC_MC.pdf} (60%)
 rename notes/images/{MC_all.pdf => MC_MC_MI.pdf} (69%)
 rename notes/images/{MC_few.pdf => MC_MI_VE.pdf} (100%)

diff --git a/README.md b/README.md
index 070dcaa..3c33f27 100644
--- a/README.md
+++ b/README.md
@@ -42,15 +42,15 @@ See this [guide](https://nixos.org/nix/manual/#chap-quick-start) if you have
 never used Nix before. Running `nix-shell` in the top-level will drop you into
 the development shell.
 
-Once ready, invoke `make` with the program you wish to build. For example
+Once ready, invoke `make` with the program you wish to build. For example:
 
     $ make ex-1/bin/main
 
-or, to build every program of an exercise
+or, to build every program of an exercise:
 
     $ make ex-1
 
-To clean up the build results run
+To clean up the build results run:
 
     $ make clean
 
@@ -82,9 +82,9 @@ The output can be redirected to `ex-1/pdf-plot.py` to generate a plot.
 Every program in `ex-2` computes the best available approximation (with a given
 method) to the Euler-Mascheroni γ constant and prints[1]:
 
-1. the leading decimal digits of the approximate value found
+1. the leading decimal digits of the approximate value found;
 
-2. the exact decimal digits of γ
+2. the exact decimal digits of γ;
 
 3. the absolute difference between the 1. and 2.
 
@@ -104,7 +104,7 @@ cases the best fit and the parameter covariance matrix are printed.
 The program then performs a t-test to assert the compatibility of the data with
 two hypothesis and print the results in a table.
 
-To plot a 2D histogram of the generated sample do
+To plot a 2D histogram of the generated sample do:
 
     $ ex-3/bin/main -i | ex-3/plot.py
 
@@ -123,7 +123,7 @@ horizontal component. It is possible to set the maximum momentum with the
 option `-p`. A χ² fit and a t-test compatibility are performed with respect
 to the expected distribution and results are printed.
 
-To plot a histogram of the generated sample do
+To plot a histogram of the generated sample do:
 
     $ ex-4/bin/main -o | ex-4/plot.py
 
@@ -133,11 +133,30 @@ and `-b`.
 
 ### Exercise 5
 
-`ex-5/main` compute the integral of exp(1) between 0 and 1 with the methods
-plain MC, MISER and VEGAS. Being reiterative routines, it takes the number of
-iterations as its only argument.  
-It prints out the obatined value and its estimated error for each method.
+`ex-5/main` compute estimations of the integral of exp(1) between 0 and 1
+with the methods plain MC, MISER and VEGAS with different number of points.
+It takes no arguments in input.  
+It prints out eight columns:
 
+1. the number of sorted points;
+
+2. the integral estimation with plain MC;
+
+3. its uncertainty;
+
+4. the integral estimation with MISER;
+
+5. its uncertainty;
+
+6. the integral estimation with VEGAS;
+
+7. its uncertainty;
+
+8. its χ².
+
+To plot the results do:
+
+    $ ex-5/main | ex-5/plot.py
 
 ### Exercise 6
 
@@ -167,8 +186,8 @@ classified data in this order: signal then noise.
 
 To plot the result of the linear classification pipe the output to
 `ex-7/plot.py`. The program generates two figures:
-  - a scatter plot showing the Fisher projection line and the cut line
-  - two histograms of the projected data and the cut line
+  - a scatter plot showing the Fisher projection line and the cut line;
+  - two histograms of the projected data and the cut line.
 
 `ex-7/bin/test` takes a model trained in `ex-7/bin/main` and test it against
 newly generated datasets (`-i` to set the number of test iterations). The
diff --git a/ex-5/main.c b/ex-5/main.c
index bfb6014..39931b5 100644
--- a/ex-5/main.c
+++ b/ex-5/main.c
@@ -15,14 +15,18 @@ double function (double * x, size_t dim, void * params)
 
 // Function which prints out the results when called.
 //
-void results (size_t calls, double result, double error)
+void results (size_t calls, double result, double error, double chi)
   {
   if (calls != 0)
     {
-    printf ("%ld\t", calls);
+    printf ("%.0e\t", (double)calls);
     }
   printf ("%.10f\t", result);
   printf ("%.10f\t", error);
+  if (chi != 0)
+    {
+    printf ("%.3f\t", chi);
+    }
   } 
 
 
@@ -33,7 +37,7 @@ int main(int argc, char** argv)
   size_t dims = 1;          // Integral dimension
   double lower[1] = {0};    // Integration lower limit
   double upper[1] = {1};    // Integration upper limit
-  size_t calls = 5000;        // Initial number of function calls
+  size_t calls = 50;        // Initial number of function calls
   double integral, error;
 
   // Initialize an RNG.
@@ -65,26 +69,37 @@ int main(int argc, char** argv)
     gsl_monte_plain_integrate (&expo, lower, upper, dims, calls,
                                r, sMC, &integral, &error);
     gsl_monte_plain_free(sMC);
-    results(calls, integral, error);
+    results(calls, integral, error, 0);
 
     // MISER
     gsl_monte_miser_state *sMI = gsl_monte_miser_alloc (dims);
     gsl_monte_miser_integrate (&expo, lower, upper, dims, calls,
                                r, sMI, &integral, &error);
     gsl_monte_miser_free(sMI);
-    results(0, integral, error);
+    results(0, integral, error, 0);
 
     // VEGAS
     gsl_monte_vegas_state *sVE = gsl_monte_vegas_alloc (dims);
     gsl_monte_vegas_integrate (&expo, lower, upper, dims, calls,
                                r, sVE, &integral, &error);
+      
+    // In order to see the parameters of the VEGAS integration, decomment this
+    // paragraph:
+    //
+    // gsl_monte_vegas_params params; 
+    // gsl_monte_vegas_params_get(sVE, &params);
+    // params.verbose = 1;
+    // params.ostream = stderr;
+    // gsl_monte_vegas_params_set(sVE, &params);
     do
       {
       gsl_monte_vegas_integrate (&expo, lower, upper, dims,
-                                 calls/5, r, sVE, &integral, &error);
+                                 calls, r, sVE, &integral, &error);
       } while (fabs (gsl_monte_vegas_chisq (sVE) - 1.0) > 0.5);
     gsl_monte_vegas_free(sVE);
-    results(0, integral, error);
+    double chi = sVE->chisq;
+    results(0, integral, error, chi);
+
 
     // Update function calls
     printf ("\n");
diff --git a/ex-5/plot.py b/ex-5/plot.py
old mode 100644
new mode 100755
index 4c74c51..d297e70
--- a/ex-5/plot.py
+++ b/ex-5/plot.py
@@ -1,8 +1,10 @@
+#!/usr/bin/env python
+
 import matplotlib.pyplot as plt
 import numpy             as np
 import sys
 
-calls, I_MC, σ_MC, I_MI, σ_MI, I_VE, σ_VE = np.loadtxt(sys.stdin, unpack=True)
+calls, I_MC, σ_MC, I_MI, σ_MI, I_VE, σ_VE, chi = np.loadtxt(sys.stdin, unpack=True)
 exact = 1.7182818285
 
 plt.rcParams['font.size'] = 17
diff --git a/notes/docs/bibliography.bib b/notes/docs/bibliography.bib
index f73e77b..af6bf5e 100644
--- a/notes/docs/bibliography.bib
+++ b/notes/docs/bibliography.bib
@@ -106,3 +106,35 @@
   year={2005},
   journal={Matrix}
 }
+
+@article{sayah19,
+  title={Adaptive stratified Monte Carlo algorithm for numerical computation
+         of integrals},
+  author={Sayah, Toni},
+  journal={Mathematics and Computers in Simulation},
+  volume={157},
+  pages={143--158},
+  year={2019},
+  publisher={Elsevier}
+}
+
+
+@article{ridder17,
+  title={Variance reduction},
+  author={Ridder, AAN and Botev, ZI},
+  journal={Wiley StatsRef: Statistics Reference Online},
+  pages={1--6},
+  year={2017},
+  publisher={Wiley}
+}
+
+@article{lepage78,
+  title={A new algorithm for adaptive multidimensional integration},
+  author={Lepage, G Peter},
+  journal={Journal of Computational Physics},
+  volume={27},
+  number={2},
+  pages={192--203},
+  year={1978},
+  publisher={Elsevier}
+}
diff --git a/notes/images/MC_some.pdf b/notes/images/MC_MC.pdf
similarity index 60%
rename from notes/images/MC_some.pdf
rename to notes/images/MC_MC.pdf
index f995a9a9fc28f4e38a1127ebe580fa1736decc6c..2f0e06130febf694e6a6c5f1d9022dc81a2bf24e 100644
GIT binary patch
delta 4619
zcmZWrc|278_x}zNV_##4Yi-6h_s%}DhAi2)R3!UWrWvA?Zc?OFxRw%S%Ogt)Pg+q=
zo-C;em6TRnkI*8?MBh8|`~DvJ{BiHO=bm%k@ALkA?s?tuLzt$ROZFHc3Pd&{A<mdM
zxRc{cl^%Fi%Ox3TR(&1K&4ombDoTzXsPb_&aMwLmw!K2zVrJioi46mY^2tpZS3le@
zU~Zgqty0HvSZg?qMpqkW=0^8(%-{X=SG}hH<E1lZa^jENwwTt=8{;A6gF&r4aO*tH
zso8ne!;gI{U2AnUr%(GC&8X<*Ubry3jt$#R#A&ndR2E=rTnNE_bfdAC>xiTORCOzM
z&#sYsJgTLWSzQ*Hv2V)?mpqX7Zb}kP-iy1Wd%3zn-s1SPMnx+&{aeq_u5ZaBS=Suc
zPn#WA`M<%ScW0MTv3LpfbsBod%h=;51V-%p%p%-RyH|tRiJ}Kr1u0ZeYURSFkHKT_
zhbSuaR$&(pf}vgzW#$0nw)^~0fA<JocQp<X{BDtzmJQ+dwVd+YQ!J$DdQ!>b{ZF!v
zdcOO|;Bv(vkMi)37!QuVSj`iGf5K?zRrBxk?mI4Dhxb4$5mY}rh>7+-R14e}o_=>&
z%cRprx%BqRbeRr)L)nURd;G1>+=02`+Uo17yi!jI_vBjKAhZ=5P!anCEM-s}?sct;
z%MzBjc7(;0@pAGm(unz#sk6CGKw?t@PV0d^q7dHRc7LZ|dWwmK(+N)#4c+9eo03xs
zNGG)Hxu0q+tsXp1#RZbS`y{3KzW#}m8*_S8q)|FmR<NbJ)jc`I;2F7bTRV&*^^>=E
zRQFsIHdMCJqEvlup-|r(%JyZ?D{M?oNsrHQPct5-gtv**`tqK6yhu{ki+>f?8{g%H
z7SnVHCG;Ch;5c1ob``DR#THi4reE{$vzH>dS$hpjeJeGjT#*=urkrgw>wuKa+(T~m
z!xj<TjMIrWW38J-C2)0c?8w-*=Q37SMXncZrb5HH87ePDHD82RS!`HudWXF|uuZUc
zP&OgSX%voSqY}JZt;vl}5tcZuRfkAN2~tsB)jnb!=XCYC8ByO%9FpCvl2XzoOR_T@
zzd)eVkaVI#nRKEZHN3q7<Ef%xrG~QhEg3#suH})O;#+bL7Ome>Vri9oq}26H+~GUo
ze=5G>i6MGuzcsw3%kBBN-eAk99*1<xwWzd?O2k8mjy^Zb4dZ&noi>QUkv<dpGuO{%
zqr7F(Jo}>4ge8;}0IL#qJ!(8s0Bx*GQv|5Weuq*y&w`yLa3oFL9Muvt9UNuQ?pCSB
z25-1XrxbhF$AQ8x=B`(IYr||oIG3Uo{?Mp)U{oLzzt2dyOQX`ia>x70voCI5zro}3
zj8mJUBV3=h=O*XrbQK@7ETBlN&P#e2Whd?STBSkvqZEf{kh(9E<0&cR<7%<PD&x)r
zyU4fDA5(GSF@t((RtzITX`@Z8Ku>B}5?q*R&74vDVHc-amv(ebdS%Vi$MPkz)JriX
zdTJI2M^tfjLwc2}jabvdlM>HvNxEUzYGK*KQMFnanFg%!Wj8l&;nAX>-Uo%_oTYJU
z@<6qhW(QZ(!TE#dx=K--{G+!pg))a)zMO}Iwb~)pLFW%1cY+Kb<9zIKeN49_F>tK%
zsLkdJTc+Ck{)qQb?>{^Jsiu+PddIC)&8FJ+N`urEw#u`nv$lq1a|gaFT<+p9@>-$1
zIS4uba6MziDY4JSFB5$>7C+H*mUFNEzU}Pb%WKIzVj4cE&2-?+5y`ko{~PK9lagC^
zsTUhv*DmO^_$V)FJ*FiRrZJOql~k}5ZcBuRjLtUC^SZjmCToobx>lNA&l%q0|Kxs*
zh_{q!|JyI>QzE2?;?(%6`{7x$J$R^7u=1R<-MciQd-tZFiIU?KWA`+e)O}5rF~ylZ
z?Dp(ACGc%$&i7EM;Kw7{o_$>GpWAuo=0y}O3%9Ddh3ws06S$*p$VA9;I-<?5Qwy#X
z=&I7#7`)|_@%<x>T@!f|k5rTf=sewyrnJ7s?A&-6^@crks|#P-bT5YK-Y3+5OMQ3I
zGwG*S?Y+ajyyrdE@zQ(q8HL(S&3?7s!%>e-6?sx80%~tFwx?yiuD*VZmArTUUQy}0
z`{C2M>b8mf_TiWGX^m$`{?bb9usv*T2%C29cQEZLJM^^JH1KO*e1pR9{&VcD-(H51
z*E?p9*^bI`qph8(!FPWaH5bI+O!b)=?B|8cU(nBTiSFjvGp;5#o5!t`dZ}eJxiSi#
z>WvI|0$;nnCWZ~)$T@H%=~|c{VQ}wu?dVPc)7t&<-~Q>EoOn)b6xjL6;44e;3d;zd
z9(sLtXz#H!?gtH-B9jL8SW@u7nTiI#73<Td<lRJqS3*VKuW~z7UlSr6FX4UYJuM>k
zW#EbFo3AZ|gz%!d-tKtMBkx}nKWK@j)vSzNd8F9(j=7LWb8)a?uq`cQHjv!#eqW>D
zjeyAUu6e=cfhm^1Q6tNL*zsxO&IfUS!DA&o_u8xJtM5jwa|%nfJeV<czWV85VPlcA
zCwX+!auyKb7uJxZ^gQBPdE)1{=Qh&e>&^ZHb6*pV<&$%!!RM>J`}=nV;R|_Qy>lr*
zsbO2mhf;62dGnVxantecN9iBBRZ`~7UP#FI(1!)f3MF$6x>a`h)MS@fwU(QOPQ%<5
zjZUWJPFrPYg@5iI^PL|!!v*mccdpNV_OazgYSiWo_Zt5vs%5HTMPRfoZgkInTt3G7
z%)#+p9Rl`egsNP6GENH>izt_En)|%!QIxk8=3H%UgU<m=pqMbnVS<4|*@a>yzR(lp
z`F^3udlUH=Rx2r{7aI3H{916BSzIV9$^QGJISjE42u%g^La`n1-yokvR|`-n$VmwW
z#6U?EvBY``Knz4fX-$0-b^``c!V>C9;-MG}6*(<I!O)Ng62~wUWDK%HN+gYB$p8qr
z17>E}kmz7mU?dizAW<X+f|JS<<-g*)wjkzY>w1_xgsFcgGltQpE&5GN0HQqr(HANZ
zSpdZNWyJsxbHOq*Ltt=$ReguND4zk!3gIcR5CM4tVz0CdDMs%jO07t@0-;`SwJQdR
zP}Hk0QxnD@f2_jScdIL65I_lE|3-5a25C|vgfT%r)tFdLNWc#4?gFRdma;8?h~VT+
zO&3>d?$6kS?w|fuy^^<DoCVSpZ4ORRBR7B(S4j#{u7p$Bc`aU`uUP*?Nyj|N-mx&H
zXK$lma@{9(a<@Hinw{{Us<{N=0jX|PO&=scRSv09#v^O>q>)?N@`#PLInswCAT`>K
zNQ^20(bTz(U{({5oHZ;YL5+YsTqTW|;w_L1>V*1o{B{g70ua#TVC0%6zW$z`9|oDx
zA|Or%AaY2XfRr1QBHy(MFqzC)(hk*~d2zr(d@K>`2ONm&Di{&Euz%uJO13}Ns=@nc
z{t3mLi0>g9UWSNK^SJ$vG?mZa+OcVJSsQiw&x>ceem6Pn*&Nn+v*}5J46q>wiy&F|
zbj0Sx^IjUMjs!)SpVW#sP1xp@en=M6DKA$VwFx%9@mK!wA137(^S3gbi$0C2K?-&f
z*V+9FlU85$wBdlzQ@Ka6e%@({`z=2e+14a}^YQcn9*wKo?@T}Lnt32DsozeVS+97b
zf;hV=`;=(x_%Y3ik``I}`vrl1#<Nmpu?7{*w0E7gCmd$&^rIa6<4f%$?rS0!^K4jj
zR-$R(Rf{i19q`sJz|!LhK4KrYvM>mruN-nXUNv8wa!)n&zMM%aSvGqr5x?JmgYI>F
zrn#g<G3>8anqDBx!o)l}cjc(ys!z`XYE`!#?rzN=OtjFGOV~F1^Xbf<`5%I+R}Oe0
zk95|;B+62m4c(Yof|eqkBMJ%YVtOVocV>S3`5xPWNlvi>j1AeYf4fNM>c9Nl=GE|K
z^0d<J_}j&=I(DfjGNl#lx|KYgD*SG$-C7x^8Y=LmjB1dqJSJ!!grz62cRT+5Y<!n8
zb;aJNQ7)mbB>h$F_PatJyWqaI(|@^?r;CmIK9p#?Ci5gjsq9S7-|2_6oHX3+ivQ8Q
zK(@?SQCj6zdR(=2|NBm@!T>L=NO7L;hUZUX<R0KCYZa?ZnCERb7Ie0^is`SMyV=!Z
z<B6*>k(S-r{M>4;;KLyOW%srE)~T7TF?BaV1)opDf2w3S{W(J?zMV0BCo=Z^hhW9|
zEsCfZa>$M~c!XoviLAgAkWiyABoV#aDM=wh#4tow4?WSuFBlS%1{z?<$Q97qg%}&p
z0*JBUoRAQJ7|PnkL9Y>Us2vlr8@+_Yg+ya%bbx3D(3k+x2B4$16TkogqAS2e?Erv@
z`hx%ngaBd)fP_P87%Y+s1z{j`21Sd2|2HWx;tEI{>82>(@$ad_6+|LQ;;<D!qWl^!
zHUNqG%W?up{8{AA$099`p#XriIEI1%(yuXu!jczYu)nC2p<hW0-H=JY9wGqpuTe<c
z4v-g-q5<+E5(|z-i!6X)i916$5v*W<yoeYBpr7#qiv=|R@**l<Aj+@yh%q|@f<gew
zB5*7~SwO`u>PlHe=1Wak#NGu^7O~?13Jp+P06+^af+mD;&{zg?j-f8T*ui=d08x<+
z#z6forZfh5W?~{hVIr$c^D!V&XG(<06ee1mzkp}}yI?>hCLP6I>LY^?+URdSh>G_3
zyN`}``kRl0R$uBvpJ-?&AhOMj4wHU=3z13u9avHYBrj(LlIaUo(7^w(2gwxvmMqPs
zg8VbF)JLN7l}8u+uVfO9wk(hg^0#tHG6d0<%K=fA`#|P0f)I+aL?W19nhMf?2LY*M
z#_#(8QmOoNxP%Nyr7@SsG5H62X&j9VF2kVFAb5EolmC${O{SA5%Y77RSuUMQUGAgN
zmI2e5;Nlj;|G#$(kh-ip21Mtd6m<P7j*8Bb-+lCD=Llp%(6TrtgwCFS=koPuQov;r
zG8xeCWYBrV|M>aqe`Fx^6fP|Yf)Ge0{az0OX^h_wHChe*Nn5>6g~}I~KeOCftXKdQ
ikZ*fLB!mFe#h;r1PV8<@NGKK}(?JqeN5{rxJ@&u;<mQY3

delta 6112
zcmZWsc|26__isqazAxFXeJ0DyoqdumWZ$<zl4WcOA+pr8SWCh!Nt5zZk@!TGvXqp4
zUz4>|mK2dn5&bTCeLmlK{+K!U+~=J4`#I-4&vWk6;F)Mdp&^0=M<J_{Au~iM-7r9+
zf%j?0gQ@Hwi{kqc1=%7>(afyzax#?d6xNc~S0_rOQAJ-C7W(|0)uRqx2Bb@>7EMiO
zoJdV8`%)^8*Q_K~yNtMuyq^2boir@CI%<Y^J*$xL-uOjB&v(#hp<pa>O6s?Usiv9m
zgSkQB<U^*qbvbXXz81W_QhBd#0!N=uVuZ?hemP^sRPvOTd<Bh;3dj-;`PSym?z80K
z#^<7EM3e11ThpL77qY9k+U(|BKp>j=M34cyC}B_2v9G9-<dL?Kknn<2a~B=kW|!hV
zyybqMTk%mUNBFB+7$^7i3uT)-_~ZQEhoWZ^dhL)w=dRUFTFhq!OQUf&_gPodpJtAl
z75iA9@51FUvcbBum(o*nRXli%yi8G)E5C1m<K8BQbyww33hs$4Qs7<h=O_oitoOge
z4(1W=XyMK`J(4LJd`wL)8@96!`X0aE&m)18%wzu@W$jv2yF)zQX`aBgN_AuAs+EzX
z*`7aM&7uOnARyZP!*jE^xN24CKa<Z3eB@+zeA^f=$3ES~;}i&(HM$xy3fhIl<BfI5
zO2STyu;lbIV4%MCa2h`5O&7<GvK?Xu2`9}iL&vYP>~w2;oA|j*%`p&2o93^udYD;w
zOTHak-j)(uqQmSe?P7>@JgGwFA$hnSHX^pdzEY^5#bxbX%#Ln0x9A<FL4~!jk`k%?
z{Nb9rw^U4;n_(ZHSSO$TjKE%*e3CMmN0#0CDAuMFew(3QPX+E5!?4=BgE;+kj%phx
zo*E}q#<_acoJqf*)bnZwyxJ9}le1q2zlEMJ(e+gbarmAJDPR0j*B00F;R4YhZbx~@
z((wpcwP0Nxy2+sTtyJn%RU*?4>HX1%=lb&@kx|CDw^dqkP^XAGl`V`S5wV|(R;cyB
zFMn4Nd!82yNRpB@Mmnd~Wb)MI&&uD}&v_$?zpwp5YP_ms2BG<|d4g|Wa0<n!yCBuH
zH^(q;+F0U&XiKyLsdP$y+)dwiI8Hi1_S<r}fKW#kJ^J>!Y{^N5#qysYO*m{ZCA(g5
zo%9v0cCk|x(d4|bOIt*;75WaT<tJScRpB_|mYZ3KdfXv3;xdd*imDe9iHa39WD-wA
zYx7y^CG_0BQIREr(H^N(3Y#Bkczvqk=a1Rx?vbc#cNaqi-tpJ`)BG+vDS_8r+sUq>
z`jHMx$4PoC)6T>fJ4DBddU--G{}k*|6!o+<G70~5^HpDz?^%Y+AyB^K-qN7&ZEP7+
z`(TrAENww4V1yJg1A66Zh&u7d*1zGbw=y%~eX*2D3c1eSJ;6kIV`_CI=(X4h8{85}
z>Ip*pqmI8<Rn{e-TC%F>4cY8h+^80p&k~Ss@ywIn98yn7j@jnEt6J@?iUo4=#<|4X
z!?S8To#OnBylyaG@Kv8=#J@e&NO?Y=A1v_OZnALs-R!LjVbSS5_2Xji&b0OkVeLLe
z>Xm)N-yaN}T1J13IMhUcdNRFJBI2ApLsITnRHdY8Ny2OqcSK`gbZ}ss^x<oDT}*Ag
z&Q-f+E9ov_5dul(iB_!M9shh_v%s$~i_ZysYe_6uMU6FDO`5HcDHEIx;^otwnBIWu
znkH^1Hc*lfUn(JCE%g0s{A?rOXUVyHTrpuR;xdoWY4LZ#@j9p3`GX0yuieI}KOe;%
zz+5z!3O4LwaxHmAefeWfR@JgiDUJt>x~Xq^`3T+YgU|E;Di80v8uqph<=8nky3#a|
zLC~*@w8!SYd`D{^DQ^*<VVJ*WI-a-E7xy^E;rBRKe5}~5p`+4o-$1i~ZJcfXVGH3M
z6|`R;zK-R%3}1g8@UgS*er~_fkl?F!u2~1wXA0da56HbP!%2^*6Axz{#X>VKRn<~2
zOmkuBYz$#6!;{tPrgu?d2%@j87i5y2@L>Er&QobsAD`RKs;zpslyCeos#u3~y;&`Y
z$*Qwk=G0yz#>B~T_E#yk+~auJ>t#6-uS~HFQVX|aWfyx@%M4RMCUSRGcwPwIe(0_g
zL$ycRS~k=P$*{SRLG9nO=jUk&zR8cLoXLLlPHLEEPF>!}o&tSF|L>*11@6V6xE@VS
zf7C*4ceuFn&Wl}pZPlKbqkKjIJKYYq^U8x77nUdc7+zd2?iuTU9=y#yNIft#jil+G
zb~MFB|HG%}Sx7iv^|i0||5D_8MI!ignPU-%e8$;g0l$BE(7h>T2Gm0g-YyA8n=02B
zI-f$4w8On;Mmi5}x{M0;YG2*`(YT)PZeMEKqmq`oy(&L)jgIfJ(bbZ<^r^6;o|sdd
z!hyGG<hk$BtnWm;Jya<@q8Mu%=6u!n$*>?E0j7OQQL7D8d4dUhq#q<1*ZO}g<K<l8
zvn7SiU4lC_3D@ylRjx+{j>`A<n7`>*9h9V(axMN~;%8XiZ`&@du0$A-_gih+H)`GL
zREBuw-A*$%U}V`NB3BKcA9BVFA6-h2WYI%9WwE-vmkT%0tnr}+ZdWBPZ}Sh$|1Ko^
zNWZQ7=zJ;1BNGwb#oM5xYB|l%0deFr+wq$sW+TDQxyOU-AH=<Q#hdg!RyC}>I+)T=
z&qLD`x^G7=tSt9`o+Z3mP{cx~S>)5dU4=pWDrB6*N-HHK2EH?DJCSNR>X+39b+`<C
znJyY8S$*}rUk^&NpLdj-9pQf6b-{XC%c=bH2Wss*s>7E~&LeXscO1{yPTo2EsHlnB
z)S#OhXU7^?diWdn@oVbQvA*o3<7M(Oo`OR?IUgUrPY2yt5^pGVJ9qsIq0cZCxQDfw
zY_DKZw~9NVQOqJF{$RpI$%nZ5^IPBjMJM9Z8>^gSNiR>b(fmDvLxO0oN7tkPH(xdo
z1TeswfW=Tq|0*fa83{#5h|LCpg?MF<#Soh@0!mlrgnmg$6ti;B5C{_F&xuEnp={0z
z2qJ{$3gsifj;sI;jengGp$ug*^qK1_f>0dCExw&V-hkS-4?r6NAZ6VE9ROf4f23dl
zyT;~!onVP;OlV~Cn&S=TCm4jpO3r5(R^7<Ek9mKPPcJ4|rK}9`IH^94?B<`OJZse8
z)}N|X@z5E&d%}>CrD2J84A@&{5L{j5ffE*7!402`nVC54fj4|F-8CLG?tqHWxH>3y
z3Q`hMEB?*D3jrw$BO#1n0CYzLMaSd+f;P1@^smKmR4xntdbTYszb@STeCFV97vsGm
zCfu<i4sZFw6crxMJ&dUk+AHot{9^Qq<CkQ~L~TdB;O{P?OhJ^z(K9Bsx+$G{Bi_TZ
zBXfInms8jd-5Ygp(;qIt2Yh%{w>rABCqq9kpKU(bJ=jB*?&Hj_*Ra?b(ZSh;I=6jk
zI9i<N4kbb`b#;)tzqISq%fZ&nTrSCLe9ek&zTRdz>C(knl8^3Dn~=k~*gdn>qCV=4
zL*$liW=Z_b-s%HI8S!7NZywx#$fR~x<)=^7cvEFR@$OSLOSNAX3P|q)?v|$QuV+c<
z6u3ItapKZv%@Jm*QwdZfss`PZ;!%Bwx5Zc#mTcpTO`BkwIFQDMian^M5r3tHg}((+
z%yAk|_2Dwn<Fenq?T#+Td)SO0;uAwc4T7PNr8o-m6sm(3#8GrCmh@Ll^^%Q~n5B7i
z78%d=Z><JqCOw=UrmLy$?eZ&>NA}~&nfW8<8%-|zS9v^EYZbgXCW;CEVU+GyYcgea
zQ1<da0%JBeS(92ut3KRHh}3P*&SY`uYv=Suc%@F9ojI#InJ7qVZz?!2-u5-SpXaGh
z7+tDAJnp-}yU`amPunSd?Um<Qeb^oPN}eYB2{dsX(cEJ?)ZW5V^u}+u-IVw!Yo&dd
zqOr)S>Mt01eogL?CCeyhQuapyFGYKWBOT#Wh96wMHD`0z9YakC7tOTFExP1br{1%U
z^!n_xlJmG>nLF9vj`G*OgW5%RK&K^;(0e&Px^b$LHtSHe^%T=TzN>BFjEuhbt3#}x
z@^76H;!74HH`W?+SY8aOlL|jOzxPb_$8n!u8ly42XLDp$o|XEDeQwAgSQg^lP^}u@
zSCG;T*6&{POqCn(dwO7l0MYp6QVsgc>n4<6J8t^~FCh&hcvN1E><*1`NY?I#cBg%M
zA4fl~Qdq%inl5%kMwC*DC3G}oj$t!JWVQR00mGY*%>@$;aFS|XD&fgvC-Q{Rr=DEw
z$6mU3)s_WYwrk?B+O=IdNEx2$TXN+ImEP0HB+vScp_tfInZ*%x2IFf*HA9!>(5^A%
zmwZX05zSadHBO^jzT?}|U7vi9IhH{*dTedcGvi$J;I)^6>+nQju@#oH{a7ISdMmG7
z<Ot}>FBzA<ILn(@TlLe~kN04CYq>(tr4{6$fs{c1<U!MG1;{a}y!Q!${Xc&H8tv%N
zOeUQyS}(e95=d(j=C2yV6=N*n+TG@c3&ybo%NoZB=%~uF<8EoGUS**IaV?#lJk3}l
zXobP!Sst3SUCeUBT-0(fPS&PatWXa6Y%{mvL$CXOEr~3)5w<`;p;D;gHj$kOXif@Q
z{7X~_0mVoo={O4ZFI;siV}BNqb*KaNSnJN5X;9|9(b*BY)Ad%HJT^3Q8F=|B$31z?
z>C||f$*FHa&TlQhJn{4OnvkAVCEoo(4RWc7eHGM|;Z|fW-h^%}v7u(!jb!W0#~}H+
zx2e22!!q_%{P(<YZcTq-Owfx+4<`Cyqh5cldYvv&0|#wh*X0z!soNG0{Kh^b0sG$k
zn5+=Dc+RwKsHgWi6Vh`BZpB4kNzIe!o=?0z=TzgI;A-LIG_?JMXD83VXv#%i+)S6~
zF$dlwn(802EwQDRqTHrlQS&2MtA#GuZFl2HvyNBCe#OIb6Nc%NfwE(nL)vt9@awS$
zY36>5=_tFDITp2M1oM-X(bS6_k#Cd}4-S1YQJm}Y$r0x-cq-|(a!%9F`1-*Hr8qr-
z+$tpNseWz#V^P$)MUQOVuoA>eS^hq6c^(eF=j}6?FDL552Sx=;Ccf%=XX{G#Ch@8x
zAJ&aC1PgP1<+l}yXL2h<29#V9z1-TG8!gxy8j)Mp;V#h|mbtw23!%C#OB7m_Rx92s
zCX9g6WKmFoG!LXKVGVtT<y;99s8S9ENl4Z}b_ysu4!)26x@{~gkJVu~J`z}j=WhJC
zg1;2}pn0g8Q%g<Ib@z6yyebc)@`KZ<why1b*N>YPF#V}<`0g*RZ^`ebBnPA`tY4-@
zim<lS?z@vjHBF0r%yI_JVVY1Lv#`iIrL-d(!+)lt?zN^?u^h%HaLzjUNP#2WErDyF
zvJkacskLA?%BImaN{v{)lO^a_ac7QrhQOoYvvag~4}IUF2Q4(zm&$LuM4|;J;>`sT
zVkXaUdCJTyGP^~7A8@m`N9Xa1zBG*6X|=uE&z}#2*Lyg9PVAL^13}YVf1svGyZU_1
z4<h**+A!ROC5z{h@G4iwKx0v=g^rxSmy@N>;${6#^JbaP$nwwtyx!vMH7gAT+w-5d
z9gr43nP-o!5>6TLjcxla>O8QAZHKLc)Xf8d3g-;T3m(~QJyxYX)yu>07IwP3N1TRY
z6xE=eQUTBfBnrMpEFduys#pZ@LqHX%wM=**D9}8-hxW=GL?8%j{RUJmi(w*w(7bFj
zteAsyRHN9;M8K_gB_&8NSh;q7HVAS<)+jF~0%85dJJbRK54R;!TmV|f+cn4{D8L<V
zDG1PZC0#lJ2cSKDLTLax*yjZNLkp$(18AC$mv;z&_V@9p0Ym~DS{opeU}Jzt0nmB?
z2?Wsk0137<1W2&n6d+>&v;{!M0%$7$Z4Ho#0NM_qz^_~Z3T)yIfEd_@#s*>t0NMwB
zA_HhY00c2mt`a8#fXbBI5jeP8GJxEa75K25n;jh7F*P8taJK|96+`Hc2#8gXgRTu=
ziGLa-9RR*L*BT_49Codnf+c+byU|D)0=SJv$`ru;X{2DXxOEgZm@p3a$HE%G;r~c>
z0PauM2J^>lV7LRg4GbE9Tf?BkAb;Y7TL)tM6G!|8mJfj6!14p|e^!L{2k;wM>yd?9
zsCB9VaFp@ufNbbMT98kG2Y}x|4F>QVsDCnp-&h@-tOEdkV}ZeeuI@B|umKzb5Y~a&
z(8mG725=|<=Xo9L1VGq8ji3bu0C0(9Vx>4XR0psCNP?QMBfMx+5CAt!Okpj)I149;
zU;=T)4{_V!e?0uJg<v8Op|AK1xaqw_K+_2XxEL^N1Mmmf)&&qx0AX62b!%ic>4+3K
z0Go6q5RSrT-C7p@(!meJHCsrXNTL(BdL~ieaQtNf5-?jC5OCmT1{geM3j+*}3>WX;
z{xCQ^2~Gt(=>IJpL%_hhWz$3~feeRwZTuS|43>o1G7pC#Zy}9a+dO~Uf>_GmIy(Fw
z#Nx4AqX4dbWvq?=fC=E*XT>HRb_+0uK!Er2=6R&8c_3m4^ew9+f~3u?F+|+jM*AC>
zNFc$rx=Bagsv}{^n-_;6;jtSS^#4{zB4f5VBySBNhP*Xe7&1IIlc8X@=AD9D+mHW4
zdM*AG{CYXU(tl(qMEDr|--&RV;F+5Vf*9=DnXoqgwFNP_E&f0Z32x&y&!d1C`er4F
zgI#T!2;#_F$bbaO=CFbwc}tE#kg}yd;3`@>^VSBy`R~Rj5OMIQFhz_&*f`0o0s=yS
l^@9n%A$%g>^J?Q`$2BPAL=erB4J44r1U7kjJ@W%>{}1HQVuk<!

diff --git a/notes/images/MC_all.pdf b/notes/images/MC_MC_MI.pdf
similarity index 69%
rename from notes/images/MC_all.pdf
rename to notes/images/MC_MC_MI.pdf
index 6363da3c01c690f3ecce4fa3db7e7ae3fe5b0ea3..3c3b86f50c8ed4a5db6618d6a07b3b225d9782cc 100644
GIT binary patch
delta 4305
zcmZu!c|25IAAd(lk<wx<xe_yzW#-=5DU$5_64gtH!KjSJ*vn{INx~6IMNyVYWhhcg
zmNqh|XjK`~;_+mOq`Y_VJnw_&{&Rl!d(ZiOf4}edoHL&(kHh5iu3BOw$uLO|0d|;B
z-ez0&^72Kuhpu$zgb-_29x?4Jd!%1ChbR{<DUuU;)khuou&B6t`ci*?e_`~^io$0p
zb%&G=YBpq=7gugldqe)MFN`p#`{p(}hWz#!#m0=0n}=1#*oO|=51BMB&3)fk(P@8S
zxl{GWvR^;0cn?RGjVyR<xNjs1O1X=0sv3E1S8X5bF16C;+*r9836r10-O5bruVlR)
zQs5Xr7bcX@ZSN@OpHITIKI(Fjaq8}`s-`^sjqI*zq*o@if!s4J2LB{|;LJ8<`=RQ=
zqc7Dhyji6QG$UCzUKL;e?dY0Y*+xfcSJh{qYc2?51>%CY)^uiW&2+Gnjkm&G`*4i6
zQ}Hml2WWl2nEAekh^OT)uy)3;Rc7i8??(zAov6F~k$uJHt#)<t`rJnw%#>ve&)Nhb
zE$3Vu7N)WW#+(lf6)wHx`xzga6EB(Na3J~w>$1GVoolO^SK=KO>N<)7U9P=bQfMI}
z>H0a`i9Gyp>xq+57FRDV%}GF~Us#;Q;>E4|RTop>D!<nkI-nH2om!)lkv_Yy5ISVF
z(?;);vM^csd%(4;loUzT@~1+^5{fl7Tf=&uWniE0EPc*1ve~^CDYD%UAJh{)tP*eS
zO*@r8d%cZWg+Nl3K`PRC7SPeoPpsFiQE(3X(H5KN?lqQHyeV&!FxfE*32D5(X$w!7
ztSyziV(p9A#Mk6%KZVRN-R_4tm21(-w0b8}kj49hKiOppT;JV+sn*4?zp=Sl#IBYK
z|FUIW`RA_tKMj9O$=1kO>6*6XQkVi)RHy#DEgsxMb5RD<U3!?M+1Ty68hQ=MG+_$k
zqtKeW0M5u!AbYzMS2&<alZ8TFRGj=hexm!cVQQOeO_Ln1=CaF|`7-*P`+4zm!<$uv
z!0}kHYILYbwnVhJ^se84(RKPVVe+?mi@HyH8=I%?|LS6SiI3)*?HVZA*@$ZKz*hHI
z^Sc@%TGu|DreEvs9U14XY|d`Xm5a7Sd&|c>)t5JayUyIKd_%*LH$1ZosIs}S%ep+a
z8Jic)b2}TlU3woQP}}&YO!{h_?gt(_+iRJ&GFv#N+gq%ynp9r!6K#|Tg}(OGY;2B1
z_6}y4NS5m>axT$0{Z;o`wuZ*yg2Ul)lB&w`vAfQZP#=_=<YCiQ*n=2dWseQ9!C#R)
z6k7yj?)$RQX*@=97+Lgi&VaAW-c&0~)WXDs)YSpK#ka8v+=jU?1z$3EeUV7UNhq3y
z9R8gC->jb%ig`of;aR0umNm6D$bOD2-l_X^bs6?$Eo67|Kw@KtibmBFB_^p(zvLo6
z(Wz1W^=P$EzQej0QQD1=`QPh}gI(9L_od0b2sAYR@See&_j~b~zRsw`sxu61t-{>R
z8ddGDEbpsmK0VFP-Me+W0jv6q!93pCI*#X|6ve=VYJixNzZ`9nD8ed#Jpm9Eiv-hz
zm`jg8oWJ>bV3#ni<NT8q65IA&gP(Sc98sn9>Em*luZO)=mLwIg^dl@zzL|w@raF%r
zCSUwMDA9io@YWpD%4)1|=)ElB54@|)xR7+L<Bu1%wiZ=oH01(mnnj(xTHT#vca9#e
z{3cUldVug+qsk73?g|ZM-1r+UPi7nFk}IRvZrwBX{+Oqo+S}WeE9|p++Ke+8l8#GW
z9MdSM{*HSyq7*RPL>10^v)+siUnWL8ueC_(muQjI;?C!}XWH?Ltm|D$&Ev}Jt<(O{
zQm@`vhFQH+*FJG6qT=Y%7`LOa<PWv-jRloaAC*jc*7j|da^?<b7tB5B0XI7@juGh-
z;2+NE>y+X?^nWgqe`y=kpK_zaZRAhK*3=a@i%)xc9nbzzn0hGD7|F@d?!4{uj%dR_
zayq;yP6@A8#`7zl71e!uyOsvDTic-Zu0;M;^qU^8TP^9aTG^c2&S${=S*@`JzdC<m
z&M4M~>j^ar47o=?l=DirHxIWC#y=a3fA8A4hG#W=&OTf>?fHVDHb$22xzPTl&tpBc
zTbYW#TY5TbcBL!p?21So{`7rR=ojTxaZtm)U}Yh&FyZN}uW0!dF{qW;K34gysU8`5
z@I%eT)aVdP%J1+uiQm?`v$lL(vii9D-QKgGH8#0#EK)q2lRnGD>(eKGgCikH1unPh
zRXog(WtEz6?Oa=Kc)styFKZugG5C{XoMs-vwzq%yHrTsdW=n=jXXJ0|d)$`jj|p3b
z*6cZ1VK93oju(P4*d}@=<<d-$lN`uTSKD}ykgIXd*?OD&>7Jr=1H*{pyCQfEK3AmJ
z=Hz)jZ!f-9yX}VQ<trk~wtI}U1?@`V#i*M6!L09TR??*V1RfQJ2T~q4f0Qk0AAa@b
zr)SKv_R+{sM~W&Rj~w)mTvi$+)nDIxjyTt^Bv{tiwQrQC<CYn<rzP1{6cpH>++MnN
z%R2YIY|FE$c!B7aI>ncH>Y$|hsF{IPsZDR5f?UVy6v@1%9rhzSei!t%DRujbD1WMY
zb?>i}@hGd3yk6kn7nLFDcJ1mnX+mH7<{i(PPnBeRI4f&E!taTStfC%YKG^LNM#Z{h
z-dEE#+PBbC#KX>hvdVD%-mcBLPOtY>M2r+y3j1ZLX6*t-a>cyC-SIhYV|#cbL&Nuf
zR*dc5wEokH>s-&>Ve4}Qr#f5wvu@`<?3Io0^4sfQPSEda*?dh#`=i%_4GV2Dl`<0X
zH=1R2^b_pQ|7TOxEq<K^{P6bEML)gDg!j|(%0b+`knDQ_ynrMZbIC@mcmG?Nl;c&8
zIQ6%7n>MvPugiJkMC|h2v&GlWJT@gbEBy+`^|kt?)_v3iFE0oSqXl^#TXR?0B*w={
z*)A_%>u760eQ~<TM_wb;z_??t**Axw2(R-yt$m&ha{Yg7wCl*RJ<H9GercCg6mZI^
z1JH8LlKV@#aJ9S@r6aE*`S|05_==<7Up{&FBT;Xt#WcTC-_xz=w|i4VZ>$w70J7Pi
z-HzQe(pl}y4Svwwb|b@t@1x*%M7^rPzwKSkszl5Nm;F0UVs*C}60xoOMRD$3>z6=X
zs-dZll9nNkOo|N^dnJZkQ`;U^t%<yTxP89zqK>_*zgc|T65(|!)WoNz)7QVpbi<dl
zjq)9pQZa^7A`N5Q^5%gN7oR^H78+GX-I)FE)qCE`z!$#K%pl(oE{Ey8Qyj6k1xmv(
zM1;fvLL$-s2I*kY0xU3}?}C9ruo7NSJKq`uPR+*&=nMQX7%IqGNXF2>(}hPc6yPNn
zwwR12ae@pD|4u*w-Exf>vcQu9O9(zK`GTR)CVAGZfe4lmOdnrCU_&rtA_{>J;<PLX
z1?*D737Qn<V!&Jl9Ka~CK|VT0W92}n@;cC_hzAh75O}Eg1GW+#^r+N<%gQ*xt7Wq>
zKuZM&?&>TSD64v7K&=W+KvZ+Y0Kqc6pb(ma0qx6if;M#p42W8R6DVjXVt|wy9=NWG
z1hEhv#Np+@Tuo<iMI8_JYubRGm3Z(&^B&M!g##U0!9a+P(O5^|u89*2V0U6blO_&b
z<O?daa02<&TQR^HivyY3MDQ4k2LswUzzdHT+|}V=z*jW*ubvVHaJ03-8Nwq_uZ;sq
z!~l??gG0Ls2P8c_xC{4TNT|wc3m}3?fH%j1!}dW1B7=7bOPYpa5S0NDOo1n9xh4aB
z_l6kgyAK2tVTiy)2MR=Bp#vQv?0{gHC=D)<f-wl{l?Cu5Yb_R={!hOU1oev!8p+2n
zFzVUz**vkM;%y<e%LGCw6DFBL2z4T|gb=Z#P7uKsLMF|0fRITuogmVLnFMEuG)^Xs
z)<{AoR`@_9(nQ3BNE2I!;KvF>UnXh%*t}2<M4F^xL8M8l9S~`pN}3P^ktV6wK}?7|
zNfylHu-Lv3d7?c+2t=M_3}J>bgCO!GYX~a>B9HS*6GDQ$eb5Ax#9W9xLAs8>4d5`*
z50lK?a5hAqY%>%hPjZJrlyPQhLIgyaB#mNn*boH)t@>*7C;7|}2&MuFgKoi~fr12B
zWV9AF*WTzPh73L%5dz2zv4G+S3WB7^0~q~}gBCa)qoLJI!${)uGaVD_F&(4AVDCB`
z9zrM1Oicb4j7mk7n39A@M#QIh8b%b0F%=`x#7B85Mj^psjEACscS%I(;*B+RBEk@B
zIu#?45Y)&i7>v-w-b}+N|AG<es55`f6Y~>^Fo}+qF&_TE7z&9xorFlCQm2y;DRlAW
zC4T;~mLcAbQ!y%$F=HN;M-uPqX%ne4F&c43Cqx<<nQk4CMiU>ysrgWcXQ)aiQO7Z~
z|NmD_r~J#}f9aS`qt75_5NE2ufN3*SWgv`y$Up8sk;tGjW~jhGOPtmTOhm-j&6Mt7
zA_<<}9eUl(JOxCC_zIf37Dk-J$8;J<X3R(kQ)l=I!?YQG!ej=W4p<wO|8=uCvDsW`
l@*;3%MKK`?S^TRQ!u95GBREVyX_x|&sM4C6ruJsi{{#1XXU6~l

delta 5177
zcmZuzc_37M_cvLyP8dt5i)e<JxwDToVTMrF9!r!pyNt1<Wki-j(G5j~r?R#&Xwg$D
zB)ce*?Wu^8EENyQ?@oDt?}PXLbI<+U?>V1y&YaJgTg1M15^>uc!B3zPwTaL!L^#{I
zDOl;*hW@rcnr)qY9bNlhp1CHcd4g9cY0DEQJ}04y<{lu(Z0GeO$9A=F{3p^L5Rs2F
zS%)Sr+`P%?_RP*-x+bz1w0jXdH&eO9Slaz#!E^Jx=8t=yg$(;6wGx#7IP$>7pgudn
zeu@yZdFf@iVYuwA$sSqjTlbOEY0TWHA$C2gsmoDt{R5RqP06Vc%ZofMyXC#yvI@gP
z%CO_%MLu`Z3I<vE9>Yirb9xD)%;y4MWR~?p!0dXq-O}0HJ`LUlrt31ETGP_MQfB*y
zD|_7=&N8)L>|5$_{lPDB?&o~rgz20eTTH>j*vJ{J^#@6siuL`cAS2uOy|C#FV9%a(
zcTp1!+MCGu{^`?15JjUqH11lQndg-h#Adq$i_3n+Rdb&+P9FA_j6c4^Tv^UrbnZO9
zY~&>UUQBeJrPz|8!wWamV+cbW@v#SoJWH;sOWoLm_@+>*a@c|C+p|5-K1D;|h-HSQ
zUx5wZx12=B0ZqHRVxmnxV58h5=)IF|Mo8r?8tlQ=m7FajqzU}Fr&J6Sd2x2T0}7pK
zVf9!|pI3(&tFK{RPA|_>wUBo{rEG6O&Qm=gPnG@d1;6Vc$}mr525#R?_emH#uXKo4
z>q49;0yK~C7SUmvj0xwbrLtkGgxBHi=yVxE$!Y&mBj{4-*q&|^?Xgpf-Rt3J#L)I@
zxR<_Fp*_9Sr{{dRRkh4lFKIFVV@(P^#p#_u@!KB?H~E!lrWe!Dfu60gr&A_(W;_q2
zG}~}`4e}ql+(anJqsp%Ct@S@(#GqRk$Cuo7m1|GWDLD4GheRpAzEgmcnJqUnS{;?B
z`~ew9D=}eVCHs%V6|q0)*;lx|4t^4QdT-Ypxo1yp2NJF9fO#b^!+I~1Xpj6m_>RC=
zuP#-^XL}Uk5Kd+yfhJ3y$sCvyOcMIqdK>NqJwij-g<os5L1^2}*y#h$1z7J(G8ujc
z-{A68uf!qqN3+j}ak?z^x?-h0B#Qa<<w-!7Ol>*^?g3ldW~S(MLq(pdhy>MrKn-Sn
zLKUubvXg7JDHZMAch}zCh!O0+km>rWB2P`Wf^2pE`Nr3=rxmwjXwl^eeR&j*x28H{
zU^|0uV5<8fRv(VwW#W-G?&<;eg5NA{v^fwqSQ`q2BH(n?t|ZIdzV{wq$gd>7K}v1I
z8@3(yoLy!mo9V?HjW6Nh2<c7>=-_Sr>4#-M)5`J!=<SzX-DX{6lT5%Z=Dy-jk?Cm#
zKPT#16057Rz1`nNKNSdsUk(&WQLdG)NH)4i6j5rb;mdJcXzLSD<ZDguSW019H%eUN
zb5P3zF*$N^`pSaw^6_P%0S!5Lk+*DPBg@fJkv<f~NZ{}eq^i>2=W2^=8aa8aZ9;oa
zrGK>yiY1R@j7CNLgEu6m%EV8re6TAtuCY)OY&UyZAhW%ke*Ad^UsAwKYm?CAV$Xfm
zROCp0TA22NLGB%cr90_+E4u3so;Dz~ynB*$D5Kpu15GWTNDjr%b$uBql%_kgO}23&
zB-Yo@{U}>HyrlWWdj3|AcFE1Fbi`nLjC4#UxlnR{S6Tr589JXG`EuWp*{Rvw5QB%+
z_s2sX_SnQ!Dw`E)Jw79JdtK{ZuZ@=$1|ljqNbNRxYV<m?x>w%gZA&L_foW%@Bt!3=
z%N@DKnCQ%FR%f<pF+@9yy*d+y+`qEDZ`9#Jt%g$?)c&mRb+E5<@ifcb1jg?_vW2qc
z;66Z_aY@4ORqy+XVa;HJ?$oG}b0evLea;<wn&_%Z{qCRm)-}$oJJNA(sf%P{Y=jbX
z|FF>!n_tDG9^ADtO$+yE;!>RZm5++|XT*kV)5sUY_VeygPDyNieEeeT%=*j1BkW1l
z4?~))bA_!4?X%*SuIvld3vW%zG8|P@#Ob?GKMJisoTBqZQ=Z`W{j{(c!RPUnZ$qNW
zJjOX6pNccsr~M;G(T=u7mvwe>4vz5dc28jZyf+v}D-hJnNP8u&ucmDHzWDk0vpJ4F
zD)%#HV(ET#>Gd(BJg$ChYrSkxbMh9pnPlPE!DI5DtEwJwe!g4K!P3|dkDl-OQhrx0
zLPFATuoe4Pz3e~$wEz=IWW5(N-l=wCxaBB&Dq;T1u{U#a0%v{FrVP)y*$bAa@v^fQ
zCG8k=p}c=~>^HDVOinwC$6sjZWL>H>*2wm2bJ|#7wBD%FAi!_O(9ZOvxU#ZA1NPJY
zGd+73?j_NLp1L%5PXC?W0!(My)3Y!)46L3N2NWf%WVG%N$SMt8cXU{h?DtfdQ^<Lx
zzQaD>HQ=))xj5sLx8bu91O35?^6(Fjj`kSMw~}g%6+2jMkL%8iE4KUtpqoaISFDdd
zyJYLYPpa!_y!SPw;Cg6uw!!Okr$I}0Zc{>)vs;^2TrQ{o`FqutY8KNX^1O44LO+yj
zctfIE^XbABP2<Vdy8(TUpRP#UpuhWBm{UJAAs>=86`LMbQH=GiLjwAz#wuc}KC`PA
zCTDIf<``mU>gr01va*^?KiyYsCb)7Ue2X^tqKxC(SyuToC2nu0PW7cqZQFT^y-Qxi
zQB^O0i>Q2WNu8L{Y(Yv+un{h7f6rN}+XHRn{nXei*3q1(_5#zBO)`6i>c?69UwC%N
zHfxcBed$b1LMY1P@gagElRry;&?)d`vzb$uVdb0D0oi@pvvr<=d2Wi;F)fpInJ)#a
zkGCcU)?1=GN2JK{r)uyIALZLcv1J+ix^?|2mp{I_`1<_$yDmzbUkXJp4oXm#CfG{{
zXLHY3AoaoS&$_YwuUz(iPkOEAkUd|H^14vt_vXO!IlWKeKG;sDP`P)#<$@P;kecd$
zp<3pd#qSHOUzA8s*%%#3!^b*&Jv!9V<t!=K@rgKEqob>NSLlhgtCY@jC2Ax4!lr#S
zv3Y6q8<7Dj`nB}uCm9{;M~6$>LUe<8ZS$)>@48l1R53O{9@EP^HY~kugECRPKe!IJ
ze`9dy*1dz~d8r}6xvZ-jl*V>3_TU48zxUO;bcX1X-0&K8cKz>Ej_O?{s@FbiGp#8q
zefQ4BO@AXq@k4FCdqTya3lwUbWFMht&5Z4kj=xxy_(x*!<Vf`X+j4LF_}1S}*z#vf
zX{ya%sK!xR?FrT9hVxH2sV<J8GM=8&zNVyG&l8F!<a9+mu6jOFe%VzU@7viejmQl;
z@$pXbj^qi25sxlg+3q)_EMV$IcT~B;HQv-w>qwp3FRHGBA5{8Xj=e7Y)=`qi4j(BV
z{2+G}(Eb)_+t|=qb7W?5q5I1&$!}8&h|O2!Ca;a&+Sntfg^#n)sk#x^gPnLw{-ics
zKJalE)&C~~DJd}TDn&W;;4QjCFsMGby}egDdrHk~Yhs2|f|yeb>CBg@@mSRKSN|?E
zx^qCRL8HL+DjqpR?KS(aJjoiR_YO3ex&0N8WFg(dpA{8)wz~oEa+TL%)a-x`uC;2&
zr`Dmo^80<$cbN{TcMqiub8~kRbY4!~aVC1OqmF+2xF}f`QoV$5PMW<W>>cbC#$tNA
z2XU<s4}W101aJV?z>`4|bW8+<@FMomyeNiqOT-d^03mSGR*tHe4+24fPKgr{ROpd-
zGJ*_QNJI#cU{_nfzyOL7M?pwX3lc$QAxNAS5M0OUlKg~#ghfQ9%>j`1>l`)&z~g>f
z&?HV8+9HGYrcz;d;pL2?E|OBJ#cFUuCWi~|ZkVIS%$~GOUwngf(Y|~5h13UHO@o01
zW3*Pw%lPktD05oAfeq0)M5oFu?9l@+f-Gv0(04Tc^GJ*r(fqw~$1ro)2@`FQ)vIs}
z+8`qZRm-3`3aIr6h$xGOEM-EVPFXY?Pb9A3x3V!0uEmK~FQ9(Z@MK=EJGA@M<=&-z
z7CLg45(#onZzUtQs+Uie$KRLHQQSwK-q|DiV^hUQZCes*se`O?<A`<8am!kh^mfyJ
z-#*p;Nez>QbYZvR0nZjj-wk5OSWn&00m05}M#**IsWi_pFSe?mtCZ=p+4kr*u?|ce
z-+bQ*MdUSFH0tcnUWx9?j`3H$&6kVBH(iiy+}h&ryQ}cLGADo2onpu5=)NMn#)O@`
zpMJv|Y7@_{(^8GT`mf8flcw!*_wI1BtX*F{KJJ)gwJlKb#WDN&+Uc1aly@Nom6<#0
z`A@b>XAQJPo%v$O*YD_30o|0-g%(ti(2UG>Xj~ox<!MPk0kQ|6L<J1AAX^V1lrWr+
za=Zu#w+Rhx$4GLvZgfXL!JE-`6x@pL5jZ#-ob_!U4?-A$KTOpkUfNh&RfLeYT|0~4
z5lC!mZ$~!5Q)Y`J68Xz5TT=5qaH-p!EYgf<){C^gemuMFC9zNrpv$?fV1t0HmC+oG
zk{|-=gAK<`G6={<1<lFaEQ^42R58#i;1KizzOPqDLP^SYki`}Z#HV5j4XI(ELn=?8
zAayk4sd^Ba*@}j%vw<wpXbx*j5CW3JpkXI3h@*k#WT?9#ps$)5kP=!O`mhZHvCvnb
zNG%L!JBEpXWVAIneHyX|1PSf|8dBWW{Hp~>EngUcWdV085>^7pZHx9aG)5AIGQ_3D
z0j#mFJJXsO;u!|uNzi_*pDPK#8UbVq{01P?0IVrM0Rb!npb!A8IY5DND}agvu+{(-
z4`6KptQ|ll1K2$P4SwYg&|nu&0K~yQ-ohXb4(10xQ2}fq0D?FGlExiG5MV)?K`l6S
zN&K=X|08fb0pbUrA@Hy)75Ug)$#T`Mr^ihRPyVe(V*pS2ZJ7di>Tiny;Foo44iJ`g
zYXuOv>jKWu8X){45{A<y5PoBJ0D<`1+5-@N?*=$~!ZL>^K%fBZU(VhDVVT6oF9Lq}
z#UqUM2?=Kc#APZ!fVfN*2oQhQgbfCW%Tyu3-T<uKU(pVFGyOum0OBu3VQd&cTxJaO
zj_`(EmMaVMivoyl%e-L+-96z5%LTFk(lRT{m+9>dkbW@>V_A_Q0BM;!93U-oM*yT{
z<|u%)TxPU46Fy3Ywo*2WF2AJX0gwXOQbq-^R<KFrkf;jl*1>m6nh|%Ip4Q2RmnlXU
z>q`c?xpOZd7!dwtfFv9RPIVOq;^D+rVR$l})@qCj=dcPR;7Cx7t|gm*qpl4?Ai<qj
z<v{))WK`;EG8~b_{n1zW;)n$9V6DUmG<cYCI5v#_f0Z~AiOL06U<42cw|o_b=g#LU
zj0EdqB?c16<Tb?+i2ug8{5Twlr@^XQ<xHf2Ybqg;xx>l5{*^Rk4Ka=k;@0?*@!aWH
z=}TslC~K-ClWA+(N~UrLZ>2K@M}x<76-M|U7!_U>tIyH6()rc@|5idJaQ~232Bebl
zYm%YTxLRC!js|b?H5hTNFAZ+tzs~*cKaC8pyVU_{YuXRuaOBm41h2t0GXw8L+UjnB
zI0|<|tw<UKK{kFZi06I-tvW~qSO0b(NL@44upqhXWu-5kKwA?JPvR~__)`7vAc1%?
zZS~NAWIU`>M?=N`Y%SZ65EcN3<L+ELzi4<jEq^h&Gg(ngZy#ZhOdt`2)znOP(}n*J
DawxAt

diff --git a/notes/images/MC_few.pdf b/notes/images/MC_MI_VE.pdf
similarity index 100%
rename from notes/images/MC_few.pdf
rename to notes/images/MC_MI_VE.pdf
diff --git a/notes/sections/4.md b/notes/sections/4.md
index f1972dc..6c5f081 100644
--- a/notes/sections/4.md
+++ b/notes/sections/4.md
@@ -209,7 +209,8 @@ $$
   P^{\text{oss}}_{\text{max}} = 10.005 \pm 0.018 \with \chi_r^2 = 0.071
 $$
 
-where $\chi_r^2$ is the reduced $\chi^2$, proving a good convergence.  
+where $\chi_r^2$ is the $\chi^2$ per degree of freedom, proving a good
+convergence.  
 In order to compare $P^{\text{oss}}_{\text{max}}$ with the expected value
 $P_{\text{max}} = 10$, the following compatibility t-test was applied:
 
diff --git a/notes/sections/5.md b/notes/sections/5.md
index 7bfc1f2..cba8c15 100644
--- a/notes/sections/5.md
+++ b/notes/sections/5.md
@@ -1,10 +1,7 @@
-# Exercize 5
+# Exercise 5
 
-The following integral must be evaluated:
-
-$$
-  I = \int\limits_0^1 dx \, e^x
-$$
+The following integral must be evaluated comparing different Monte Carlo
+techniques.
 
 \begin{figure}
 \hypertarget{fig:exp}{%
@@ -24,6 +21,8 @@ $$
   % Plot
   \draw [domain=-2:7, smooth, variable=\x,
          cyclamen, ultra thick] plot ({\x},{exp(\x/5)});
+  % Equation
+  \node [above] at (2.5, 2.5) {$I = \int\limits_0^1 dx \, e^x$};
 \end{tikzpicture}
 \caption{Plot of the integral to be evaluated.}
 }
@@ -33,14 +32,14 @@ whose exact value is 1.7182818285...
 
 The three most popular Monte Carlo (MC) methods where applied: plain MC, Miser
 and Vegas. Besides this popularity fact, these three method were chosen for
-being the only ones implemented in the GSL library.
+being implemented in the GSL libraries `gsl_monte_plain`, `gsl_monte_miser` and
+`gsl_monte_vegas` respectively.
 
 
 ## Plain Monte Carlo
 
-When the integral $I$ over a $n-$dimensional space $\Omega$ of volume $V$ of a
+When an integral $I$ over a $n-$dimensional space $\Omega$ of volume $V$ of a
 function $f$ must be evaluated, that is:
-
 $$
   I = \int\limits_{\Omega} dx \, f(x)
   \with V = \int\limits_{\Omega} dx
@@ -48,14 +47,12 @@ $$
 
 the simplest MC method approach is to sample $N$ points $x_i$ evenly distributed
 in $V$ and approx $I$ as:
-
 $$
   I \sim I_N = \frac{V}{N} \sum_{i=1}^N f(x_i) = V \cdot \langle f \rangle
 $$
 
-with $I_N \rightarrow I$ for $N \rightarrow + \infty$ for the law of large
-numbers. Hence, the sample variance can be extimated by the sample variance:
-
+with $I_N \rightarrow I$ for $N \rightarrow + \infty$ because of the law of
+large numbers, whereas the sample variance can be estimated as:
 $$
   \sigma^2_f = \frac{1}{N - 1} \sum_{i = 1}^N \left( f(x_i) - \langle f
   \rangle \right)^2 \et \sigma^2_I = \frac{V^2}{N^2} \sum_{i = 1}^N
@@ -67,31 +64,55 @@ Unlike in deterministic methods, the estimate of the error is not a strict error
 bound: random sampling may not uncover all the important features of the
 integrand and this can result in an underestimate of the error.
 
-In this case, $f(x) = e^{x}$ and $\Omega = [0,1]$.
+In this case $f(x) = e^{x}$ and $\Omega = [0,1]$, hence $V = 1$.
 
 Since the proximity of $I_N$ to $I$ is related to $N$, the accuracy of the
-method is determined by how many points are generated, namely how many function
-calls are exectuted when the method is implemented. In @tbl:MC, the obtained
-results and errors $\sigma$ are shown. The estimated integrals for different
-numbers of calls are compared to the expected value $I$ and the difference
-'diff' between them is given.  
-As can be seen, the MC method tends to underestimate the error for scarse
-function calls. As previously stated, the higher the number of function calls,
-the better the estimation of $I$. A further observation regards the fact that,
-even with $50'000'000$ calls, the $I^{\text{oss}}$ still differs from $I$ at
-the fifth decimal digit.
+method lies in how many points are generated, namely how many function calls
+are executed when the reiterative method is implemented. In @fig:MC and
+@fig:MI, results obtained with the plain MC method are shown in red. In
+@tbl:MC, some of them are listed: the estimated integrals $I^{\text{oss}}$ are
+compared to the expected value $I$ and the differences diff between them are
+given.  
 
--------------------------------------------------------------------------
-                   500'000 calls     5'000'000 calls    50'000'000 calls
------------------ ----------------- ------------------ ------------------
-$I^{\text{oss}}$     1.7166435813      1.7181231109       1.7183387184
+![Estimated values of $I$ obatined by Plain MC technique with different
+number of function calls; logarithmic scale; errorbars showing their
+estimated uncertainties. As can be seen, the process does a sort o seesaw
+around the correct value.](images/MC_MC.pdf){#fig:MC}
 
-$\sigma$             0.0006955691      0.0002200309       0.0000695809
+---------------------------------------------------------------------------
+calls               $I^{\text{oss}}$        $\sigma$            diff
+------------------ ------------------ ------------------ ------------------
+\SI{5e5}{}            1.7166435813       0.0006955691       0.0016382472
 
-diff                 0.0016382472      0.0001587176       0.0000568899
--------------------------------------------------------------------------
+\SI{5e6}{}            1.7181231109       0.0002200309       0.0001587176
 
-Table: MC results with different numbers of function calls. {#tbl:MC}
+\SI{5e7}{}            1.7183387184       0.0000695809       0.0000568899
+---------------------------------------------------------------------------
+
+Table: Some MC results with three different numbers of function calls.
+       Differences between computed and exact values are given in
+       diff. {#tbl:MC}
+
+As can be seen, $\sigma$ is always of the same order of magnitude of diff,
+except for very low numbers of function calls. Even with \SI{5e7}{} calls,
+$I^{\text{oss}}$ still differs from $I$ at the fifth decimal digit, meaning that
+this method shows a really slow convergence. In fact, since the $\sigma$s
+dependence on the number $C$ of function calls is confirmed:
+$$
+  \begin{cases}
+    \sigma_1 = \SI{6.95569e-4}{} \longleftrightarrow C_1 = \SI{5e6}{} \\
+    \sigma_2 = \SI{6.95809e-5}{} \longleftrightarrow C_1 = \SI{5e8}{}
+  \end{cases}
+$$
+
+$$
+  \thus
+  \frac{\sigma_1}{\sigma_2} = 9.9965508 \sim 10 = \sqrt{\frac{C_2}{C_1}}
+$$
+
+if an error of $\sim 1^{-n}$ is required, a number $\propto 10^{2n}$ of
+function calls must be executed, meaning that for $\sigma \sim 1^{-10}
+\rightarrow C = \SI{1e20}{}$, which would be impractical.
 
 
 ## Stratified sampling
@@ -101,12 +122,12 @@ partitioned into subpopulations. Stratification, indeed, is the process of
 dividing the primary sample into subgroups (strata) before sampling random 
 within each stratum.  
 Given the mean $\bar{x}_i$ and variance ${\sigma^2_x}_i$ of an entity $x$
-sorted with simple random sampling in each strata, such as:
-
+sorted with simple random sampling in the $i^{\text{th}}$ strata, such as:
 $$
   \bar{x}_i = \frac{1}{n_i} \sum_j x_j
 $$
 
+and:
 $$
   \sigma_i^2 = \frac{1}{n_i - 1} \sum_j \left( x_j - \bar{x}_i \right)^2
   \thus
@@ -115,13 +136,12 @@ $$
 
 where:
 
-  - $j$ runs over the points $x_j$ sampled in the $i^{\text{th}}$ stratum
-  - $n_i$ is the number of points sorted in it
-  - $\sigma_i^2$ is the variance associated with the $j^{\text{th}}$ point
+  - $j$ runs over the points $x_j$ sampled in the $i^{\text{th}}$ stratum,
+  - $n_i$ is the number of points sorted in it,
+  - $\sigma_i^2$ is the variance associated with the $j^{\text{th}}$ point.
 
-then the mean $\bar{x}$ and variance $\sigma_x^2$ estimated with stratified
+then, the mean $\bar{x}$ and variance $\sigma_x^2$ estimated with stratified
 sampling for the whole population are:
-
 $$
   \bar{x} = \frac{1}{N} \sum_i N_i \bar{x}_i \et
   \sigma_x^2 = \sum_i \left( \frac{N_i}{N} \right)^2 {\sigma_x}^2_i
@@ -138,7 +158,7 @@ result will have a smaller error in estimation with respect to the one otherwise
 obtained with simple sampling.  
 For this reason, stratified sampling is used as a method of variance reduction
 when MC methods are used to estimate population statistics from a known
-population.
+population [@ridder17].
 
 
 ### MISER
@@ -147,36 +167,32 @@ The MISER technique aims to reduce the integration error through the use of
 recursive stratified sampling.  
 As stated before, according to the law of large numbers, for a large number of
 extracted points, the estimation of the integral $I$ can be computed as:
-
 $$
   I= V \cdot \langle f \rangle
 $$
 
-
 Since $V$ is known (in this case, $V = 1$), it is sufficient to estimate
 $\langle f \rangle$.
 
 Consider two disjoint regions $a$ and $b$, such that $a \cup b = \Omega$, in
-which $n_a$ and $n_b$ points were uniformely sampled. Given the Monte Carlo
-estimates of the means $\langle f \rangle_a$ and $\langle f \rangle_b$ of those
-points and their variances $\sigma_a^2$ and $\sigma_b^2$, if the weights $N_a$
-and $N_b$ of $\langle f \rangle_a$ and $\langle f \rangle_b$ are chosen unitary,
-then the variance $\sigma^2$ of the combined estimate $\langle f \rangle$:
-
+which $n_a$ and $n_b$ points are respectively uniformly sampled. Given the
+Monte Carlo estimates of the means $\langle f \rangle_a$ and $\langle f
+\rangle_b$ of those points and their variances $\sigma_a^2$ and $\sigma_b^2$, if
+the weights $N_a$ and $N_b$ of $\langle f \rangle_a$ and $\langle f \rangle_b$
+are chosen unitary, then the variance $\sigma^2$ of the combined estimate
+$\langle f \rangle$:
 $$
   \langle f \rangle = \frac{1}{2} \left( \langle f \rangle_a
                     + \langle f \rangle_b \right)
 $$
 
 is given by:
-
 $$
   \sigma^2 = \frac{\sigma_a^2}{4n_a} + \frac{\sigma_b^2}{4n_b}
 $$
 
 It can be shown that this variance is minimized by distributing the points such
 that:
-
 $$
   \frac{n_a}{n_a + n_b} = \frac{\sigma_a}{\sigma_a + \sigma_b}
 $$
@@ -184,43 +200,48 @@ $$
 Hence, the smallest error estimate is obtained by allocating sample points in
 proportion to the standard deviation of the function in each sub-region.  
 The whole integral estimate and its variance are therefore given by:
-
 $$
   I = V \cdot \langle f \rangle \et \sigma_I^2 = V^2 \cdot \sigma^2
 $$
 
-When implemented, MISER is in fact a recursive method. With a given step, all
-the possible bisections are tested and the one which minimizes the combined
-variance of the two sub-regions is selected. The variance in the sub-regions is
-estimated with a fraction of the total number of available points. The remaining
-sample points are allocated to the sub-regions using the formula for $n_a$ and
-$n_b$, once the variances are computed.  
+When implemented, MISER is in fact a recursive method. First, all the possible
+bisections of $\Omega$ are tested and the one which minimizes the combined
+variance of the two sub-regions is selected. In order to speed up the
+algorithm, the variance in the sub-regions is estimated with a fraction of the
+total number of available points (function calls) which is default set to 0.1.
+The remaining points are allocated to the sub-regions using the formula for
+$n_a$ and $n_b$, once the variances are computed.  
 The same procedure is then repeated recursively for each of the two half-spaces
-from the best bisection. At each recursion step, the integral and the error are
-estimated using a plain Monte Carlo algorithm. After a given number of calls,
-the final individual values and their error estimates are then combined upwards
-to give an overall result and an estimate of its error.
+from the best bisection. When less than 32*16 (default option) function calls
+are available in a single section, the integral and the error in this section
+are estimated using the plain Monte Carlo algorithm.  
+The final individual values and their error estimates are then combined upwards
+to give an overall result and an estimate of its error [@sayah19].
 
-Results for this particular sample are shown in @tbl:MISER.
+![Estimations $I^{\text{oss}}$ of the integral $I$ obtained for the three
+  implemented method for different values of function calls. Errorbars
+  showing their estimated uncertainties.](images/MC_MC_MI.pdf){#fig:MI}
 
--------------------------------------------------------------------------
-                   500'000 calls     5'000'000 calls    50'000'000 calls
------------------ ----------------- ------------------ ------------------
-$I^{\text{oss}}$     1.7182850738      1.7182819143       1.7182818221
+Results for this particular sample are shown in black in @fig:MI and some of
+them are listed in @tbl:MI. Except for the first very little number of calls,
+the improvement with respect to the Plain MC technique (in red) is appreciable.
 
-$\sigma$             0.0000021829      0.0000001024       0.0000000049
+---------------------------------------------------------------------------
+calls               $I^{\text{oss}}$        $\sigma$            diff
+------------------ ------------------ ------------------ ------------------
+\SI{5e5}{}            1.7182850738       0.0000021829       0.0000032453
 
-diff                 0.0000032453      0.0000000858       000000000064
--------------------------------------------------------------------------
+\SI{5e6}{}            1.7182819143       0.0000001024       0.0000000858
 
-Table: MISER results with different numbers of function calls. Be careful:
-       while in @tbl:MC the number of function calls stands for the number of
-       total sampled poins, in this case it stands for the times each section
-       is divided into subsections. {#tbl:MISER}
+\SI{5e7}{}            1.7182818221       0.0000000049       0.0000000064
+---------------------------------------------------------------------------
 
-This time the error, altough it lies always in the same order of magnitude of
-diff, seems to seesaw around the correct value, which is much more closer to
-the expected one.
+Table: MISER results with different numbers of function calls. Differences
+       between computed and exact values are given in diff. {#tbl:MI}
+
+This time the convergence is much faster: already with 500'000 number of points,
+the correct results differs from the computed one at the fifth decimal digit.
+Once again, the error lies always in the same order of magnitude of diff.
 
 
 ## Importance sampling
@@ -230,16 +251,14 @@ probability distribution $f$ itself, so that the points cluster in the regions
 that make the largest contribution to the integral.
 
 Remind that $I = V \cdot \langle f \rangle$ and therefore only $\langle f
-\rangle$ must be estimated. Then, consider a sample of $n$ points {$x_i$}
-generated according to a probability distribition function $P$ which gives
-thereby the following expected value:
-
+\rangle$ must be estimated. Consider a sample of $n$ points {$x_i$} generated
+according to a probability distribution function $P$ which gives thereby the
+following expected value:
 $$
   E [x, P] = \frac{1}{n} \sum_i x_i
 $$
 
 with variance:
-
 $$
   \sigma^2 [E, P] = \frac{\sigma^2 [x, P]}{n}
   \with \sigma^2 [x, P] = \frac{1}{n -1} \sum_i \left( x_i - E [x, P] \right)^2
@@ -252,13 +271,11 @@ are evenly distributed in $\Omega$. The idea is to sample points from a
 different distribution to lower the variance of $E[x, P]$, which results in
 lowering $\sigma^2 [x, P]$. This is accomplished by choosing a random variable
 $y$ and defining a new probability $P^{(y)}$ in order to satisfy:
-
 $$
   E [x, P] = E \left[ \frac{x}{y}, P^{(y)} \right]
 $$
 
 which is to say:
-
 $$
   I = \int \limits_{\Omega} dx f(x) = 
       \int \limits_{\Omega} dx \, \frac{f(x)}{g(x)} \, g(x)=
@@ -266,7 +283,6 @@ $$
 $$
 
 where $E \, \longleftrightarrow \, I$ and:
-
 $$
   \begin{cases}
     f(x) \, \longleftrightarrow \, x              \\
@@ -274,42 +290,32 @@ $$
   \end{cases}
   \et
   \begin{cases}
-    w(x) \, \longleftrightarrow \, \frac{x}{y}    \\
-    g(x) \, \longleftrightarrow \, y = P^{(y)}
+    g(x) \, \longleftrightarrow \, y = P^{(y)}    \\
+    w(x) = \frac{f(x)}{g(x)} \, \longleftrightarrow \, \frac{x}{y}
   \end{cases}
 $$
 
 Where the symbol $\longleftrightarrow$ points out the connection between the
-variables. This new estimate is better than the former if:
-
+variables. The best variable $y$ would be:
 $$
-  \sigma^2 \left[ \frac{x}{y}, P^{(y)} \right] < \sigma^2 [x, P]
+  y = \frac{x}{E [x, P]} \, \longleftrightarrow \, \frac{f(x)}{I}
+  \thus \frac{x}{y} = E [x, P]
 $$
 
-The best variable $y$ would be:
-
-$$
-  y^{\star} = \frac{x}{E [x, P]} \, \longleftrightarrow \, \frac{f(x)}{I}
-  \thus \frac{x}{y^{\star}} = E [x, P]
-$$
-
-and even a single sample under $P^{(y^{\star})}$ would be sufficient to give its
-value. Obviously, it is not possible to take exactly this choice, since $E [x,
-P]$ is not given a priori.  
-However, this gives an insight into what importance sampling does. In fact,
-given that:
-
+and even a single sample under $P^{(y)}$ would be sufficient to give its value:
+it can be shown that under this probability distribution the variance vanishes.
+Obviously, it is not possible to take exactly this choice, since $E [x, P]$ is
+not given a priori. However, this gives an insight into what importance sampling
+does. In fact, given that:
 $$
   E [x, P] = \int \limits_{a = - \infty}^{a = + \infty}
-             a P(x \in [a, a + da])
+             da \, a P(x \in [a, a + da])
 $$
 
-the best probability change $P^{(y^{\star})}$ redistributes the law of $x$ so
-that its samples frequencies are sorted directly according to their weights in
-$E[x, P]$, namely:
-
+the best probability $P^{(y)}$ redistributes the law of $x$ so that its samples
+frequencies are sorted directly according to their weights in $E[x, P]$, namely:
 $$
-  P^{(y^{\star})}(x \in [a, a + da]) = \frac{1}{E [x, P]} a P (x \in [a, a + da])
+  P^{(y}(x \in [a, a + da]) = \frac{1}{E [x, P]} a P (x \in [a, a + da])
 $$
 
 In conclusion, since certain values of $x$ have more impact on $E [x, P]$ than
@@ -320,37 +326,82 @@ frequently. As a consequence, the estimator variance will be reduced.
 ### VEGAS
 
 
-The VEGAS algorithm is based on importance sampling. It aims to reduce the
-integration error by concentrating points in the regions that make the largest
-contribution to the integral.
+The VEGAS algorithm of Lepage is based on importance sampling. It aims to
+reduce the integration error by concentrating points in the regions that make
+the largest contribution to the integral.
 
 As stated before, in practice it is impossible to sample points from the best
-distribution $P^{(y^{\star})}$: only a good approximation can be achieved. In
-GSL, the VEGAS algorithm approximates the distribution by histogramming the
-function $f$ in different subregions. Each histogram is used to define a
-sampling distribution for the next pass, which consists in doing the same thing
-recorsively: this procedure converges asymptotically to the desired
-distribution. It follows that a better estimation is achieved with a greater
-number of function calls.  
-The integration uses a fixed number of function calls. The result and its
-error estimate are based on a weighted average of independent samples, as for
-MISER.  
-For this particular sample, results are shown in @tbl:VEGAS.
+distribution $P^{(y)}$: only a good approximation can be achieved. In GSL, the
+VEGAS algorithm approximates the distribution by histogramming the function $f$
+in different subregions with a reiterative method [@lepage78], namely:
 
--------------------------------------------------------------------------
-                   500'000 calls     5'000'000 calls    50'000'000 calls
------------------ ----------------- ------------------ ------------------
-$I^{\text{oss}}$     1.7182818354      1.7182818289       1.7182818285
+  - a fixed number of points (function calls) is evenly generated in the whole
+    region;
+  - the volume $V$ is divided into $N$ intervals with width $\Delta x_i =
+    \Delta x \, \forall \, i$, where $N$ is limited by the computer storage
+    space available and must be held constant from iteration to iteration.
+    Default $N = 50$;
+  - each interval is then divided into $m_i + 1$ subintervals, where:
+      $$
+        m_i = K \frac{\bar{f}_i \Delta x_i}{\sum_j \bar{f}_j \Delta x_j}
+      $$
 
-$\sigma$             0.0000000137      0.0000000004       0.0000000000
+    where $j$ runs over all the intervals and $\bar{f}_i$ is the average value
+    of $f$ in the interval. Hence, $m_i$ is therefore a measure of the
+    "importance" of the interval with respect to the others: the higher
+    $\bar{f}_i$, the higher $m_i$. The stiffness $K$ is default set to 1.5;
+  - as it is desirable to restore the number of intervals to its original value
+    $N$, groups of the new intervals must be amalgamated into larger intervals,
+    the number of subintervals in each group being constant. The net effect is
+    to alter the intervals sizes, while keeping the total number constant, so
+    that the smallest intervals occur where $f$ is largest;
+  - the new grid is used and further refined in subsequent iterations until the
+    optimal grid has been obtained. By default, the number of iterations is set
+    to 5.
 
-diff                 0.0000000069      0.0000000004       0.0000000000
--------------------------------------------------------------------------
+At the end, a cumulative estimate of the integral $I$ and its error $\sigma_I$
+are made based on weighted average:
+$$
+  I = \sigma_I^2 \sum_i \frac{I_i}{\sigma_i^2}
+  \with
+  \sigma_I^2 = \left( \sum_i \frac{1}{\sigma_i^2} \right)^{-1}
+$$
 
-Table: VEGAS results with different numbers of
+where $I_i$ and $\sigma_i$ are are the integral and standard deviation
+estimated in each interval $\Delta x_i$ of the last iteration using the plain
+MC technique.  
+For the results to be reliable, the chi-squared per degree of freedom
+$\chi_r^2$ must be consistent with 1. While performing the reiterations, if
+the $\chi_r^2$ value exceed 1.5, the cycle breaks.
+
+Clearly, once again a better estimation is achieved with a greater number of
+function calls. For this particular sample, the most accurate results are shown
+in @fig:MI_VE and some of them are listed in @tbl:VEGAS.
+
+----------------------------------------------------------------------------------------------
+calls               $I^{\text{oss}}$        $\sigma$            diff            $\chi_r^2$
+------------------ ------------------ ------------------ ------------------ ------------------
+\SI{5e5}{}            1.7182818281       0.0000000012       0.0000000004           1.457
+                                                                                                   
+\SI{5e6}{}            1.7182818284       0.0000000000       0.0000000001           0.632
+                                                                                                   
+\SI{5e7}{}            1.7182818285       0.0000000000       0.0000000000           0.884
+----------------------------------------------------------------------------------------------
+
+Table: Some VEGAS results with different numbers of
        function calls. {#tbl:VEGAS}
 
-This time, the error estimation is notably close to diff for each number of
-function calls, meaning that the estimation of both the integral and its
-error turn out to be very accurate, much more than the ones obtained with
-both plain Monte Carlo method and stratified sampling.
+As can be appreciated in @fig:MI_VE, the VEGAS algorithm manages to compute
+the integral value in a most accurate way with respect to MISER. The $\chi_r^2$
+turns out to be close enough to 1 to guarantee a good estimation of $I$,
+goodness which is also confirmed by the very small difference between estimation
+and exact value, as shown in @tbl:VEGAS: with a number of \SI{5e7}{} of function
+calls, the difference is smaller than \SI{1e-10}{}.
+
+![Only the most accurate results are shown in order to stress the
+  differences between VEGAS (in gray) and MISER (in black) methods
+  results.](images/MC_MI_VE.pdf){#fig:MI_VE}
+
+In conclusion, between a plain Monte Carlo technique, stratified sampling and
+importance sampling, the last one turned out to be the most powerful mean to
+obtain a good estimation of the integrand.
diff --git a/notes/todo b/notes/todo
index 2742125..95a722e 100644
--- a/notes/todo
+++ b/notes/todo
@@ -1,5 +1,3 @@
 - cambiare simbolo convoluzione
 - aggiungere citazioni e referenze
-- rifare grafici senza bordino
 - leggere l'articolo di Lucy
-- riscrivere il readme del 5