R Notebook

library(binom)
library(collapsibleTree)
library(dbplyr)
library(devtools)
library(dplyr)
library(EnvStats)
library(ggformula)
library(ggplot2)
library(ggpubr)
library(htmltools)
library(ISLR)
library(knitr)
library(lawstat)
library(markdown)
library(mosaic)
library(mdsr)
library(mosaicData)
library(nycflights13)
library(plyr)
library(purrr)
library(rmarkdown)
library(stringi)
library(tibble)
library(tidyr)
library(tidyselect)
library(tinytex)
library(yaml)

Q1 a.

0.4 * 0.4

[1] 0.16

\[ P(X=2) = 0.1600 \] b.

(1-0.4) * (1-0.4)

[1] 0.36

\[ P(X=0) = 0.3600 \] c.

choose(2,1)*0.4*(1-0.4) + choose(2,2)*0.4*(1-0.4)

[1] 0.72

\[ P(X\geq1) = 0.7200 \] d. \[ P(1 \leq X \leq n) = 0.95 \] \[ Then \; P(X=0) = 1 - P(1 \leq X \leq n) = 0.05 \]

\[ P(X=0)=0.4^n = 0.05 \] \[ Then \; \ln{0.4^n} = n\ln{0.4} = \ln{0.05} \] \[ Then \; n = \ln{0.05} / \ln{0.4} \]

log(0.05) / log(0.4)

[1] 3.269412

\[ n = 3.2694 \approx 4 \] Q2

Step 1: Create a series of vectors to hold output

nsims = 1000  # the number of simulations
outcome = numeric(nsims) #create empty vector that will be filled with numeric outcomes

Step 2: Run the Simulation

for(i in 1:nsims){ #we are going to perform the body of the loop 1000 times
  outcome[i] = sample(c(1,2,3,4,5,6), 1, replace=FALSE) #stores the outcome of the ith toss in position i of outcome
  }  #close the for-loop
simresult = data.frame(outcome) #creates a data frame with two columns
head(simresult,3)  #shows the first three rows of the data frame, outcome of each die toss in column 2

tail(simresult,3)  #shows the last three rows of the data frame, die toss outcome in column 2

Step 3: Visualize the Simulation with ggplot2

sum(simresult$outcome)  #computes the sum of the 1000 tosses appearing in the outcome variable of simresult data frame

[1] 3580

ggplot(data = simresult, aes(x = outcome, fill = outcome)) + geom_histogram(col = 'black', fill = "red", binwidth = 1)

Step 4: Using Simulation to Compute a Probability

nsims = 3000  # the number of simulations
outcome = numeric(nsims) #create empty vector that will be filled with numeric outcomes
fivesix = numeric(nsims) #create empty vector that will be filled with outcomes of five or six 

for(i in 1:nsims){ #we are going to perform the body of the loop 1000 times
  outcome[i] = sample(c(1,2,3,4,5,6), 1, replace=FALSE) #stores the outcome of the ith toss in position i of outcome
  fivesix[i] = if (outcome[i] == 5 || outcome[i] == 6) 1 else 0 #stores the outcome of five or six 
  }  #close the for-loop
simresult = data.frame(outcome) #creates a data frame with two columns
head(simresult,3)  #shows the first three rows of the data frame, outcome of each die toss in column 2

simresultfivesix = data.frame(fivesix)
head(simresultfivesix, 10)

nsims = 3000  # the number of simulations
outcome1 = numeric(nsims) #create empty vector that will be filled with numeric outcomes
toss1 = numeric(nsims) #create empty vector that will be filled with outcomes of toss1
toss2 = numeric(nsims) #create empty vector that will be filled with outcomes of toss2
toss3 = numeric(nsims) #create empty vector that will be filled with outcomes of toss3

for(i in 1:nsims){ #we are going to perform the body of the loop 1000 times
  toss1[i] = sample(c(1,2,3,4,5,6), 1, replace=FALSE) #stores the outcome of the ith toss1 in position i of outcome
  toss2[i] = sample(c(1,2,3,4,5,6), 1, replace=FALSE) #stores the outcome of the ith toss2 in position i of outcome
  toss3[i] = sample(c(1,2,3,4,5,6), 1, replace=FALSE) #stores the outcome of the ith toss2 in position i of outcome
  outcome1[i] = if (toss1[i]+toss2[i]+toss3[i] >= 14 ) 1 else 0 #stores the outcome of five or six 
  }  #close the for-loop

simresult1 = data.frame(outcome1) #creates a data frame with two columns
head(simresult1,10)  #shows the first three rows of the data frame, outcome of each die toss in column 2

tail(simresult1, 10)

sum(simresult1$outcome1)  #computes the sum of the 3000 tosses appearing in the outcome variable of simresult1 data frame

[1] 495

495 / 3000

[1] 0.165

\[ P(Sum \geq 14) = 0.165 \]

Q3

1. Compute the probability that your hand will have neither ♠s nor ♣s.

(choose(10,1)*choose(10,4) + choose(10,2)*choose(10,3) + choose(10,3)*choose(10,2) + choose(10,4)*choose(10,1) + choose(10,5)) / choose(20,5)

[1] 0.9837461

\[ P(♠s \; or \; ♣s) = 0.9837461 \]

2. Compute the probability your hand will consist of a 10, Jack, Queen, King, and Ace of the same suit. One example of such a hand is: (10♡,J♡,Q♡,K♡,Ace♡).

 4 / choose(20,5)

[1] 0.0002579979

\[ P(suits) = 0.0002579979 \approx 0.0003 \]

3. Compute the probability that you get a three-of-a-kind. For example, (10♡, 10♢, 10♠, J♣, K♡).

choose(4,3)*choose(16,2)*5 / choose(20,5)  #calculate probability all three of a kind hands, times 5 because there are 10, Jack, Queen, King, Ace, 5 kinds

[1] 0.1547988

\[ P(three of a kind) = 0.1547988 \approx 0.1548 \]

4. Compute probability that one observes two Aces and two ♢’s.

choose(4,2)*choose(5,2)*choose(11,1) / choose(20,5)

[1] 0.04256966

\[ P(two Aces and two ♢s) = 0.04256966 \approx 0.0426 \] Q4

1. The executive has called the office of the lobby-group to say they have missed their connecting flight. Compute the probability that the executive called from Chicago (or is flying Airline UA).

\[ P(AA) = 0.15 \;\; P(UA) = x \;\; P(D) = 3x \] P(AA) + P(UA) + P(D) = 1, Then 0.15 + x + 3x = 1

(1-0.15) / 4 #solving x

[1] 0.2125

\[ P(UA) = 0.2125 \]

0.2125 * 3 #determine P(D)

[1] 0.6375

\[ P(D) = 0.6375 \]

(0.2125*0.3) / ((0.2125*0.3) + (0.15*0.15) + (0.6375*0.1)) #calculate P(Chicago) 

[1] 0.425

\[ P(Chicago) = 0.425 \] b. The executive has not missed a connecting flight and made it to IAD. Compute the probability that they flew on a Delta flight

(0.6375*0.9) / ((0.2125*0.7) + (0.15*0.85) + (0.6375*0.9))

[1] 0.675

\[ P(Minneapolis-St.Paul) = 0.675 \] c. Provide a statement that interprets the meaning of the probability computed in (b) in the context of these data.

Ans: The executive will most likey flying with Detal airline because the probability is over 50% and Detal airline has the best chance of getting on the connecting flight and arriving in IAD.

Q5.

a.Using R Studio, create a display that shows the probability distribution of this particular random variable X. Refer to the various code provided for examples appearing in both Probability Module 4 and Review Exercise 5 from Thursday, September 10th. For values of x, use xvalues = 0:15.

pdf1 <- function(x) # create a probablity distribution function of 2 / 3^(x+1)
{
  2 / (3^(x+1))
}

curve(pdf1, 0, 15, xlab="Values of X", ylab="Values of pdf1(x)", main="Distribution of X", col='red')

2. Compute P(X>3)

integrate(pdf1, lower=3, upper=15)$value

[1] 0.022475

\[ P(X>3) = 0.022475 \] c. Compute the mean or expected value of X, E(X) or μX. (Hint: In computing E(X), change the upper limit on xvalues from 15 to 50…)

pdf1mean <- function(x) #create a funiction of x*pdf1(x) so that it can be integrated using integrate commande with lower=0 and upper=50
{
  x*pdf1(x)
}

integrate(pdf1mean, lower=0, upper=50)$value

[1] 0.552357

\[ E(X) \; or \; \mu{x} = 0.552357 \] d. Compute the standard deviation of X, SD(X) or σX. (use the same values of x as you did in part (c))

pdf1sd <- function(x) #create a function of x^2 * pdf1(x) so that it can be integrated using integrate commande with lower=0 and upper=50
{
  x*pdf1mean(x)
}

sqrt(integrate(pdf1sd, lower=0, upper=50)$value - (integrate(pdf1mean, lower=0, upper=50)$value)^2)

[1] 0.8369323

\[ SD(X) \; = \; 0.8369323 \]

e.Consider the interval (μx−σx, μx+σx). Compute P(μx−σx <X< μx+σx)

0.552357 - 0.8369323

[1] -0.2845753

\[ \mu{x} - σx = -0.2845753 \]

integrate(pdf1, lower=(0), upper=(0.552357+0.8369323))$value

[1] 0.4749383

Q6.

pdf2 <- function(n,x) #create the function for the probability of number of canadians disagree with internet companies 
{
  choose(n,x)*(0.6^x)*(0.4^(n-x))
}

pdf2(40,30) #calculate exactly 30 will disagree with internet companies handing over private information to authorities

[1] 0.01964986

p2835 = 0           #define a variable stores P(28 =< x <= 35)
for (i in 28:35)    #for loop to calculate P(x = 28), P(x = 29)...P(x = 35) and sum them up
{
  p2835 = p2835 + pdf2(40,i) #sum up P(x = 28), P(x = 29)...P(x = 35) and stores the total in p2835
}
p2835

[1] 0.1284811

p15 = 0
for (i in 15:25)    
{
  p15 = p15 + pdf2(i,15) 
}
p15

[1] 0.868853

\[ P(15th \; person \; disagree) = 0.868853 \]

OddPersonOut <- function(x)
{
  0.6875^(x-1)*0.3125
}

OddPersonOut(10)

[1] 0.0107227

\[ P(x=10) = 0.0107227 \]

OddPersonOutMean <- function(x)
{
  x * OddPersonOut(x)
}

integrate(OddPersonOutMean, lower=1, upper=50)$value

[1] 3.059875

\[ E(X) = 3.059875 \; trials \]

Ans: when x >=50, expected trails to obersve “odd person out” converges to 3.059875, then yes, 10 trails is more than expected to observe “Odd person out”

Q8.

pnorm(1.91, 1.89, 0.05) - pnorm(1.83, 1.89, 0.05)

[1] 0.5403521

\[ P(1.83 \leq x \leq 1.91) = 0.5403521 \]

qnorm(0.9, 1.89, 0.05)

[1] 1.954078

Interpretation: 90% of 2-liter bottles have 1.954078 liter or less soft drink filled.

1 - pnorm(2, 1.89, 0.05) #calculate the probability of 2-liter botters being overfilled

[1] 0.01390345

\[ P(x > 2) = 0.01390345 \]

pnorm(1.85, 1.89, 0.05) #calculate the probability of 2-liter bottle has less than 1.85 liters of soft-drink filled

[1] 0.2118554

1 - pnorm(1.85, 1.89, 0.05) #calculate the probability of a 2-liter bottle has more than 1.85 liters of soft-drink filled

[1] 0.7881446

P510 = 0 
for (i in 5:10)
{
  P510 = P510 + choose(50,i)*(0.2118554^i)*(0.7881446^(50-i))
}
P510

[1] 0.4892489

\[ P( Between \; 5 \; and \; 10 \; bottles \; have \; less \; 1.85 \; liters \; of \; soft \; drink) = 0.4892489 \]

Q9.

gss = read.csv("http://people.ucalgary.ca/~jbstall/DataFiles/GSS2002.csv")

head(gss, 4)

ggplot(data=gss, aes(x = Religion, fill=Religion)) + geom_bar(position="dodge", na.rm=TRUE) + theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

counts=as.data.frame(sort(table(gss$Religion)))
head(counts, 4)

ggplot(data = gss, aes(x = Education, fill = Race)) + geom_bar(position = "dodge", na.rm = TRUE)

Interpretation: White people are majority population in each education level, then followed by black people.

ggplot(data=gss, aes(x = Marital, fill = Education)) + geom_bar(position = "dodge", na.rm=TRUE) 

The plot above shows that each marital category has similar distribution in each education level. Majority of population has HS, then Bachelors, Left HS, Graduate,…, only exception is Widowed, it has Left HS come in as second.

ggplot(data=gss, aes(x = Education, fill = Marital)) + geom_bar(position = "dodge", na.rm=TRUE)

Second plot shows each education level has similar distribution in each marital category. Majority people in each education level are married and followed by Never Married, Divorced, Widowed, etc.

ggplot(data=gss, aes(x = Gender, fill=Politics)) + geom_bar(position = "dodge", na.rm=TRUE) 

It seems like Female and Male has similar distribution across different politics. Majority of them are Moderate, and followed by Conservative, Slightly conservative, Slightly liberal. Wemon has higher probability than men to be conservative(including Extremely conservative, Conservative, slightly conservative)

library(ISLR)
head(Default, 4) #shows the first four lines of the Default data set

ggplot(data=Default, aes(x = balance, y = income, color = student)) + geom_point(size = 2, position = "jitter") + xlab("Balance") + ylab("Income") + ggtitle("Monthly Income vs Monthly Credit Card Balance")

ggplot(data=Default, aes(x= student, y=balance, color = student)) + geom_boxplot(na.rm=TRUE) + xlab("Student Status") + ylab("Monthly Balance")

mean(~balance, data=Default, na.rm=TRUE)   #calculate the mean of balances from the data frame Default

[1] 835.3749

median(~balance, data=Default, na.rm=TRUE)  #calculate the mean of balances from the data frame Default

[1] 823.637

sd(~balance, data=Default, na.rm=TRUE)

[1] 483.715

balancedata = c(Default$balance)  #create a vector that stores all balances from data frame Default
quantile(balancedata, 0.05)   #calculate 5th percentiles of all balances

       5% 
0.4246597 

quantile(balancedata, 0.95)

     95% 
1665.963 

Q11.

courier <- read.csv('Data602Assignment1Question11.csv', header = TRUE)
courier = courier %>% mutate(Courier = "Local Courier")
courier

ggplot(data = courier, aes(x = Courier , y = Delivery_time)) + geom_violin(col = "blue", fill = "red", na.rm = TRUE) + xlab("Courier") + ylab("Delivery time") + ggtitle("Voilin plot of this courier delivery time")

mean(~Delivery_time, data = courier, na.rm = TRUE)

[1] 5.6875

median(~Delivery_time, data = courier, na.rm = TRUE)

[1] 5.775

sd(~Delivery_time, data = courier, na.rm = TRUE)

[1] 1.580369

qdata(~Delivery_time, c(0.25, 0.75), data = courier)

  25%   75% 
4.750 6.565 

qdata(~Delivery_time, 0.99, data = courier)

   99% 
7.9778 

qnorm(0.99, 5, 1.5)

[1] 8.489522

Because the amount of delivery time can be modeled by the Normal distribution with a mean delivery time of 5 hours and a standard deviation of 1.2 hours, so 99th percentile is 8.489522 hours, which means all 99% of packages will be delivered within 8.489522 hours or less. Anything is over this delivery time will be refunded.