1. (Bayesian). A new test for multinucleoside-resistant (MNR) human immunodeficiency virus type 1 (HIV-1) variants was recently developed. The test maintains 96% sensitivity, meaning that, for those with the disease, it will correctly report “positive” for 96% of them. The test is also 98% specific, meaning that, for those without the disease, 98% will be correctly reported as “negative.” MNR HIV-1 is considered to be rare (albeit emerging), with about a .1% or .001 prevalence rate. Given the prevalence rate, sensitivity, and specificity estimates, what is the probability that an individual who is reported as positive by the new test actually has the disease? If the median cost (consider this the best point estimate) is about $100,000 per positive case total and the test itself costs $1000 per administration, what is the total first-year cost for treating 100,000 individuals?
#P(Disease/Positive) = P(Positive/Disease) *P(Disease)/P(Positive)

# Given values
sens <- 0.96 #With disease, report positive (True positive)
spec <- 0.98 #Without disease, report negative (True negative). P(false positive) = 1-0.98 = 0.02
prev <- 0.001 # frequency
cost_per_pos <- 100000
cost_per_test <- 1000
n <- 100000


# what is the probability that an individual who is reported as positive by the new test actually has the disease? 
# Probability (positive disease | positive test)
# Bayes Theorem P(A|B) = [P(B|A) * P(A)] /P(B)


#B=Positive test
#A1=Actual disease positive
#A2=Actual disease negative

#P(A1 | B)

#P(A1 | B) =  P(B|A1)P(A1) / P(B|A1)P(A1)+P(B|A2)P(A2)
#P(B|A1)=.96 [sens]
#P(B|A2)=1-.98=.02 [1-spec]
#P(A1)=.001 [prev]
#P(A2)=.999 [1-prev]

Prob_Positive = (sens*prev) / ((sens*prev) + (1-spec)*(1-prev))


# true positives
true_pos = n*prev*sens
total_cost = (true_pos*cost_per_pos) + (n * cost_per_test)


cat("The probability of someone who tested positive, actually having the disease is: ",Prob_Positive,'%\n')
## The probability of someone who tested positive, actually having the disease is:  0.04584527 %
cat("Total first-year cost for treating", n, "individuals: $", total_cost)
## Total first-year cost for treating 1e+05 individuals: $ 109600000
2. (Binomial). The probability of your organization receiving a Joint Commission inspection in any given month is .05. What is the probability that, after 24 months, you received exactly 2 inspections? What is the probability that, after 24 months, you received 2 or more inspections? What is the probability that your received fewer than 2 inspections? What is the expected number of inspections you should have received? What is the standard deviation?
# What is the probability that, after 24 months, you received exactly 2 inspections?
p_2<- dbinom(2, size=24, prob=0.05)  # 0.2232381
cat("Probability of exactly 2 inspections:", p_2 ,"\n")
## Probability of exactly 2 inspections: 0.2232381
# What is the probability that, after 24 months, you received 2 or more inspections? 
p_2_or_more <- pbinom(1, size=24, prob=0.05, lower.tail = FALSE) # 0.3391827
cat("Probability of 2 or more inspections:", p_2_or_more, "\n")
## Probability of 2 or more inspections: 0.3391827
#Alternatively,  1 - pbinom(1, size = 24, prob = 0.05)
#Alternatively,  1-(dbinom(1,24,.05)+dbinom(0,24,.05))




# What is the probability that your received fewer than 2 inspections? 
p_fewer_than_2 <- pbinom(1, size = 24, prob = 0.05) # 0.6608173
cat("Probability of fewer than 2 inspections:", p_fewer_than_2, "\n")
## Probability of fewer than 2 inspections: 0.6608173
#Alternatively, dbinom(1,24,.05)+dbinom(0,24,.05)




# What is the expected number of inspections you should have received? 
expected <- 24 * 0.05 #1.2
cat("Expected number of inspections:", expected, "\n")
## Expected number of inspections: 1.2
# What is the standard deviation?
sd <- sqrt(24 * 0.05 * (1 - 0.05))
cat("Standard deviation:", sd, "\n")
## Standard deviation: 1.067708
# Alternatively, use library(binom)
# p = 0.05 #prob of recive inspection
# n = 24 # time periods
#prob of 2 inspections after 24 months
#prob_2 <- dbinom(2,size=n,prob=p)
#prob 2 or more after 24 months
#prob_2_more <- 1-pbinom(1, size=n, prob=p)
#prob fewer than 2 inspections after 24months
#prob_less_2 <- pbinom(1, size=n, prob=p)
#expected number of inspections
#expected <-n*p
#sd <- sqrt(n*p*(1-p))
3. (Poisson). You are modeling the family practice clinic and notice that patients arrive at a rate of 10 per hour. What is the probability that exactly 3 arrive in one hour? What is the probability that more than 10 arrive in one hour? How many would you expect to arrive in 8 hours? What is the standard deviation of the appropriate probability distribution? If there are three family practice providers that can see 24 templated patients each day, what is the percent utilization and what are your recommendations?
#10 people / hr = lambda
my_lambda <- 10

# What is the probability that exactly 3 arrive in one hour?
p_exactly_3 <- dpois(3, lambda=my_lambda)  #0.007566655
cat("Probability exactly 3 patients arrive in one hour is: ", p_exactly_3, "\n")
## Probability exactly 3 patients arrive in one hour is:  0.007566655
# What is the probability that more than 10 arrive in one hour? 
p_more_10 <- ppois(10, lambda=my_lambda, lower=FALSE) # 0.4169602
cat("Probability more than 10 patients arrive in one hour is: ", p_more_10, "\n")
## Probability more than 10 patients arrive in one hour is:  0.4169602
# Alternatively, # prob_more_10 <- 1- ppois(10, lambda = 10) 


# How many would you expect to arrive in 8 hours? 
exp_8hrs <- my_lambda * 8
cat("Expected number of patients in 8 hours is: ", exp_8hrs, "\n")
## Expected number of patients in 8 hours is:  80
# What is the standard deviation of the appropriate probability distribution? 
sd <- sqrt(my_lambda) # 3.162278
cat("Standard deviation is: ", sd, "\n")
## Standard deviation is:  3.162278
# If there are 3 family practice providers that can see 24 templated patients each day, what is the percent utilization and what are your recommendations?


# 24 people x 3 providers = 72 in 8 hours  # total_cap
total_cap <- 24*3
# 10 people x 8 hrs = 80 which is how many we expect to arrive
exp_px <- exp_8hrs
percent_util <- (exp_px/total_cap)*100
cat("Percent utilization is: ", percent_util, "\n")
## Percent utilization is:  111.1111
4. (Hypergeometric). Your subordinate with 30 supervisors was recently accused of favoring nurses. 15 of the subordinate’s workers are nurses and 15 are other than nurses. As evidence of malfeasance, the accuser stated that there were 6 company-paid trips to Disney World for which everyone was eligible. The supervisor sent 5 nurses and 1 non-nurse. If your subordinate acted innocently, what was the probability he/she would have selected five nurses for the trips? How many nurses would we have expected your subordinate to send? How many non-nurses would we have expected your subordinate to send?
nurses<-15
workers<-30
total_trips<-6
trips_nurses<-5

#Prob selecting 5 nurses
prob_5_nurse<-dhyper(trips_nurses,nurses,workers-nurses,total_trips) #0.07586207
cat("Probability of randomly selecting 5 nurses is: ", prob_5_nurse, "\n")
## Probability of randomly selecting 5 nurses is:  0.07586207
#Alternatively, dhyper(5,15,15,6,log=FALSE)


#Expected number of nurses,  E(X)=KM/N
exp_nurse<-(nurses/workers)*total_trips #3
cat("The expected number of nurses is:", exp_nurse, "\n")
## The expected number of nurses is: 3
#Expected number of non-nurses
exp_not_nurse<-total_trips-exp_nurse #3
cat("The expected number of non-nurses is:", exp_not_nurse, "\n")
## The expected number of non-nurses is: 3
# With such a small chance of randomly choosing 5 numbers, it seems the subordinate is actually biased.
5. (Geometric). The probability of being seriously injured in a car crash in an unspecified location is about .1% per hour. A driver is required to traverse this area for 1200 hours in the course of a year. What is the probability that the driver will be seriously injured during the course of the year? In the course of 15 months? What is the expected number of hours that a driver will drive before being seriously injured? Given that a driver has driven 1200 hours, what is the probability that he or she will be injured in the next 100 hours?
p_crash <- .001

# What is the probability that the driver will be seriously injured during the course of the year?
p_injured_year <- pgeom(1200, p_crash)
cat("Probability of being seriously injured during the year: ", p_injured_year, "\n")
## Probability of being seriously injured during the year:  0.6992876
#Alternatively, p_injured_year <- 1 - (1 - p_crash)^1200



# In the course of 15 months? #15 months is 1500 hrs
p_injured_15_months <- pgeom(1500, p_crash)
cat("Probability of being seriously injured in 15 months: ", p_injured_15_months, "\n")
## Probability of being seriously injured in 15 months:  0.7772602
#Alternatively, p_injured_15_months <- 1 - (1 - p_crash)^(15 * 30 * 24)




# What is the expected number of hours that a driver will drive before being seriously injured?   E[X]=1/p
expected_hours <- 1 / p_crash
cat("Expected number of hours to drive before being seriously injured: ", expected_hours, "\n")
## Expected number of hours to drive before being seriously injured:  1000
# Given that a driver has driven 1200 hours, what is the probability that he or she will be injured in the next 100 hours?
p_injured_in_100_hours <- sum(dgeom(c(1201:1300), p_crash)) # 0.02863018
cat("Probability of being injured in the next 100 hours, if driven 1200 hours: ", p_injured_in_100_hours, "\n")
## Probability of being injured in the next 100 hours, if driven 1200 hours:  0.02863018
#Alternatively, p_injured_in_100_hours <- (1 - (1 - p_crash)^(1200 + 100)) - p_crash
6. You are working in a hospital that is running off of a primary generator which fails about once in 1000 hours. What is the probability that the generator will fail more than twice in 1000 hours? What is the expected value?
# This resembles a Poisson distribution
prob_fail_twice <- 1- ppois(2, lambda = 1) #0.0803014
cat("Probability the generator will fail more than twice: ", prob_fail_twice, "\n")
## Probability the generator will fail more than twice:  0.0803014
#E[X]=lambda =1
E_x <- 1
cat("The expected number of failures in 1000 hours is: ", E_x, "\n")
## The expected number of failures in 1000 hours is:  1
7. A surgical patient arrives for surgery precisely at a given time. Based on previous analysis (or a lack of knowledge assumption), you know that the waiting time is uniformly distributed from 0 to 30 minutes. What is the probability that this patient will wait more than 10 minutes? If the patient has already waited 10 minutes, what is the probability that he/she will wait at least another 5 minutes prior to being seen? What is the expected waiting time?
# We're told this is a uniform distribution between 0 and 30 mins. Probability of waiting more than 10 mins is based on the remaining time (10,30), divided by the total range (0,30). Probability of waiting at least another 5 mins is based on first prob. Expected wait time is average of the upper and lower limits

lower <- 0
upper <- 30



#Prob more than 10 min
prob_10_more= (upper-10)/(upper) 
cat("Probability of waiting more than 10 minutes: ", prob_10_more, "\n") #0.6666667 
## Probability of waiting more than 10 minutes:  0.6666667
# Alternatively, 1-punif(10,0,30)




#Prob of 5 more mins, if already waited 10 mins
prob_5_more=(20-5)/(20-0)
cat("probability of waiting at least another 5 minutes after waiting 10 minutes: ", prob_5_more, "\n") # 0.75 
## probability of waiting at least another 5 minutes after waiting 10 minutes:  0.75
# Alternatively, # Bayes Theorem; conditional P(A|B)/P(B)
#P(A&B)/PA=PA/PB
#PA= 1-CDF P(15)
#PB =CDF P(10)
#PA<-1-punif(15,10,30)
#PB<-1-punif(10,10,30)
#PA/PB # 0.75



#Expected wait
expected_min <- (upper-lower)/2  #15
cat("Expected waiting time: ", expected_min, "minutes.")
## Expected waiting time:  15 minutes.
8. Your hospital owns an old MRI, which has a manufacturer’s lifetime of about 10 years (expected value). Based on previous studies, we know that the failure of most MRIs obeys an exponential distribution. What is the expected failure time? What is the standard deviation? What is the probability that your MRI will fail after 8 years? Now assume that you have owned the machine for 8 years. Given that you already owned the machine 8 years, what is the probability that it will fail in the next two years?
# Geometric distribution, and we're already told the expected failure time of 10 years. The standard deviation of failure is always the same (10 years).



# What is the expected failure time? 
print("Expected Failure Time: 10 years.")
## [1] "Expected Failure Time: 10 years."
# What is the standard deviation? 
# Std dev= Var(Lambda^2)^.5 = lambda = 0.1
print("The standard deviation is always equal to the mean: 10 years.")
## [1] "The standard deviation is always equal to the mean: 10 years."
#Probability MRI fails after 8 years 
P_fail_8yrs <- pexp(8, rate=1/10, lower.tail = FALSE) #0.449329
cat("Probability the MRI will fail after 8 years is: ", P_fail_8yrs, "\n")
## Probability the MRI will fail after 8 years is:  0.449329
# Alternatively, prob_fail_8yrs <- 1 - pexp(8,.1)



#Probability it fails in next 2 years given I owned it for 8 years
P_fail_2yrs <- pexp(8, rate=1/10, lower.tail = FALSE) - pexp(10, rate=1/10, lower.tail = FALSE)
cat("Probability that MRI will fail in the next 2 years after 8 years is", P_fail_2yrs,"\n")
## Probability that MRI will fail in the next 2 years after 8 years is 0.08144952
# Alternatively, prob_fail_2yrs <- exp(-0.1*2)