HW 6
Ali Nemati
Last compiled on 19 April, 2022 , 00:23:38
Answer question 1
a
A fair coin is flipped 100 times. Find the probability that the number of times it shows heads is exactly equal to 40
options(warn=-1) # remove warning with -1, show waring with 0
############################################################################
############################################################################
### ###
### Q1-a ###
### ###
############################################################################
############################################################################
set.seed(1984) # My birth date LOL
goal <- 100 # set goal
# initiate with null
heads <- rep("Na", goal) # create heads with na
i<- 0 # initiate i
dt <- c("heads", "tails") # vector dt
#################################################################
## while i<= goal => True. ##
#################################################################
while (i<= goal) { # while i<=1000
# print(i)
flips <- sample(dt, goal, replace=TRUE) # sample takes a sample of the specified
#size from the elements of x using either with or without replacement.
#print(flips)
heads[i] <- length(which(flips == "heads")) / goal # len/1000 if which true
i = i+1 # counter
}
flips[1:5]## [1] "heads" "heads" "tails" "tails" "heads"plot(density(as.numeric(heads)))
head(heads)## [1] "0.5" "0.49" "0.47" "0.5" "0.49" "0.45"sum(length(which(heads == 0.4))) / goal ## [1] 0.01# sum of all result from for/1000 if which TRUe
# display data ggplot
df <- data.frame(flips, heads)
library(ggplot2)
# Note: stat = "identity" tells ggplot to use the values of the y-axis variable (p) as the height of the bars in our histogram (as opposed to counting the number of occurances of those values)
ggplot(df, aes(x = flips,y = heads)) +
geom_bar(stat = "identity", fill = "purple") +
labs(x = "Number of Heads", y = "Probability in 100 flips (p)") +
theme_minimal()
b.
An urn contains 10 balls. Two red, three white, and five black. All 10 are drawn one at a time without replacement. Find the probability that both the first and last balls drawn are black.
############################################################################
############################################################################
### ###
### Q1-b ###
### ###
############################################################################
############################################################################
# initiate vac
urn <- c(rep("Black", 5), rep("White", 3), rep("Red", 2))
# time loop should run
goal <- 10000
# create our data to save data.frame
data <- structure(list(col1 = character(),
col2 = character(),
black = integer()),
class = "data.frame")
datai<-0 # inintiate i = 0
drawn <- 10 # inintiate drawn =10
set.seed(1984) # My birth date LOL
#################################################################
## while i<= goal => True. ##
#################################################################
while(i<= goal) {# while i<= 1000 True
sample_urn<- sample(urn, drawn, replace = FALSE) # All 10 are drawn one at a time without replacement
data[i, 1:2] <- sample_urn[c(1, drawn)]
data[i, 3] <- ifelse(data[i, "col1"] == "Black" &
data[i, "col2"] == "Black", 1, 0)
i =i+1
}
tail(data)plot(density(data$black))
data['is_black_prob']<-data['black'] / goal # data / 10000
sum(data['is_black_prob']) # our result 0.217## [1] 0.2236c.
############################################################################
############################################################################
### ###
### Q1-c ###
### ###
############################################################################
############################################################################
set.seed(1984) # My birth date LOL
# make an n_runs x n_flips matrix of 0s and 1s
# add function getting nsim, nflips, nheads, prob
#################################################################
## add function getting nsim, nflips, nheads, prob. ##
#################################################################
coins <- function(nflips = 100, nsim = 1000, nheads = 5,prob = 0.53){
flips <- matrix(rbinom(nflips * nsim, 1, prob), nrow = nflips)
#print(flips)
##################################################################
## To determine the probability, divide the winnings ##
## by nsim if they exceed the amount. ##
##################################################################
flips <- ifelse(flips==0, -1,1) # if 0 => -1, else 1
#
# result
return (sum(colSums(apply(flips, 2, cumsum) >= nheads) > 0) / nsim)
}
# run test if
coins(nflips = 100, nsim = 1000) # our result 0.801## [1] 0.801# run test if
coins(nflips = 100, nsim = 2000, nheads = 5,prob = 0.55)# our result is 0.8735## [1] 0.8735Question 2
a
############################################################################
############################################################################
### ###
### Q2-a ###
### ###
############################################################################
############################################################################
set.seed(1984) # My birth date LOL
b_median <- function(data, num = 30) {
### Generate 30 bootstrap samples.
resamples <- lapply(1:num, function(i) sample(data, replace=T))
### Calculate the median for each bootstrap sample
r.median <- sapply(resamples, median)
### sqrt of r.median
std.err <- sqrt(var(r.median))
# return as list with std.err
list(std.err=std.err, resamples=resamples, medians=r.median)
}
#################################################################################
## In the same manner as in the example above, generate the data to be used. ##
#################################################################################
data1 <- round(rnorm(100, 5, 3))
##################################################################
## saving the results of the function median in the object b. ##
##################################################################
b <- b_median(data1, 30)
#################################################################
## displaying the first of the 30 bootstrap samples. ##
#################################################################
b$resamples[1]## [[1]]
## [1] 8 10 0 3 4 12 5 7 5 0 4 -1 5 9 5 2 0 3 3 0 4 7 1 2 6
## [26] 8 6 1 7 6 3 8 5 7 1 1 5 9 5 9 5 4 7 -1 8 9 5 7 3 6
## [51] 5 5 5 2 13 8 6 0 9 9 6 4 0 5 4 6 6 9 3 8 5 8 -1 9 4
## [76] 5 4 3 0 4 9 1 8 9 0 3 4 -1 5 3 5 3 4 5 4 7 1 8 0 2##################################################################
## displaying the standard error. ##
##################################################################
b$std.err## [1] 0.4722166##################################################################
## displaying the histogram of the distribution of medians. ##
##################################################################
plot(density(b$medians))
hist(b$medians)
b
Write a simulation with 1,000 iterations that evaluates the coverage properties of this CI (i.e., in what percent of simulations was the true median of the chisq distribution included in the 95% CI
set.seed(1984) # My birth date LOL
############################################################################
############################################################################
### ###
### Q2-b ###
### ###
############################################################################
############################################################################
# initiate with null
result <- NULL
goal <- 1000
n<- goal/10
i<-0
#################################################################
## while i<= goal => True. ##
#################################################################
while(i<=goal){
############################################################################
## Chisquare: The (non-central) Chi-Squared Distribution. Density, ##
## distribution function, quantile function and random generation for ##
## the chi-squared (chi^2) distribution with df degrees of freedom and ##
## optional non-centrality parameter ncp. ##
############################################################################
y_rchisq <- rchisq(n=n, df=2)
# Calc mean
mean_y_rchisq <- mean(y_rchisq)
# initiate with null
result_resample <- NULL
for(j in 1:goal){
y_resample_y_rchisq <- sample(y_rchisq, replace=TRUE, size=n)
result_resample[j] <- 1/mean(y_resample_y_rchisq)
}
##################################################################
## The generic function quantile produces sample quantiles ##
## corresponding to the given probabilities. The smallest ##
## observation corresponds to a probability of 0 and the ##
## largest to a probability of 1. ##
## . ##
##################################################################
# The percentile method is often used to provide an approximate 95% confidence
# interval for the population parameter. The percentile method is
#not as accurate as
# the BCa method, but it is very quick and easy to calculate. Assume that the vector
# "invmean_resample" contains a large number of bootstrap
#replicate estimates of the statistic of
# interest. By the percentile method, a 95% boostrap confidence interval for the
# parameter is obtained as
LowerBond <- quantile(result_resample, probs=c(0.025))
UpperBond <- quantile(result_resample, probs=c(0.975))
result[i] <- ifelse(LowerBond <= 0.5 & UpperBond >= 0.5, 1, 0)
i=i+1
}
plot(density(result))
mean(result)## [1] 0.925Ref
The best way to predict the future is to create it.” Abraham Lincoln.
End of Document.