Exercise 3.2 (p. 94)

The standard Laplace distribution has density \(\mathbf{f(x) = \frac12e^{-|x|},\,x \in \mathbb{R}}\). Use the inverse transform method to generate a random sample of size 1000 from this distribution. Use one of the methods shown in this chapter to compare the generated sample to the target distribution.

\(f_X(x)=\begin{cases}\dfrac12e^x&,\,x<0\\[12pt]\dfrac12e^{-x}&,\,x\geq0\end{cases}\)

\(F_X(x)=\begin{cases}\dfrac12e^x&,\,x<0\\[12pt]1-\dfrac12e^x&,\,x\geq0\end{cases}\)

Case 1: If \(x<0\)

\[\begin{equation} \begin{split} P(X\leq x) & = \int_{-\infty}^{x} {\frac{1}{2}e^t dt}\\ & = \lim_{c\to -\infty}{\int_{c}^{x} {\frac{1}{2}e^t dt}}\\ & = \lim_{c\to -\infty}{\frac{1}{2}e^x|_{c}^{x}}\\ & = \lim_{c\to -\infty}{(\frac{1}{2}e^x-\frac{1}{2}e^c)}\\ & = \frac{1}{2}e^x \end{split} \end{equation}\] Thus \[\begin{equation} \begin{split} U &=\frac{1}{2}e^x\\ 2U &=e^x\\ ln(2U) &=x \end{split} \end{equation}\] Since \(x<0\), \[\begin{equation} \begin{split} ln(2U)&<0\\ 2U&<1\\ U&<\frac{1}{2} \end{split} \end{equation}\]

Case 2: If \(x\geq0\)

\[\begin{equation} \begin{split} P(X\leq x) & = \int_{-\infty}^{0} {\frac{1}{2}e^t dt}+\int_{0}^{x} {\frac{1}{2}e^{-t} dt}\\ & = \frac{1}{2}e^0+\frac{-1}{2}e^{-x}|_0^x\\ & = \frac{1}{2}-\frac{1}{2}e^{-x}+\frac{1}{2}e^0\\ & = 1-\frac{1}{2}e^{-x} \end{split} \end{equation}\]

Thus \[\begin{equation} \begin{split} U &=1-\frac{1}{2}e^{-x}\\ \frac{1}{2}e^{-x} &=1-U\\ e^{-x} &=2-2U\\ -x&=ln(2-2U)\\ x&=-ln(2-2u)\\ \end{split} \end{equation}\] Since \(x<0\), \[\begin{equation} \begin{split} -ln(2-2U)&\geq0\\ ln(2-2U)&\leq0\\ 2-2U&\leq1\\ U&\geq\frac{1}{2} \end{split} \end{equation}\]

\(F^{-1}_X(u)=\begin{cases}\ln(2u)&,\,u<\dfrac12\\[12pt]-\ln[2(1-u)]&,\,u\geq\dfrac12\end{cases}\) 

n <- 1000
u <- runif(n)
x.vect <- numeric(n) # inititating vector space in R's memory

for (i in 1:length(u)) {
  if (u[i] < .5) {
    x.vect[i] <- log(2*u[i])
  }
  
  else {
    x.vect[i] <- -1*log(2*(1-u[i]))
  }
}

hist(x.vect, freq = F, ylim = c(0, 0.5), xlab = "x",main= expression(f(x)==frac(1,2)*e^{-abs(x)}))

# curve draws a function of x; if add = T it uses x-values on the previous
# plot's x-axis
curve(.5*exp(-abs(x)), col = "red", add = T)

# lines takes a vector of x & y ('density' provides both vectors)
lines(density(x.vect), lty = 2, lwd = 2)

Alternate R-Code

n <- 1000
u <- runif(n) 
x<-seq(0,0,length.out=1000)
for (i in 1:1000) {
 if (u[i]<0.5) {x[i]=log(2*u[i])
} else {
          x[i]=-log(2-2*u[i])}}

hist(x,prob = TRUE, main=expression(f(x)==frac(1,2)*e^{-abs(x)}), ylim=c(0,0.5), xlim=c(-7,7))#Plot random x values

y<-sort(x) #use sample x values to plot theoretical curve
lines(y,((1/2)*exp(-abs(y))))

Exercise 3.3 (p. 94)

The Pareto(a, b) distribution has cdf \(\mathbf{F(x)=1-\left(\frac{b}{x}\right)^a\ ,\quad x\geq b > 0\ ,\quad a>0}\). Derive the probability inverse transformation and use the inverse transform method to simulate a random sample from the Pareto(2, 2) distribution. Graph the density histogram of the sample with the Pareto(2, 2) density superimposed for comparison.

For a Pareto(a, b) distribution,

\(f_X(x)=\dfrac{ab^a}{x^{a+1}}\,,\quad x\geq b > 0\,,\quad a>0\ \) and

\(F^{-1}_X(u)=\dfrac{b}{\sqrt[a]{1-u}}\,,\quad 1>u\geq 0\,,\quad a,\,b>0\,.\)

Then, for Pareto(2, 2,),

\(f_X(x)=\dfrac{8}{x^3}\,,\quad x\geq2\,,\)

\(F_X(x)=1-\dfrac{4}{x^2}\,,\quad x\geq2\ ,\) and

\(F^{-1}_X(u)=\dfrac{2}{\sqrt{1-u}}\,,\quad 1>u\geq 0\,.\) 

n <- 1000

# runif does not usually generate either 0 or 1 so no specific arguments
# need to be given to ensure that there are no u values that equal 1
u <- runif(n)
x.pareto <- 2/sqrt(1-u)

hist(x.pareto, freq = F, xlab = "x", main = "Pareto(2,2)")
curve(8/(x)^3, col = "red", add = T)
lines(density(x.pareto), lty = 2, lwd = 2)

Alternate R Code for 3.3

n <- 1000 
u <- runif(n) 
x<-0
for(i in 1:1000){
x[i]= (2/sqrt(1-u[i]))}
hist(x,prob = TRUE, main="Pareto(2,2)")
y<-sort(x) #use sample x values to plot theoretical curve
lines(y,8/y^3)

Exercise 3.11 (p. 95)

Generate a random sample of size 1000 from a normal location mixture. The components of the mixture have \(\mathbf{\mathcal{N}(0,\, 1)}\) and \(\mathbf{\mathcal{N}(3,\, 1)}\) distributions with mixing probabilities \(\mathbf{p_1}\) and \(\mathbf{p_2}\) = \(\mathbf{1-p_1}\). Graph the histogram of the sample with density superimposed, for \(\mathbf{p_1}\) = 0.75. Repeat with different values for \(\mathbf{p_1}\) and observe whether the empirical distribution of the mixture appears to be bimodal. Make a conjecture about the values of \(\mathbf{p_1}\) that produce bimodal mixtures.

\(p_1 = 0.75\): 

n<-1000
x1<-rnorm(1000,0,1)
x2<-rnorm(1000,3,1)
p1<-.75
p2<-1-p1
u<-runif(n)
k<-as.integer(u>p2)
x<-k*x1+(1-k)*x2
par(mfcol=c(3,3))
hist(x, prob=TRUE, ylim=c(0,.4), main=bquote(X==0.75*X[1]+0.25*X[2]))
lines(density(x))


p1<-.5
p2<-1-p1
u<-runif(n)
k<-as.integer(u>p2)
x<-k*x1+(1-k)*x2
hist(x, prob=TRUE, ylim=c(0,.4), main=bquote(X==0.5*X[1]+0.5*X[2]))
lines(density(x))

p1<-.25
p2<-1-p1
u<-runif(n)
k<-as.integer(u>p2)
x<-k*x1+(1-k)*x2
hist(x, prob=TRUE, ylim=c(0,.4), main=bquote(X==0.25*X[1]+0.75*X[2]))
lines(density(x))

p1<-.9
p2<-1-p1
u<-runif(n)
k<-as.integer(u>p2)
x<-k*x1+(1-k)*x2
hist(x, prob=TRUE, ylim=c(0,.4), main=bquote(X==0.9*X[1]+0.1*X[2]))
lines(density(x))

p1<-.6
p2<-1-p1
u<-runif(n)
k<-as.integer(u>p2)
x<-k*x1+(1-k)*x2
hist(x, prob=TRUE, ylim=c(0,.4), main=bquote(X==0.6*X[1]+0.4*X[2]))
lines(density(x))

p1<-.1
p2<-1-p1
u<-runif(n)
k<-as.integer(u>p2)
x<-k*x1+(1-k)*x2
hist(x, prob=TRUE, ylim=c(0,.4), main=bquote(X==0.1*X[1]+0.9*X[2]))
lines(density(x))

From graphs one might conjecture that mixture is bimodal for p-values between 0.2 <p <0.8. Conditions for which the density function of a mixture of two normal distributions is bimodal are investigated in detail by Eisenberger (1964).

\cite Isidore Eisenberger (1964) Genesis of Bimodal Distributions, Technometrics, 6:4, 357-363, DOI: 10.1080/00401706.1964.10490199

For fixed values of the variances \(\sigma_{1}^2\) and \(\sigma_{2}^{2}\) of the normal distributions if the difference between the means is sufficiently small, the distribution of the mixture will be unimodal, independent of the proportions \(p\) and \((1 – p)\), \(0 < p < 1\). If the difference exceeds a critical value which depends on \(\sigma_{1}^2\) and \(\sigma_{2}^2\) the bimodality property then depends on p. Values of p sufficiently close to zero and one always exist for which the distribution is unimodal.

For Graduate Students

From Acceptance-Rejection Method

Part1: Prove E(N)=c.

N~Geo(p) so E(N)= \(\frac{1}{p}\) where \(p\)=probability of success

\[\begin{equation} \begin{split} P(\text{accept for any iteration}) &=\sum_y P(U<\frac{f_x(y)}{cg_y(y)}|Y=y)P(Y=y)\\ &=\sum_y \frac{f_x(y)}{cg_y(y)}*g_y(y)\\ &=\frac{\sum_y f_x(y)}{c}\\ &=\frac{1}{c} \end{split} \end{equation}\]

(Continuous case:

\[\begin{equation} \begin{split} &\int_D{P(U<\frac{f_x(y)}{cg_y(y)}|Y=y)P(Y=y)dy}\\ &\int_D{(\int_0^{\frac{f_x(y)}{cg_y(y)}}{1 dx})(g_y(y))dy}\\ &\int_D{(\frac{f_x(y)}{cg_y(y)})g_y(y)dy }\\ &\frac{1}{c}\int_D{f_x(y)dy}=\frac{1}{c}*1=\frac{1}{c} \end{split} \end{equation}\]

)

So E(N)= \(\frac{1}{p}=\frac{1}{\frac{1}{c}}=c\)

Part2. Prove \(G_y(y)=F_x(y)\), for accepted y’s.

\[\begin{equation} \begin{split} g_y(y)&=P(Y=y|U<\frac{f_x(y)}{cg_y(y)})\\ &=\frac{P(Y=y \cap U<\frac{f_x(y)}{cg_y(y)}} {P(U<\frac{f_x(y)}{cg_y(y)})}\\ &=\frac{P(U<\frac{f_x(y)}{cg_y(y)}|Y=y)P(Y=y)}{P(U<\frac{f_x(y)}{cg_y(y)})}\\ &=\frac{\frac{f_x(y)}{cg_y(y)}g_y(y)}{\frac{1}{c}}\\ &=\frac{{f_x(y)}{c}}{\frac{1}{c}}=f_x(y) \end{split} \end{equation}\]