From 1eae6ac7dc189236e763b111075c8fba20fb27c3 Mon Sep 17 00:00:00 2001
From: Juan Marín Noguera <juan@mnpi.eu>
Date: Sat, 29 Mar 2025 22:32:55 +0100
Subject: 3.4.1 Numerical Distributions

---
 3.4.1.lyx | 1677 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 1677 insertions(+)
 create mode 100644 3.4.1.lyx

(limited to '3.4.1.lyx')

diff --git a/3.4.1.lyx b/3.4.1.lyx
new file mode 100644
index 0000000..2cd7261
--- /dev/null
+++ b/3.4.1.lyx
@@ -0,0 +1,1677 @@
+#LyX 2.4 created this file. For more info see https://www.lyx.org/
+\lyxformat 620
+\begin_document
+\begin_header
+\save_transient_properties true
+\origin unavailable
+\textclass book
+\begin_preamble
+\input defs
+\end_preamble
+\use_default_options true
+\maintain_unincluded_children no
+\language english
+\language_package default
+\inputencoding utf8
+\fontencoding auto
+\font_roman "default" "default"
+\font_sans "default" "default"
+\font_typewriter "default" "default"
+\font_math "auto" "auto"
+\font_default_family default
+\use_non_tex_fonts false
+\font_sc false
+\font_roman_osf false
+\font_sans_osf false
+\font_typewriter_osf false
+\font_sf_scale 100 100
+\font_tt_scale 100 100
+\use_microtype false
+\use_dash_ligatures true
+\graphics default
+\default_output_format default
+\output_sync 0
+\bibtex_command default
+\index_command default
+\float_placement class
+\float_alignment class
+\paperfontsize default
+\spacing single
+\use_hyperref false
+\papersize default
+\use_geometry false
+\use_package amsmath 1
+\use_package amssymb 1
+\use_package cancel 1
+\use_package esint 1
+\use_package mathdots 1
+\use_package mathtools 1
+\use_package mhchem 1
+\use_package stackrel 1
+\use_package stmaryrd 1
+\use_package undertilde 1
+\cite_engine basic
+\cite_engine_type default
+\biblio_style plain
+\use_bibtopic false
+\use_indices false
+\paperorientation portrait
+\suppress_date false
+\justification true
+\use_refstyle 1
+\use_formatted_ref 0
+\use_minted 0
+\use_lineno 0
+\index Index
+\shortcut idx
+\color #008000
+\end_index
+\secnumdepth 3
+\tocdepth 3
+\paragraph_separation indent
+\paragraph_indentation default
+\is_math_indent 0
+\math_numbering_side default
+\quotes_style english
+\dynamic_quotes 0
+\papercolumns 1
+\papersides 1
+\paperpagestyle default
+\tablestyle default
+\tracking_changes false
+\output_changes false
+\change_bars false
+\postpone_fragile_content false
+\html_math_output 0
+\html_css_as_file 0
+\html_be_strict false
+\docbook_table_output 0
+\docbook_mathml_prefix 1
+\end_header
+
+\begin_body
+
+\begin_layout Standard
+\begin_inset ERT
+status open
+
+\begin_layout Plain Layout
+
+
+\backslash
+exerc1[10]
+\end_layout
+
+\end_inset
+
+If 
+\begin_inset Formula $\alpha$
+\end_inset
+
+ and 
+\begin_inset Formula $\beta$
+\end_inset
+
+ are real numbers with 
+\begin_inset Formula $\alpha<\beta$
+\end_inset
+
+,
+ how would you generate a random real number uniformly distributed between 
+\begin_inset Formula $\alpha$
+\end_inset
+
+ and 
+\begin_inset Formula $\beta$
+\end_inset
+
+?
+\end_layout
+
+\begin_layout Standard
+\begin_inset ERT
+status open
+
+\begin_layout Plain Layout
+
+
+\backslash
+answer 
+\end_layout
+
+\end_inset
+
+By taking a uniformly distributed real number 
+\begin_inset Formula $U$
+\end_inset
+
+ between 0 and 1 and returning 
+\begin_inset Formula $\alpha+(\beta-\alpha)U$
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Standard
+\begin_inset ERT
+status open
+
+\begin_layout Plain Layout
+
+
+\backslash
+rexerc3[14]
+\end_layout
+
+\end_inset
+
+Discuss treating 
+\begin_inset Formula $U$
+\end_inset
+
+ as an integer and computing its 
+\emph on
+remainder
+\emph default
+ 
+\begin_inset Formula $\bmod k$
+\end_inset
+
+ to get a random integer between 0 and 
+\begin_inset Formula $k-1$
+\end_inset
+
+,
+ instead of multiplying as suggested in the text.
+ Thus (1) would be changed to
+\begin_inset Formula 
+\begin{align*}
+ & \mathtt{ENTA\ 0}; &  & \mathtt{LDX\ U}; &  & \mathtt{DIV\ K},
+\end{align*}
+
+\end_inset
+
+with the result appearing in register X.
+ Is this a good method?
+\end_layout
+
+\begin_layout Standard
+\begin_inset ERT
+status open
+
+\begin_layout Plain Layout
+
+
+\backslash
+answer 
+\end_layout
+
+\end_inset
+
+While in theory it's as good as the other one,
+ less significant bits tend to be less random in many random number generation methods,
+ so we prefer the method in the text.
+ This is specially important if we use a linear congruential method where 
+\begin_inset Formula $k$
+\end_inset
+
+ is not relatively prime to 
+\begin_inset Formula $m$
+\end_inset
+
+.
+ In addition,
+ if 
+\begin_inset Formula $k$
+\end_inset
+
+ is not much smaller than 
+\begin_inset Formula $m$
+\end_inset
+
+,
+ numbers lower than 
+\begin_inset Formula $m\bmod k$
+\end_inset
+
+ will have a larger probability of appearing,
+ whereas in the method of the text a similar effect is observed but the numbers with higher probability are more evenly distributed.
+\end_layout
+
+\begin_layout Standard
+\begin_inset ERT
+status open
+
+\begin_layout Plain Layout
+
+
+\backslash
+rexerc5[21]
+\end_layout
+
+\end_inset
+
+Suggest an efficient way to compute a random variable with the distribution 
+\begin_inset Formula $F(x)=px+qx^{2}+rx^{3}$
+\end_inset
+
+,
+ where 
+\begin_inset Formula $p\geq0$
+\end_inset
+
+,
+ 
+\begin_inset Formula $q\geq0$
+\end_inset
+
+,
+ 
+\begin_inset Formula $r\geq0$
+\end_inset
+
+,
+ and 
+\begin_inset Formula $p+q+r=1$
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Standard
+\begin_inset ERT
+status open
+
+\begin_layout Plain Layout
+
+
+\backslash
+answer 
+\end_layout
+
+\end_inset
+
+Assume this formula for the distribution is valid for 
+\begin_inset Formula $x\in[0,1]$
+\end_inset
+
+,
+ as 
+\begin_inset Formula $F(0)=0$
+\end_inset
+
+ and 
+\begin_inset Formula $F(1)=1$
+\end_inset
+
+.
+ Now,
+ 
+\begin_inset Formula $F(x)=x(rx^{2}+qx+p)$
+\end_inset
+
+,
+ or more precisely,
+\begin_inset Formula 
+\[
+F(x)=\max\{0,\min\{1,x\}\}\cdot\begin{cases}
+0, & x\leq-\tfrac{q}{2r};\\
+\min\{1,rx^{2}+qx+p\}, & x\geq-\tfrac{q}{2r},
+\end{cases}
+\]
+
+\end_inset
+
+so we may take 
+\begin_inset Formula $U_{1}$
+\end_inset
+
+ and 
+\begin_inset Formula $U_{2}$
+\end_inset
+
+ uniformly distributed between 0 and 1 and return
+\begin_inset Formula 
+\[
+\max\left\{ U_{1},-\frac{q}{2r}+\sqrt{\left(\frac{q}{2r}\right)^{2}-\frac{p}{r}+\frac{U_{2}}{r}}\right\} .
+\]
+
+\end_inset
+
+Note that we take the 
+\begin_inset Quotes eld
+\end_inset
+
+
+\begin_inset Formula $+$
+\end_inset
+
+
+\begin_inset Quotes erd
+\end_inset
+
+ sign in the solution to the quadratic equation to choose the solution in the right side of the parabola.
+ The method in the book is more efficient.
+\end_layout
+
+\begin_layout Standard
+\begin_inset ERT
+status open
+
+\begin_layout Plain Layout
+
+
+\backslash
+rexerc7[20]
+\end_layout
+
+\end_inset
+
+(A.
+ J.
+ Walker) Suppose we have a bunch of cubes of 
+\begin_inset Formula $k$
+\end_inset
+
+ different colors,
+ say 
+\begin_inset Formula $n_{j}$
+\end_inset
+
+ cubes of color 
+\begin_inset Formula $C_{j}$
+\end_inset
+
+ for 
+\begin_inset Formula $1\leq j\leq k$
+\end_inset
+
+,
+ and we also have 
+\begin_inset Formula $k$
+\end_inset
+
+ boxes 
+\begin_inset Formula $\{B_{1},\dots,B_{k}\}$
+\end_inset
+
+ each of which can hold exactly 
+\begin_inset Formula $n$
+\end_inset
+
+ cubes.
+ Furthermore 
+\begin_inset Formula $n_{1}+\dots+n_{k}=kn$
+\end_inset
+
+,
+ so the cubes will just fit into the boxes.
+ Prove (constructively) that there is always a way to put the cubes into the boxes so that each box contains at most two different colors of cubes;
+ in fact,
+ there is a way to do it so that,
+ wherever box 
+\begin_inset Formula $B_{j}$
+\end_inset
+
+ contains two colors,
+ one of those colors is 
+\begin_inset Formula $C_{j}$
+\end_inset
+
+.
+ Show how to use this principle to compute the 
+\begin_inset Formula $P$
+\end_inset
+
+ and 
+\begin_inset Formula $Y$
+\end_inset
+
+ tables required in (3),
+ given a probability distribution 
+\begin_inset Formula $(p_{1},\dots,p_{k})$
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Standard
+\begin_inset ERT
+status open
+
+\begin_layout Plain Layout
+
+
+\backslash
+answer 
+\end_layout
+
+\end_inset
+
+We prove this by induction on 
+\begin_inset Formula $k$
+\end_inset
+
+.
+ For 
+\begin_inset Formula $k=1$
+\end_inset
+
+ this is trivial.
+ For 
+\begin_inset Formula $k>1$
+\end_inset
+
+,
+ if there is a 
+\begin_inset Formula $j$
+\end_inset
+
+ such that 
+\begin_inset Formula $n=n_{j}$
+\end_inset
+
+,
+ we pack the 
+\begin_inset Formula $n_{j}$
+\end_inset
+
+ cubes of color 
+\begin_inset Formula $C_{j}$
+\end_inset
+
+ into 
+\begin_inset Formula $B_{j}$
+\end_inset
+
+;
+ otherwise there exists 
+\begin_inset Formula $i$
+\end_inset
+
+ and 
+\begin_inset Formula $j$
+\end_inset
+
+ such that 
+\begin_inset Formula $n_{i}<n<n_{j}$
+\end_inset
+
+,
+ so we pack the 
+\begin_inset Formula $n_{i}$
+\end_inset
+
+ cubes of color 
+\begin_inset Formula $C_{i}$
+\end_inset
+
+ and also 
+\begin_inset Formula $n-n_{i}$
+\end_inset
+
+ cubes of color 
+\begin_inset Formula $C_{j}$
+\end_inset
+
+ into 
+\begin_inset Formula $B_{i}$
+\end_inset
+
+.
+ Either way we are left with 
+\begin_inset Formula $k-1$
+\end_inset
+
+ colors and 
+\begin_inset Formula $k-1$
+\end_inset
+
+ boxes.
+\end_layout
+
+\begin_layout Standard
+To construct the tables,
+ we note that this proof doesn't require 
+\begin_inset Formula $n$
+\end_inset
+
+ and the 
+\begin_inset Formula $n_{j}$
+\end_inset
+
+s to be integers,
+ they can be arbitrary nonnegative nonnegative reals as long as the conditions are met,
+ so we can let 
+\begin_inset Formula $n=\frac{1}{k}$
+\end_inset
+
+ and 
+\begin_inset Formula $n_{j}=p_{j}$
+\end_inset
+
+ (the probability of event 
+\begin_inset Formula $x_{j}$
+\end_inset
+
+) for each 
+\begin_inset Formula $j$
+\end_inset
+
+.
+ Then,
+ if after proceeding as in this proof,
+ we find that 
+\begin_inset Formula $B_{i}$
+\end_inset
+
+ has 
+\begin_inset Formula $p$
+\end_inset
+
+ cubes of color 
+\begin_inset Formula $C_{i}$
+\end_inset
+
+ and 
+\begin_inset Formula $n-p$
+\end_inset
+
+ cubes of some other color 
+\begin_inset Formula $C_{j}$
+\end_inset
+
+ (where 
+\begin_inset Formula $p\in[0,n]$
+\end_inset
+
+),
+ we set 
+\begin_inset Formula $P=pk$
+\end_inset
+
+ and 
+\begin_inset Formula $Y=j$
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Standard
+\begin_inset ERT
+status open
+
+\begin_layout Plain Layout
+
+
+\backslash
+rexerc10[HM24]
+\end_layout
+
+\end_inset
+
+Explain how to calculate auxiliary constants 
+\begin_inset Formula $P_{j}$
+\end_inset
+
+,
+ 
+\begin_inset Formula $Q_{j}$
+\end_inset
+
+,
+ 
+\begin_inset Formula $Y_{j}$
+\end_inset
+
+,
+ 
+\begin_inset Formula $Z_{j}$
+\end_inset
+
+,
+ 
+\begin_inset Formula $S_{j}$
+\end_inset
+
+,
+ 
+\begin_inset Formula $D_{j}$
+\end_inset
+
+,
+ 
+\begin_inset Formula $E_{j}$
+\end_inset
+
+ so that Algorithm M delivers answers with the correct distribution.
+\end_layout
+
+\begin_layout Standard
+\begin_inset ERT
+status open
+
+\begin_layout Plain Layout
+
+
+\backslash
+answer 
+\end_layout
+
+\end_inset
+
+Let 
+\begin_inset Formula $p_{1},\dots,p_{31}$
+\end_inset
+
+ be the probabilities as in the text and 
+\begin_inset Formula $p_{0}=0$
+\end_inset
+
+.
+ We may calculate probabilities 
+\begin_inset Formula $p_{16},\dots,p_{30}$
+\end_inset
+
+ by a numeral method like Simpson's rule,
+ and this gives us values that are all lower than 
+\begin_inset Formula $\frac{1}{32}$
+\end_inset
+
+.
+ 
+\end_layout
+
+\begin_layout Standard
+First,
+ calculate 
+\begin_inset Formula $p_{1},\dots,p_{31}$
+\end_inset
+
+ as defined in the text,
+ obtaining that 
+\begin_inset Formula $p_{16},\dots,p_{31}<\frac{1}{32}$
+\end_inset
+
+ (we may use some numerical method like the Simpson's rule),
+ and set 
+\begin_inset Formula $p_{0}=0$
+\end_inset
+
+.
+ Then,
+ operate like in Exercise 7 for these probabilities to obtain tables 
+\begin_inset Formula $(P_{j},Y'_{j})_{j=0}^{31}$
+\end_inset
+
+,
+ ensuring that 
+\begin_inset Formula $p_{15},\dots,p_{31}$
+\end_inset
+
+ are considered first so that all the 
+\begin_inset Formula $Y'_{j}\in\{1,\dots,15\}$
+\end_inset
+
+.
+ This gives us the 
+\begin_inset Formula $P_{j}$
+\end_inset
+
+,
+ and then we set 
+\begin_inset Formula $Y_{j}=\frac{Y'_{j}-1}{5}$
+\end_inset
+
+ and 
+\begin_inset Formula $Z_{j}=\frac{1}{5(1-P_{j})}$
+\end_inset
+
+ for 
+\begin_inset Formula $0\leq j\leq31$
+\end_inset
+
+.
+ We also set 
+\begin_inset Formula $S_{j}=\frac{j-1}{5}$
+\end_inset
+
+ for 
+\begin_inset Formula $1\leq j\leq16$
+\end_inset
+
+ and 
+\begin_inset Formula $Q_{j}=\frac{1}{5P_{j}}$
+\end_inset
+
+ for 
+\begin_inset Formula $1\leq j\leq15$
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Standard
+Now,
+ the wedge functions are
+\begin_inset Formula 
+\[
+f_{j+15}(x)\coloneqq\sqrt{\frac{2}{\pi}}\left(\text{e}^{-x^{2}/2}-\text{e}^{-j^{2}/50}\right)\bigg/p_{j+15}
+\]
+
+\end_inset
+
+for 
+\begin_inset Formula $1\leq j\leq15$
+\end_inset
+
+.
+ For 
+\begin_inset Formula $1\leq j\leq5$
+\end_inset
+
+,
+ the curve is concave,
+ so we set
+\begin_inset Formula 
+\begin{align*}
+b_{j} & =-f'_{j+15}(S_{j+1})=\frac{j}{5p_{j+15}}\sqrt{\frac{2}{\pi}}\text{e}^{-j^{2}/50}, & a_{j} & =f_{j+15}(S_{j}).
+\end{align*}
+
+\end_inset
+
+For 
+\begin_inset Formula $6\leq j\leq15$
+\end_inset
+
+ the curve is convex and we set 
+\begin_inset Formula $b_{j}=5f_{j+15}(S_{j})$
+\end_inset
+
+ (the slope) and 
+\begin_inset Formula $a_{j}=f_{j+15}(x)+(x-S_{j})b$
+\end_inset
+
+,
+ where 
+\begin_inset Formula $x\in[S_{j},S_{j+1}]$
+\end_inset
+
+ is such that 
+\begin_inset Formula $f'_{j+15}(x)=-b_{j}$
+\end_inset
+
+.
+ Now,
+\begin_inset Formula 
+\begin{multline*}
+f'_{j+15}(x)=-\sqrt{\frac{2}{\pi}}\frac{x\text{e}^{-x^{2}/2}}{p_{j+15}}=-\frac{5}{p_{j+15}}\sqrt{\frac{2}{\pi}}(\text{e}^{-(j-1)^{2}/50}-\text{e}^{-j^{2}/50})=-b_{j}\iff\\
+x\text{e}^{-x^{2}/2}=5(\text{e}^{-(j-1)^{2}/50}-\text{e}^{-j^{2}/50}).
+\end{multline*}
+
+\end_inset
+
+We know this value exists by Lagrange's mean value theorem,
+ although we need to calculate it numerically (for example,
+ by Newton's method,
+ using that 
+\begin_inset Formula $\frac{\text{d}}{\text{d}x}(x\text{e}^{-x^{2}/2})=(1-x^{2})\text{e}^{-x^{2}/2}$
+\end_inset
+
+).
+ Then,
+ for for 
+\begin_inset Formula $1\leq j\leq15$
+\end_inset
+
+,
+ we set 
+\begin_inset Formula $D_{j+15}=a_{j}/b_{j}$
+\end_inset
+
+ and 
+\begin_inset Formula $E_{j+15}=\sqrt{\frac{2}{\pi}}\frac{e^{-j^{2}/50}}{b_{j}p_{j+15}}$
+\end_inset
+
+,
+ so
+\begin_inset Formula 
+\[
+E_{j+15}=\begin{cases}
+\frac{5}{j}, & 1\leq j\leq5;\\
+\frac{5}{\text{e}^{(2j-1)/50}-1}, & 6\leq j\leq15.
+\end{cases}
+\]
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+\begin_inset Note Note
+status open
+
+\begin_layout Plain Layout
+\begin_inset Note Comment
+status open
+
+\begin_layout Plain Layout
+\begin_inset listings
+inline false
+status open
+
+\begin_layout Plain Layout
+
+use POSIX;
+\end_layout
+
+\begin_layout Plain Layout
+
+\end_layout
+
+\begin_layout Plain Layout
+
+my (@D,
+ @E,
+ @P,
+ @Q,
+ @S,
+ @Y,
+ @Z);
+\end_layout
+
+\begin_layout Plain Layout
+
+\end_layout
+
+\begin_layout Plain Layout
+
+sub normal {
+\end_layout
+
+\begin_layout Plain Layout
+
+    my $rand  = rand;
+\end_layout
+
+\begin_layout Plain Layout
+
+    my $sign  = POSIX::floor( $rand * 2 );
+\end_layout
+
+\begin_layout Plain Layout
+
+    my $piece = POSIX::floor( $rand * 64 ) % 32;
+\end_layout
+
+\begin_layout Plain Layout
+
+    my $dev   = ( $rand * 64 ) % 1;
+\end_layout
+
+\begin_layout Plain Layout
+
+    my $abs;
+\end_layout
+
+\begin_layout Plain Layout
+
+    if ( $dev >= $P[$piece] ) {
+\end_layout
+
+\begin_layout Plain Layout
+
+        $abs = $Y[$piece] + $dev * $Z[$piece];
+\end_layout
+
+\begin_layout Plain Layout
+
+    }
+\end_layout
+
+\begin_layout Plain Layout
+
+    elsif ( $piece <= 15 ) {
+\end_layout
+
+\begin_layout Plain Layout
+
+        $abs = $S[$piece] + $dev * $Q[$piece];
+\end_layout
+
+\begin_layout Plain Layout
+
+    }
+\end_layout
+
+\begin_layout Plain Layout
+
+    elsif ( $piece < 31 ) {
+\end_layout
+
+\begin_layout Plain Layout
+
+        $piece -= 16;
+\end_layout
+
+\begin_layout Plain Layout
+
+        my ( $u,
+ $v );
+\end_layout
+
+\begin_layout Plain Layout
+
+        do {
+\end_layout
+
+\begin_layout Plain Layout
+
+            ( $u,
+ $v ) = ( rand,
+ rand );
+\end_layout
+
+\begin_layout Plain Layout
+
+            ( $u,
+ $v ) = ( $v,
+ $u ) if $u > $v;
+\end_layout
+
+\begin_layout Plain Layout
+
+            $abs = $S[$piece] + $u / 5;
+\end_layout
+
+\begin_layout Plain Layout
+
+        } while ( $v < $D[$piece]
+\end_layout
+
+\begin_layout Plain Layout
+
+            || $v <= $u +
+\end_layout
+
+\begin_layout Plain Layout
+
+            $E[$piece] * ( exp( ( $S[ $piece + 1 ]**2 - $abs**2 ) / 2 ) - 1 ) );
+\end_layout
+
+\begin_layout Plain Layout
+
+    }
+\end_layout
+
+\begin_layout Plain Layout
+
+    else {
+\end_layout
+
+\begin_layout Plain Layout
+
+        do { $abs = sqrt( 9 - 2 * log rand ) } while ( $abs * rand ) >= 3;
+\end_layout
+
+\begin_layout Plain Layout
+
+    }
+\end_layout
+
+\begin_layout Plain Layout
+
+    return $sign ?
+ -$abs :
+ $abs;
+\end_layout
+
+\begin_layout Plain Layout
+
+}
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Plain Layout
+\begin_inset listings
+inline false
+status open
+
+\begin_layout Plain Layout
+
+local D = ...
+\end_layout
+
+\begin_layout Plain Layout
+
+local E = ...
+\end_layout
+
+\begin_layout Plain Layout
+
+local P = ...
+\end_layout
+
+\begin_layout Plain Layout
+
+local Q = ...
+\end_layout
+
+\begin_layout Plain Layout
+
+local S = ...
+\end_layout
+
+\begin_layout Plain Layout
+
+local Y = ...
+\end_layout
+
+\begin_layout Plain Layout
+
+local Z = ...
+\end_layout
+
+\begin_layout Plain Layout
+
+\end_layout
+
+\begin_layout Plain Layout
+
+local function normal ()
+\end_layout
+
+\begin_layout Plain Layout
+
+    local rand = math.random()
+\end_layout
+
+\begin_layout Plain Layout
+
+    local sign = math.floor(rand * 2)
+\end_layout
+
+\begin_layout Plain Layout
+
+    local piece = math.floor(rand * 64) % 32
+\end_layout
+
+\begin_layout Plain Layout
+
+    local dev = (rand * 64) % 1
+\end_layout
+
+\begin_layout Plain Layout
+
+    local abs
+\end_layout
+
+\begin_layout Plain Layout
+
+    if dev >= P[piece] then
+\end_layout
+
+\begin_layout Plain Layout
+
+        abs = Y[piece] + dev * Z[piece]
+\end_layout
+
+\begin_layout Plain Layout
+
+    elseif piece <= 15 then
+\end_layout
+
+\begin_layout Plain Layout
+
+        abs = S[piece] + dev * Q[piece]
+\end_layout
+
+\begin_layout Plain Layout
+
+    elseif piece < 31 then
+\end_layout
+
+\begin_layout Plain Layout
+
+        local u,
+ v
+\end_layout
+
+\begin_layout Plain Layout
+
+        repeat
+\end_layout
+
+\begin_layout Plain Layout
+
+            u,
+ v = math.random(),
+ math.random()
+\end_layout
+
+\begin_layout Plain Layout
+
+            if u > v then u,
+ v = v,
+ u end
+\end_layout
+
+\begin_layout Plain Layout
+
+            abs = S[piece-15] + u/5
+\end_layout
+
+\begin_layout Plain Layout
+
+        until v >= D[piece-15] or
+\end_layout
+
+\begin_layout Plain Layout
+
+            v > u + E[piece-15] * (math.exp((S[piece-14]^2 - abs^2) / 2) - 1)
+\end_layout
+
+\begin_layout Plain Layout
+
+    else
+\end_layout
+
+\begin_layout Plain Layout
+
+        repeat
+\end_layout
+
+\begin_layout Plain Layout
+
+            abs = math.sqrt(9 - 2 * math.log(math.random()))
+\end_layout
+
+\begin_layout Plain Layout
+
+        until abs * rand >= 3
+\end_layout
+
+\begin_layout Plain Layout
+
+    end
+\end_layout
+
+\begin_layout Plain Layout
+
+    if sign == 0 then return abs else return -abs end
+\end_layout
+
+\begin_layout Plain Layout
+
+end
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Plain Layout
+\begin_inset listings
+inline false
+status open
+
+\begin_layout Plain Layout
+
+from math import exp,
+ floor,
+ log,
+ sqrt
+\end_layout
+
+\begin_layout Plain Layout
+
+from random import random
+\end_layout
+
+\begin_layout Plain Layout
+
+\end_layout
+
+\begin_layout Plain Layout
+
+D,
+ E,
+ P,
+ Q,
+ S,
+ Y,
+ Z = [],
+ [],
+ [],
+ [],
+ [],
+ [],
+ []
+\end_layout
+
+\begin_layout Plain Layout
+
+\end_layout
+
+\begin_layout Plain Layout
+
+\end_layout
+
+\begin_layout Plain Layout
+
+def normal():
+\end_layout
+
+\begin_layout Plain Layout
+
+    rand = random()
+\end_layout
+
+\begin_layout Plain Layout
+
+    sign = floor(rand * 2)
+\end_layout
+
+\begin_layout Plain Layout
+
+    piece = floor(rand * 64) % 32
+\end_layout
+
+\begin_layout Plain Layout
+
+    dev = (rand * 64) % 1
+\end_layout
+
+\begin_layout Plain Layout
+
+    if dev >= P[piece]:
+\end_layout
+
+\begin_layout Plain Layout
+
+        result = Y[piece] + dev * Z[piece]
+\end_layout
+
+\begin_layout Plain Layout
+
+    elif piece <= 15:
+\end_layout
+
+\begin_layout Plain Layout
+
+        result = S[piece] + dev * Q[piece]
+\end_layout
+
+\begin_layout Plain Layout
+
+    elif piece < 31:
+\end_layout
+
+\begin_layout Plain Layout
+
+        piece = piece - 16
+\end_layout
+
+\begin_layout Plain Layout
+
+        while True:
+\end_layout
+
+\begin_layout Plain Layout
+
+            u,
+ v = random(),
+ random()
+\end_layout
+
+\begin_layout Plain Layout
+
+            if u > v:
+\end_layout
+
+\begin_layout Plain Layout
+
+                u,
+ v = v,
+ u
+\end_layout
+
+\begin_layout Plain Layout
+
+            result = S[piece] + u/5
+\end_layout
+
+\begin_layout Plain Layout
+
+            if v >= D[piece]:
+\end_layout
+
+\begin_layout Plain Layout
+
+                break
+\end_layout
+
+\begin_layout Plain Layout
+
+            if v < u + E[piece] * (exp((S[piece+1]**2 - abs**2)/2) - 1):
+\end_layout
+
+\begin_layout Plain Layout
+
+                break
+\end_layout
+
+\begin_layout Plain Layout
+
+    else:
+\end_layout
+
+\begin_layout Plain Layout
+
+        while True:
+\end_layout
+
+\begin_layout Plain Layout
+
+            result = sqrt(9 - 2*log(random()))
+\end_layout
+
+\begin_layout Plain Layout
+
+            if result * random() < 3:
+\end_layout
+
+\begin_layout Plain Layout
+
+                break
+\end_layout
+
+\begin_layout Plain Layout
+
+    if sign:
+\end_layout
+
+\begin_layout Plain Layout
+
+        result = -result
+\end_layout
+
+\begin_layout Plain Layout
+
+    return result
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+\begin_inset ERT
+status open
+
+\begin_layout Plain Layout
+
+
+\backslash
+rexerc16[HM22]
+\end_layout
+
+\end_inset
+
+(J.
+ H.
+ Ahrens.) Develop an algorithm for gamma deviates of order 
+\begin_inset Formula $a$
+\end_inset
+
+ when 
+\begin_inset Formula $0<a<1$
+\end_inset
+
+,
+ using the rejection method with 
+\begin_inset Formula $cg(t)=t^{a-1}/\Gamma(a)$
+\end_inset
+
+ for 
+\begin_inset Formula $0<t<1$
+\end_inset
+
+,
+ and with 
+\begin_inset Formula $cg(t)=\text{e}^{-t}/\Gamma(a)$
+\end_inset
+
+ for 
+\begin_inset Formula $t\geq1$
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Standard
+\begin_inset ERT
+status open
+
+\begin_layout Plain Layout
+
+
+\backslash
+answer 
+\end_layout
+
+\end_inset
+
+The density function for the gamma distribution is
+\begin_inset Formula 
+\begin{align*}
+f(x) & \coloneqq\frac{1}{\Gamma(a)}x^{a-1}\text{e}^{-x}, & x & \geq0.
+\end{align*}
+
+\end_inset
+
+Let
+\begin_inset Formula 
+\[
+g(t)\coloneqq\begin{cases}
+t^{a-1}, & 0<t<1;\\
+\text{e}^{-t}, & t\geq1.
+\end{cases}
+\]
+
+\end_inset
+
+Then,
+\begin_inset Formula 
+\begin{align*}
+G_{0}(x)\coloneqq\int_{0}^{x}g(t)\text{d}t & =\begin{cases}
+\int_{0}^{x}t^{a-1}\text{d}t=\left[\frac{t^{a}}{a}\right]_{t=0}^{x}=\frac{x^{a}}{a}, & 0<t<1;\\
+a^{-1}+\int_{1}^{x}\text{e}^{-t}\text{d}t=\frac{1}{a}-\left[\text{e}^{-t}\right]_{t=1}^{x}=a^{-1}+\text{e}^{-1}-\text{e}^{-x}, & t\geq1.
+\end{cases}
+\end{align*}
+
+\end_inset
+
+If 
+\begin_inset Formula $M\coloneqq\lim_{x\to\infty}G(x)=a^{-1}+\text{e}^{-1}$
+\end_inset
+
+,
+ the real probability distribution function is 
+\begin_inset Formula $G(x)\coloneqq G_{0}(x)/M$
+\end_inset
+
+.
+ Thus,
+ a possible method could be the following:
+ Generate to independent deviates 
+\begin_inset Formula $U$
+\end_inset
+
+ and 
+\begin_inset Formula $V$
+\end_inset
+
+ uniformly distributed between 0 and 1.
+ If 
+\begin_inset Formula $U<a^{-1}/M$
+\end_inset
+
+,
+ set 
+\begin_inset Formula $X\coloneqq\sqrt[a]{aMU}=\sqrt[a]{1+\frac{a}{\text{e}}}\sqrt[a]{U}$
+\end_inset
+
+ and 
+\begin_inset Formula $Y\coloneqq\Gamma(a)f(X)/g(X)=\text{e}^{-X}$
+\end_inset
+
+.
+ Otherwise set 
+\begin_inset Formula $X\coloneqq-\ln M\cdot\ln(1-U)$
+\end_inset
+
+ and 
+\begin_inset Formula $Y\coloneqq\Gamma(a)f(X)/g(X)=X^{a-1}$
+\end_inset
+
+.
+ If 
+\begin_inset Formula $V\geq Y$
+\end_inset
+
+,
+ reject 
+\begin_inset Formula $X$
+\end_inset
+
+ and repeat all the steps.
+ Otherwise return 
+\begin_inset Formula $X$
+\end_inset
+
+.
+\end_layout
+
+\begin_layout Standard
+\begin_inset ERT
+status open
+
+\begin_layout Plain Layout
+
+
+\backslash
+rexerc17[M24]
+\end_layout
+
+\end_inset
+
+What is the 
+\emph on
+distribution function
+\emph default
+ 
+\begin_inset Formula $F(x)$
+\end_inset
+
+ for the geometric distribution with probability 
+\begin_inset Formula $p$
+\end_inset
+
+?
+ What is the 
+\emph on
+generating function
+\emph default
+ 
+\begin_inset Formula $G(z)$
+\end_inset
+
+?
+ What are the mean and standard deviation of this distribution?
+\end_layout
+
+\begin_layout Standard
+\begin_inset ERT
+status open
+
+\begin_layout Plain Layout
+
+
+\backslash
+answer 
+\end_layout
+
+\end_inset
+
+
+\begin_inset Formula 
+\[
+F(x)=P(X\leq x)=1-P(X>x)=1-P(X>\lfloor x\rfloor)=1-(1-p)^{\max\{\lfloor x\rfloor,0\}}.
+\]
+
+\end_inset
+
+The generating function is
+\begin_inset Formula 
+\[
+G(z)\coloneqq\sum_{k\geq1}(1-p)^{k-1}pz^{k}=zp\sum_{k\geq0}(1-p)^{k}z^{k}=\frac{zp}{1-(1-p)z}=\frac{zp}{1-qz},
+\]
+
+\end_inset
+
+where 
+\begin_inset Formula $q\coloneqq1-p$
+\end_inset
+
+.
+ Then
+\begin_inset Formula 
+\begin{align*}
+G'(z) & =\frac{p-pqz+pqz}{(1-qz)^{2}}=\frac{p}{(1-qz)^{2}}, & G''(z) & =-\frac{2pq}{(1-qz)^{3}},
+\end{align*}
+
+\end_inset
+
+so the mean is
+\begin_inset Formula 
+\[
+G'(1)=\frac{p}{(1-q)^{2}}=\frac{p}{p^{2}}=\frac{1}{p},
+\]
+
+\end_inset
+
+and the variance is
+\begin_inset Formula 
+\[
+G''(1)+G'(1)-G'(1)^{2}=-\frac{2pq}{(1-q)^{3}}+\frac{1}{p}-\frac{1}{p^{2}}=\frac{2q}{p^{2}}+\frac{1}{p}-\frac{1}{p^{2}}=\frac{2q+p-1}{p^{2}}=\frac{q}{p^{2}},
+\]
+
+\end_inset
+
+so the standard deviation is 
+\begin_inset Formula $\frac{\sqrt{q}}{p}$
+\end_inset
+
+.
+\end_layout
+
+\end_body
+\end_document
-- 
cgit v1.2.3