# plot density func
curve(4*exp(-4*x),from = 0, to = NULL)

true.mean= 1/4
true.sd= sqrt(1/4^2)
inf.fn = function(u){
  x= -(1/4)*log(1-u)
  x
}
inf.fn(u= 0.8)
## [1] 0.4023595
N= 1000
n= 500
mean.est = numeric(N)
sd.est = numeric(N)
for (i in 1:N) {
  u = runif(n,0,1)
  sim.vl= inf.fn(u)
  mean.est[i] = mean(sim.vl)
  sd.est[i] = sd(sim.vl)
}

hist(sim.vl,prob=T, nclas=30, col = "green")
y= seq(0,1.5,0.01)
lines(y,4*exp(-4*y),col="blue", lwd=2)
curve(4*exp(-4*x),from = 0, to = NULL, add=T)

mean.sim = mean(mean.est)
cbind(true.mean, mean.sim)
##      true.mean  mean.sim
## [1,]      0.25 0.2495928

CDF > p f(x)= p(x<= x) = 0 <= p <= 1

generate 5 5 random values from this distribution using the ITM???????????????

1/ Compute the CDF? x<=1 0.2 x<=2 0.35 x<=3 0.60 x<=4 1

2/ Generate U from U(0, 1).

u1 = 0.5074782 x<=2 u1 x<=3 larger x1=3

u2 =0.3067685 x<=1 u2 x<=2 large x2=2

u3 = 0.4269077 x3=3 [1] 0.6931021 x4 =4 [1] 0.08513597 x5=1

set.seed(10)
u1 = runif(1); u1
## [1] 0.5074782
u2 = runif(1); u2
## [1] 0.3067685
u3 = runif(1); u3
## [1] 0.4269077
u4 = runif(1); u4
## [1] 0.6931021
u5 = runif(1); u5
## [1] 0.08513597
p = c(0.20, 0.15,.25,0.4); x= 1:4
cdf= cumsum(p); sim.x=c()
u= c(); n=5
set.seed(10)
for (i in 1:n) {
  u[i]= runif(1)
  sim.x[i]= min(x[cdf>u[i]])
}
u
## [1] 0.50747820 0.30676851 0.42690767 0.69310208 0.08513597
sim.x
## [1] 3 2 3 4 1
#min(x[cdf>u[i]])finding x values such
#that cdf>u[i],then take the min.
p= c(0.20, 0.15,.25,0.4); x=1:4
cdf= cumsum(p) ; sim.x= c()
u= c() ; n=1000
set.seed(10)
for (i in 1:n) {
  u[i]= runif(1)
  sim.x[i]= min(x[cdf>u[i]])
}
mean(sim.x)
## [1] 2.867
table(sim.x)/length(sim.x) #Gives approximated p(X=x)
## sim.x
##     1     2     3     4 
## 0.196 0.146 0.253 0.405
mean(sim.x==1)#Gives approximated p(X=1)
## [1] 0.196
mean(sim.x==2) #Gives approximated p(X=2)
## [1] 0.146
cbind(mean(sim.x==1),mean(sim.x==2),mean(sim.x==3),mean(sim.x==4))
##       [,1]  [,2]  [,3]  [,4]
## [1,] 0.196 0.146 0.253 0.405
p2 <- c(0.1, 0.3, 0.1, 0.05, 0.25, 0.2);x= 1:6

cdf= cumsum(p2); cdf
## [1] 0.10 0.40 0.50 0.55 0.80 1.00
sim.x= c()
u= c() ; n=1000
set.seed(123)
for (i in 1:n) {
  u[i]= runif(1)
  sim.x[i]= min(x[cdf>u[i]])
}
mean(sim.x)
## [1] 3.631
table(sim.x)/length(sim.x)
## sim.x
##     1     2     3     4     5     6 
## 0.099 0.305 0.103 0.050 0.245 0.198

—- pmf —Discrete unif —equally likely outcomes—-

die discunif(1,6) Note: n=b-a+1 for DiscUinf[a, b].

trunc(2.32); floor(2.23)
## [1] 2
## [1] 2

Discrete Uniform Prob Dist.[0, 9] trunc(10*runif(1))comes from n=b-a+1=9-0+1=10

x<- 0:9
sim.x = c()
n= 5
set.seed(10)
for (i in 1:n) {
  sim.x[i]= trunc(10*runif(1))+0
  
}
sim.x
## [1] 5 3 4 6 0
mean(sim.x)
## [1] 3.6
table(sim.x)/length(sim.x) #Gives approximated p(X=x)
## sim.x
##   0   3   4   5   6 
## 0.2 0.2 0.2 0.2 0.2
x = 0:9
sim.x =c()
n=5
set.seed(10)
for (i in 1:n) {
  sim.x[i]= trunc(10* runif(1))+0
}
sim.x
## [1] 5 3 4 6 0
mean(sim.x)
## [1] 3.6
table(sim.x)/length(sim.x)
## sim.x
##   0   3   4   5   6 
## 0.2 0.2 0.2 0.2 0.2

Consider X~DiscU(7, 15), Use R to generate 5 observations from this distribution.

15-7+1
## [1] 9
n= 5
set.seed(10)
sim.x= c()
for (i in 1:n) {
  sim.x[i]= trunc(9*runif(1))+7
  
}
sim.x
## [1] 11  9 10 13  7
mean(sim.x)
## [1] 10
table(sim.x)/length(sim.x)
## sim.x
##   7   9  10  11  13 
## 0.2 0.2 0.2 0.2 0.2

——— Generating from the Geometric distribution——–

first success

  1. Let the random variable X have a Geometric distribution with p= 0.2. Generate 10 values from this distribution.
g.sim= c()
n=10
set.seed(10)
p =0.2
q = 1-p

for (i in 1:n) {
  g.sim[i]= trunc((log(runif(1))/log(q)))+1
}
g.sim
##  [1]  4  6  4  2 12  7  6  6  3  4
table(g.sim)/length(g.sim)
## g.sim
##   2   3   4   6   7  12 
## 0.1 0.1 0.3 0.3 0.1 0.1
barplot(table(g.sim)/length(g.sim))

set.seed(10)
rgeom(10,0.2)
##  [1]  0  2  2  1  1  9 11  1  6  1

approximate probabilities P(X=4)=0.3 P(X=12)=0.1 Did you get the same answer? set.seed(10); rgeom(10, 0.2)??? NO, uses diff methods

———Generating a Poisson Random Variable———

let x~ Poisson (mean=4) and P(x=0)= 0.01832 i = 0 START to compute P(x=1), P(x=2) and P(x=3).

dpois(0:3,4)
## [1] 0.01831564 0.07326256 0.14652511 0.19536681

2)Consider Pois(2) and seed (123), let’s generate 3 observations manually

cdf = p = dpois Po= at 0

set.seed(123)
x= 0:6 ; dpois(x,2)
## [1] 0.13533528 0.27067057 0.27067057 0.18044704 0.09022352 0.03608941 0.01202980
ppois(x,2)
## [1] 0.1353353 0.4060058 0.6766764 0.8571235 0.9473470 0.9834364 0.9954662
cumsum(dpois(x,2))
## [1] 0.1353353 0.4060058 0.6766764 0.8571235 0.9473470 0.9834364 0.9954662
u1 = runif(1); u1
## [1] 0.2875775
u2 = runif(1); u2
## [1] 0.7883051
u3 = runif(1); u3
## [1] 0.4089769
set.seed(123)
rpois(3,2)
## [1] 1 3 2

if u1 < cdf0 (p0) u1 =0.2875775 > cdf0 = 0.1353353 NOO then: check u1= 0.2875775 < cdf1 (p0 + p1) = 0.4060058 YESS so x1=1

u2= 0.7883051 x2=3 u3= 0.4089769 x3=2

set.seed(123) rpois(3,2) same 1 3 2

——- while loop ——–

i <- 1
while(i<6){
print(i)
i =i+1
}
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
y <- 10
while(y>3){
  y <- y -1.5
  print(y)
}
## [1] 8.5
## [1] 7
## [1] 5.5
## [1] 4
## [1] 2.5
  1. Generate 10 observations from Poisson(lambda=4) using the ITM.
set.seed(10); n=10; lambda= 4; u=runif(n); xsim=c()

for(j in 1:n){
  
  i=0   #setting i=0

  p= exp(-lambda)   #setting p(x=0)=exp(-lambda)
  cdf= p     #setting cdf=pp(x=0)

  while (u[j]>= cdf) {
    p=(lambda*p)/(i+1)  #updating p, cdf and i
    cdf= cdf+p ; i= i+1
  }
  xsim[j]=i
}
xsim
##  [1] 4 3 3 5 1 2 3 3 4 3
mean(xsim)
## [1] 3.1
table(xsim)/length(xsim)
## xsim
##   1   2   3   4   5 
## 0.1 0.1 0.5 0.2 0.1
set.seed(10)
rpois(10,4)
##  [1] 4 3 3 5 1 2 3 3 4 3

——–Generating Binomial Random Variables———

  1. Generate 5 observations from the binomial distribution with p=0.35 and n=10.

run rbinom(5, 10, 0.35) and seed=15

set.seed(15)
m=5 #Number of runs.
u=runif(m)
x.bin=c()
for(j in 1:m){
  n=10 #binomial sample size
  p=0.35 #setting p=0.35
  i=0 #setting i=0
  prb=(1-p)^n
  c=p/(1-p)
  cdf=prb #setting cdf=prb (p(x=0))
  while(u[j]>=cdf){
    prb=(c*((n-i)/(i+1)))*prb
    cdf=cdf+prb
    i=i+1
  }
  x.bin[j]=i
  }
x.bin
## [1] 4 2 6 4 3
barplot(table(x.bin)/length(x.bin),col="red")

set.seed(15)
rbinom(5, 10, 0.35)
## [1] 4 2 6 4 3

——-Generating Samples Using Acceptance-Rejection Method ARM ———— “Target distribution” “proposed distribution” = pmf

p(y)/q(y) ≤ c

Generate n observations from this distribution using the Acceptance-Rejection Method?

proposed distribution >> discrete uniform[1, 5] = q(y)= 1/5 , y= 1:5 prob.x is simply P(X=x) prob.y is simply P(Y=y) proposed distribution >> always be (1/5) because Y~DiscUnf[1, 5].

set.seed(123)
n.gen= 10
 x=1:5
prop.x=c(.05,.25,.45,.15,.10) #X-distribution
prop.y = rep((1/5),5) #y distribution proposed distribution
c = max(prop.x/prop.y); c
## [1] 2.25
sim.vector= rep(0,n.gen)

for (j in 1:n.gen) {
  u1 = runif(1)
  #Generating Y from disct. Unif[1,5] with pmf = 1/5
  y= floor(5*u1) +1 
  
  
  k= c* (prop.y[y])
  u2 = runif(1)
  
  #compute ratio= prop.x[y]/ c*prop.y[y]
  rt= prop.x[y]/ k
  while (u2>rt) {
    #To take care if we reject so it will search for the
    #next possible value to be accepted.
    u1 = runif(1)
    y= floor(5*u1)+1
    u2= runif(1) #You need to add this inside the loop

   
  }
  sim.vector[j] = y
  
}
sim.vector
##  [1] 5 3 3 2 2 4 2 4 2 4

c=2.25

  1. X takes values 1, …, 10 with probabilities 0.11, 0.12, 0.09, 0.08, 0.12,0.10, 0.09, 0.09, 0.10, 0.10.

generate 20 values

x= 1:10
set.seed(123); n.genr=20 
propx= c(0.11, 0.12, 0.09, 0.08, 0.12, 0.10, 0.09, 0.09,
0.10, 0.10)

propy= rep( (1/10) ,10)

c= max(propx/propy) ;c
## [1] 1.2
sim.vec= rep(0, n.genr)

for (j in 1:n.genr) {
  u1= runif(1)
  y= floor(u1 *10)+1
  
  k= (propy[y]) *c
  u2= runif(1)
  rt= (propx[y]/(propy[y])*c)
  
  while (u2>rt) {
    u1= runif(1)
    y= floor(10*u1)+1
    
    u2= runif(1)
    
  }
  sim.vec[j]= y
  
}
sim.vec
##  [1]  3  5 10  6  6 10  7  2  3  4  9  7  7  6  3 10  7  1  8  4