Problem 4

#data
x <- c(1,2,5)
p <- c(0.1,0.3,0.6)
df <- data.frame(x,p)
df
##   x   p
## 1 1 0.1
## 2 2 0.3
## 3 5 0.6
cdf <- data.frame(x=c(1,2,5), prob=c(0.1,.3,.6))

#cdf
plot(stepfun(cdf$x,c(0,cdf$prob)))

#simulate function 

distrib4 <- function(n, p){
U <- runif(n)
X <- rep(0,n)
w1 <- which(U <= p[1])
X[w1] <- 1
w2 <- which( (U > p[1]) & (U < sum(p[1:2])) )
X[w2] <- 2
w3 <- which( U > sum(p[1:2]) )
X[w3] <- 5
return(X)
}

X = distrib4(10000, c(.1, .3, .6))
mean(X == 1)
## [1] 0.0955
mean(X == 2)
## [1] 0.3046
mean(X == 5)
## [1] 0.5999

Problem 11 a.

The PMF = \((\prod_{i=2}^{y} \frac{1}{i}) * (\frac{y}{y+1})\)

N is the number of times we go around the while loop when Y.sim is called. So, EN is the expectation of N which is found when we run the function many times. The expectation we obtain is 12.056. To obtain the time, we use system.time for one simulation and see that the simulation takes 0.017 to run once.

So to obtain how long it takes to run we multiply the expected value by the run time for one, which would be 0.20496237.

Y.sim <- function() {
U <- runif(1)
Y <- 1
while (U > 1 - 1/(1+Y)) {
Y <- Y + 1
}
return(Y)
}

set.seed(123456)
expectation <- replicate(100000, Y.sim())
mean(expectation)
## [1] 9.21011
options(digits=10)
system.time(Y.sim())
##    user  system elapsed 
##       0       0       0
0.017*12.05661
## [1] 0.20496237

Problem 12

Part A - Poisson distribution with lambda=10

Part B - Geometric with prob of success=0.8

Part C - Code below

Part D - The proportion estimate is about 0.1361. The theoretical is e^-2 = 0.1358. We also calculate the prop confidence interval and both the estimate and theoretical are found in the the interval!

#Part A 
#Poisson distribution with lambda=10
N<-rpois(1,lambda=10)

#Part B
#Geometric with prob of success=0.8
xi<-rgeom(1,prob=0.8)

#Part C
N<-rpois(1,lambda=10)
xi<-rgeom(1,prob=0.8)

shoes <- function(){
  N <- rpois(1,10)
  Xi <- rgeom(N,prob=0.8)
  sum(Xi)
}
shoes()
## [1] 1
#Part D
set.seed(1234)
k <-replicate(100000,shoes())
prop <- length(k[k == 0])/length(k)
prop
## [1] 0.1337
#Theoretical
exp(-2)
## [1] 0.1353352832
prop + 1.96 * sqrt( (prop*(1-prop))/length(k))
## [1] 0.1358093854
prop - 1.96*  sqrt( (prop*(1-prop))/length(k)) 
## [1] 0.1315906146

Problem 13

#CDF Plot
 samplepdf <- function(x) { 
3 * (x - 1) ** 2
}
samplecdf <- function(x) {
sum(samplepdf(seq(1, x, .00001)) * .00001) 
}
plot(c(0, seq(1, 2, .01), 3), c(0, lapply(seq(1, 2, .01), samplecdf), 1), type = 'l', xlim = c(0,3), ylim = c(0,1), xlab = 'X-Values', ylab = 'Cumulative Probability')

#Simulate from CDF
X.sim <- function(S)1+(S)^(1/3)
S <- runif(100000)   
x <- X.sim(S)
max(x)
## [1] 1.999998384

Problem 14

PDF: exp(-x)/(1+exp(-x))^2

CDF: y = 1/(1 + e^-x)

Inverse CDF: x = -ln(1/y - y)

num.samples <-  1000
U           <-  runif(num.samples)
X           <- -log(1/U - U)

#pdf
pdfnorm <- function(x, mu = 0, sigma = 1){
out <- 1/(sigma*sqrt(2*pi)) * exp(-1/(2*sigma^2)*(x - mu)^2)
return(out)
}
#cdf
cdfnorm <- function(x, mu = 0, sigma = 1){
sum(pdfnorm(seq(-10*sigma,x,0.00001))*0.00001)
}
#inverse cdf
icdfnorm <- function(p, mu = 0, sigma = 1){
#Objective function
f <- function(x,p){
abs(cdfnorm(x) - p)
}
optimize(f, lower = -10, upper = 10, p=p)$minimum
}
icdfnorm <- Vectorize(icdfnorm)
u <- runif(10)
icdfnorm(u)
##  [1]  0.9083914093 -1.1866378655 -0.5398417502  1.6245927650  0.8670706010
##  [6]  1.8122145636  0.5590457789 -0.2465484052  0.3076418881 -0.0180162682
# plot
hist(X, freq=F, xlab='X', main='Generating Exponential R.V.')
curve(dexp(x, rate=2) , 0, 3, lwd=2, xlab = "", ylab = "", add = T)

Problem 17

Part A.

Part B.

Part C.

Part D.

#Inverted CDF Plot
set.seed(1234)
inverty <- function(y){sqrt(3/(3-2*y))-1}

u <- runif(10000)
y <- inverty(u)
plot(y,u)

Problem 18

#For alpha = 1 
x = c()
for (i in 1:1000){
u = runif(1,0,1)
x[i] = tan(u-0.5)
}

hist(x)

cauch_func <- function(alpha){
  f <- function(x) (1/sqrt(2*pi))*exp(-x^2/2)
  h <- function(x) alpha/(pi * (alpha^2 + x^2))
  kstar<- (exp(-1+(alpha/2)))/sqrt(2*pi)*((2*pi)/alpha)
  while(TRUE){
    X <- alpha*tan(pi*runif(1))
    Y <- runif(1,0,kstar*h(X))
    if(Y < f(X) )
     return(X)
  }
}

cauch_func(12)
## [1] -1.366987842