Problem 1

Given an LCG with parameters \(a,c,m\), prove that

which shows that the \((n+k)th\) term can be computed directly from the \(nth\) term.

Solution:

We know from D. H. Lehmer’s linear congruential generator that
\[\begin{equation} x_n \equiv ax_{n-1}+c \mod m \end{equation}\]
where \(a\) is called the multiplier, \(c\) is called the shift or increment, and \(m\) is called the modulus of the generator. The given equation is also an LCG. We can prove this by induction method. Since \(k\ge 0\) so, let \(k=0\). Then the given relation can be written as

If \(k=1\). Then the given relation can be written as
\[\begin{align*} x_{n+1}& \equiv ax_n+\frac{a-1}{a-1}c \mod m\\ &\equiv ax_n+c \mod m \end{align*}\]

If \(k=2\). Then the given relation can be written as
\[\begin{align*} x_{n+2}& \equiv a^2x_n+\frac{a^2-1}{a-1}c \mod m\\ &\equiv a^2x_n+(a+1)c \mod m\\ &\equiv a^2x_n+ac+c \mod m \\ &\equiv a(ax_n+c)+c \mod m\\ &\equiv ax_{n+1}+c \mod m \end{align*}\]

Now for any \(k=p\) where \(p\in \mathbb{N}\), \[\begin{align*} x_{n+p}& \equiv a^px_n+\frac{a^p-1}{a-1}c \mod m \\ \end{align*}\]

Now by the method of induction, the given equation would be a lcg if it holds for any \(k=p\in \mathbb{N}\) then it must hold for \(k=p+1\) where \(p\in \mathbb{N}\). Now from equation (1) \[\begin{align*} x_{n+p+1} &\equiv ax_{(n+p+1)-1}+c \mod m\\ & \equiv ax_{n+p}+c \mod m \\ & \equiv a(a^px_n+\frac{a^p-1}{a-1}c) +c \mod m\\ & \equiv a^{p+1}x_n+(a\frac{a^p-1}{a-1}+1)c \mod m\\ & \equiv a^{p+1}x_n+\frac{a^{p+1}-1}{a-1}c \mod m\\ \end{align*}\]

Which proves that \(x_{n+k}=a^kx_n+\frac{(a^k-1)}{a-1}c (\mod m)\); \((a\ge 2, k\ge0)\) is an lcg such that \((n+k)th\) term can be computed directly from the \(nth\) term.

Problem 2

(a)

If \(U\) and \(V\) are independently distributed random variables from the uniform distribution \(U(0,1)\) show that \(U+V (\mod 1)\) is also \(U(0,1)\).

Solution

Let \(Z=U+V\) where \(U\) and \(V\) are independently distributed random variables from the uniform distribution \(U(0,1)\). So the minimum value that \(Z\) can have is \(0\) and the maximum value could be \(2\). If \(f_U(u)\) is the PDF of \(U\) and \(f_V(v)\) is the PDF of \(V\) then the PDF of \(Z\) can be found from the convolution of two distribution as follows \[\begin{align*} f_Z(z)=\int_{-\infty}^{+\infty}f_U(u)f_V(z-u)du=\begin{cases} z & \text{for} \hspace{2mm} 0 < z < 1\\ 2-z & \text{for} \hspace{2mm} 1 \le z <2\\ 0 & \text{otherwise} \end{cases} \end{align*}\] Now for any \(x\in (0,1)\) \[\begin{align*} \mathbb{P}(U+V (\mod 1) \le x) &= \mathbb{P}(Z \le x)+ \mathbb{P}(1\le Z \le x+1)\\ &= \int_{0}^{x} z dz +\int_{1}^{1+x}(2-z)dz\\ &=x \end{align*}\]

which is the CDF of a random variable distributed \(U(0,1)\)

(b)

A random number generator is designed by

where \(X_0=0, Y_0=1, X_{n+1}=(9X_n+3) \mod 8\) and \(Y_{n+1}=3Y_n \mod 7\) for \(n=0,1,2,\hdots\). Calculate \(R_0,R_1,R_2, \hdots , R_5.\). What is the period of the generator \(\{R_n\}\)

Solution

rand.gen<-function(n){
  RN<-vector(length = n)
  x<-rep(n)
  y<-rep(n)
  x[1]<-0;
  y[1]<-1;
  RN[1]<-(x[1]/8+y[1]/7)%% 1
  for (i in 1:n) {
    x[i+1]<-(9*x[i]+3)%% 8
    y[i+1]<-(3*y[i]) %% 7
    RN[i+1]<-(x[i+1]/8+y[i+1]/7)%% 1
  }
  return(data.frame(X_values=x,Y_values=y,R_values=RN))
}
rand.gen(4)  
##   X_values Y_values   R_values
## 1        0        1 0.14285714
## 2        3        3 0.80357143
## 3        6        2 0.03571429
## 4        1        6 0.98214286
## 5        4        4 0.07142857

So the unique values are

##      R_values
## 1  0.14285714
## 2  0.80357143
## 3  0.03571429
## 4  0.98214286
## 5  0.07142857
## 6  0.58928571
## 7  0.39285714
## 8  0.05357143
## 9  0.28571429
## 10 0.23214286
## 11 0.32142857
## 12 0.83928571
## 13 0.64285714
## 14 0.30357143
## 15 0.53571429
## 16 0.48214286
## 17 0.57142857
## 18 0.08928571
## 19 0.89285714
## 20 0.55357143
## 21 0.78571429
## 22 0.73214286
## 23 0.82142857
## 24 0.33928571

So from the above data we can see that the period is \(24\).

Problem 3

Write a code that would implement RANDU. For debugging purpose print \(x_{1000}\) when the seed is \(x_0=1\)

(a)

Using RANDU generate \(u_1,u_2,\hdots, u_{20,002}\) where \(u=\frac{x_n}{M}\). For all triplets in your sequence, \(u_i, u_{i+1}, u_{i+2}\), in which \(0.5\le u_{i+1} \le 0.51\) plot \(u_i\) versus \(u_{i+2}\). Comment on the pattern of your scatterplot.

options(repos=c(CRAN="<something sensible near you>"))
install.packages("dplR")
## Installing package into 'C:/Users/Rakibul/Documents/R/win-library/3.6'
## (as 'lib' is unspecified)
## Warning: unable to access index for repository <something sensible near you>/src/contrib:
##   scheme not supported in URL '<something sensible near you>/src/contrib/PACKAGES'
## Warning: package 'dplR' is not available (for R version 3.6.1)
## Warning: unable to access index for repository <something sensible near you>/bin/windows/contrib/3.6:
##   scheme not supported in URL '<something sensible near you>/bin/windows/contrib/3.6/PACKAGES'
library(dplyr)
## Warning: package 'dplyr' was built under R version 3.6.2
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
N=20002
A = matrix(0, ncol=3, nrow=N)
seed <- as.double(1)

RANDU <- function() {
  seed <<- ((2^16 + 3) * seed) %% (2^31)
  round(seed/(2^31),6)
}

for (i in 1:N) {
  A[i,]<-c(RANDU(),RANDU(),RANDU())
}
B=as.data.frame(A)
C<-B %>% filter(V2>=0.5, V2<=0.51)

plot(C$V1,C$V3,xlab = "u_i",ylab = "u_(i+3)")

(b)

Generate a sequence of lenght 1002. Use a program that plots points in 3 dimensions and rotates the axes to rotate the points until you can see the 15 planes.

list.of.packages <- c("rgl","ggplot2","knitr","rglwidget")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
lapply(list.of.packages,function(x){library(x,character.only=TRUE)})
## Warning: package 'rgl' was built under R version 3.6.2
## Warning: package 'rglwidget' was built under R version 3.6.2
## The functions in the rglwidget package have been moved to rgl.
## [[1]]
## [1] "rgl"       "dplyr"     "stats"     "graphics"  "grDevices" "utils"    
## [7] "datasets"  "methods"   "base"     
## 
## [[2]]
##  [1] "ggplot2"   "rgl"       "dplyr"     "stats"     "graphics" 
##  [6] "grDevices" "utils"     "datasets"  "methods"   "base"     
## 
## [[3]]
##  [1] "knitr"     "ggplot2"   "rgl"       "dplyr"     "stats"    
##  [6] "graphics"  "grDevices" "utils"     "datasets"  "methods"  
## [11] "base"     
## 
## [[4]]
##  [1] "rglwidget" "knitr"     "ggplot2"   "rgl"       "dplyr"    
##  [6] "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [11] "methods"   "base"
N=1002
A = matrix(0, ncol=3, nrow=N)
seed <- as.double(1)

RANDU <- function() {
  seed <<- ((2^16 + 3) * seed) %% (2^31)
  round(seed/(2^31),6)
}

for (i in 1:N) {
  A[i,]<-c(RANDU(),RANDU(),RANDU())
}
B=as.data.frame(A)
plot3d(B$V1, B$V2, B$V3, type="s", size=1, lit=TRUE,col = rainbow(1000))
A snapshot of the 3D plot

A snapshot of the 3D plot