DM2HW4

Question 1:

#Actual true probability 
print("Actual True Probability")

## [1] "Actual True Probability"

pnorm(3,mean = 0,sd=1,lower.tail = FALSE,log.p = FALSE)

## [1] 0.001349898

#Setting seed to replicate results
set.seed(123)
#Sampling 10,000 observations from population using normal dist.
samp=rnorm(10000,mean=0,sd=1)
#Extracting onservations greater than 3
newdata <- samp[which(samp > 3)]
#Probability of sample observations greater than 3
print("Probability of P(Z>3)")

## [1] "Probability of P(Z>3)"

length(newdata)/length(samp)

## [1] 0.0015

#Variance of sample
print("Variance of Sample")

## [1] "Variance of Sample"

var(samp)

## [1] 0.9972751

Question 2:

#Setting up library for importance sampling
library(truncnorm)
#Setting seed to replicate results
set.seed(123)
#Sampling ONLY observations greater than 3
samp=rtruncnorm(n=10000, a=3, mean=0, sd=1)
#Probability of sample observations greater than 3
print("Probability of P(Z>3)")

## [1] "Probability of P(Z>3)"

length(samp)/length(samp)

## [1] 1

#Variance of sample
print("Variance of Sample")

## [1] "Variance of Sample"

var(samp)

## [1] 0.06956769

Question 3:

# metropolis for N(0,1) based on uniform candidates
norm<-function (n, alpha) 
{
        vec <- vector("numeric", n)
        x <- 0
        vec[1] <- x
        for (i in 2:n) {
                innov <- runif(1, -alpha, alpha)
                can <- x + innov
                aprob <- min(1, dnorm(can)/dnorm(x))
                u <- runif(1)
                if (u < aprob) 
                        x <- can
                vec[i] <- x
        }
        vec
}

#Setting seed to replicate results
set.seed(123)
normvec<-norm(5000,1)
par(mfrow=c(2,1))
plot(ts(normvec))
hist(normvec,30)

par(mfrow=c(1,1))
#Extracting onservations greater than 3
newdata <- normvec[which(normvec > 3)]
#Probability of sample observations greater than 3
print("Probability of P(Z>3)")

## [1] "Probability of P(Z>3)"

length(newdata)/length(samp)

## [1] 0.0023

#Variance of sample
print("Variance of Sample")

## [1] "Variance of Sample"

var(normvec)

## [1] 1.02479