#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