Discrete & Continuous Probability Distributions

Author

Guibril Ramde

Homework 3

A hospital wants to improve staffing by understanding patient arrivals and treatment times.

#download Arrivals csv file and Treatment csv file


library(tidyverse)

Warning: package 'ggplot2' was built under R version 4.5.2

Warning: package 'tibble' was built under R version 4.5.2

Warning: package 'tidyr' was built under R version 4.5.2

Warning: package 'readr' was built under R version 4.5.2

Warning: package 'purrr' was built under R version 4.5.2

Warning: package 'dplyr' was built under R version 4.5.2

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.0     ✔ readr     2.1.6
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.2     ✔ tibble    3.3.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.2
✔ purrr     1.2.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

arrivals_data <- read.csv("arrivals.csv")
arrivals_data

   Hour Patients
1     1       18
2     2       24
3     3       20
4     4       16
5     5       19
6     6       30
7     7       28
8     8       25
9     9       33
10   10       27
11   11       21
12   12       17
13   13       20
14   14       26
15   15       23
16   16       31
17   17       29
18   18       35
19   19       22
20   20       18

treatments_data <- read_csv("treatment_times.csv")

Rows: 15 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (2): PatientID, TreatmentMinutes

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

treatments_data

# A tibble: 15 × 2
   PatientID TreatmentMinutes
       <dbl>            <dbl>
 1         1               32
 2         2               27
 3         3               45
 4         4               19
 5         5               55
 6         6               40
 7         7               29
 8         8               38
 9         9               44
10        10               22
11        11               31
12        12               36
13        13               47
14        14               24
15        15               52

##Part A: Patient Arrivals

Compute descriptive statistics

#Summary of data
summary(arrivals_data)

      Hour          Patients    
 Min.   : 1.00   Min.   :16.00  
 1st Qu.: 5.75   1st Qu.:19.75  
 Median :10.50   Median :23.50  
 Mean   :10.50   Mean   :24.10  
 3rd Qu.:15.25   3rd Qu.:28.25  
 Max.   :20.00   Max.   :35.00

#descriptive  statistics
arrival_mean <- mean(arrivals_data$Patients)
arrival_mean

[1] 24.1

arrival_sd <- sd(arrivals_data$Patients)
arrival_sd

[1] 5.609203

arrival_var <- var(arrivals_data$Patients)
arrival_var

[1] 31.46316

Determine whether a Poisson distribution is appropriate.

##Interpretation:

The mean number of arrivals is 24.1 and the variance is approximately 31.46. For a Poisson distribution, the mean and variance should be approximately equal. Although the variance is somewhat larger than the mean, the values are reasonably close, so a Poisson distribution may be an appropriate model for the patient arrivals.

Estimate the arrival rate parameter.

lambda <- mean(arrivals_data$Patients)
lambda

[1] 24.1

##Interpretation:

The estimated arrival rate is 24.1 patients per hour, meaning the hospital receives approximately 24 patients per hour on average.

Calculate the probability of more than 30 arrivals in an hour.

proportion_arr <- mean(arrivals_data$Patients > 30)
proportion_arr

[1] 0.15

prob_for_more_than_30 <- ppois(30, lambda = lambda, lower.tail = FALSE )
prob_for_more_than_30

[1] 0.09952132

##Interpretation:

In the observed data, 15% of the hours had more than 30 patient arrivals. Using a Poisson distribution with an estimated arrival rate of 24.1 patients per hour, the probability of having more than 30 arrivals in an hour is approximately 10%. Therefore, according to the Poisson model, there is about a 10% chance that more than 30 patients will arrive in a given hour.

##Part B: Treatment Times

treatments_data

# A tibble: 15 × 2
   PatientID TreatmentMinutes
       <dbl>            <dbl>
 1         1               32
 2         2               27
 3         3               45
 4         4               19
 5         5               55
 6         6               40
 7         7               29
 8         8               38
 9         9               44
10        10               22
11        11               31
12        12               36
13        13               47
14        14               24
15        15               52

summary(treatments_data)

   PatientID    TreatmentMinutes
 Min.   : 1.0   Min.   :19.00   
 1st Qu.: 4.5   1st Qu.:28.00   
 Median : 8.0   Median :36.00   
 Mean   : 8.0   Mean   :36.07   
 3rd Qu.:11.5   3rd Qu.:44.50   
 Max.   :15.0   Max.   :55.00

Visualize the distribution.

hist(treatments_data$TreatmentMinutes,
     main = "Distribution of Treatment Times",
     xlab = "Treatment in Minutes",
     ylab = "Frequency")

##Interpretation:

The histogram shows that treatment times are concentrated around 30–40 minutes. The distribution is not perfectly bell-shaped, but it appears roughly symmetric with no obvious extreme outliers. Since the sample contains only 15 observations, it is difficult to determine normality from the histogram alone.

Estimate mean and standard deviation

# Treatment descriptive statistics
treatment_mean <- mean(treatments_data$TreatmentMinutes)
treatment_mean

[1] 36.06667

treatment_sd <- sd(treatments_data$TreatmentMinutes)

treatment_sd

[1] 11.02897

treatment_var <- var(treatments_data$TreatmentMinutes)
treatment_var

[1] 121.6381

Determine whether a Normal distribution is reasonable.

qqnorm(treatments_data$TreatmentMinutes,
       main = "Q-Q Plot of Patient Treatment Times")
qqline(treatments_data$TreatmentMinutes)

##Interpretation:

Based on the Q-Q plot, a Normal distribution appears to be reasonable for the treatment times. Most of the data points fall close to the reference line, with only small deviations at the tails. Therefore, the treatment times can reasonably be modeled using a Normal distribution.

Calculate the probability that treatment exceeds 45 minutes.

prob_that_trea_45 <- pnorm(45, mean = 36.06667, sd =  11.02897, lower.tail = FALSE )
prob_that_trea_45

[1] 0.2089736

##Interpretation:

Based on the Normal distribution model, the probability that a patient’s treatment time exceeds 45 minutes is approximately 20.9%. This means that about 21% of treatments are expected to take longer than 45 minutes. This is consistent with the histogram, where most treatment times are concentrated around 30–40 minutes.

#Part C: Simulation

Simulate one month of arrivals.

n = 720
# Simulation

set.seed(12345)

# Simulate g
simulated_arrival <- rpois(
  n = 720,

  lambda = 24.1
)

simulated_arrival

  [1] 26 27 23 21 35 26 22 19 21 25 26 20 11 22 29 25 27 31 20 15 32 21 24 17 28
 [26] 34 34 32 25 26 22 15 26 19 37 24 16 17 24 27 22 25 25 30 29 20 27 28 34 12
 [51] 24 26 28 19 15 26 29 19 24 20 32 24 22 24 22 27 23 25 25 28 24 19 33 27 22
 [76] 26 28 19 20 22 19 30 21 21 23 30 24 27 16 21 20 21 19 17 26 23 25 24 26 24
[101] 30 22 31 27 26 24 28 32 20 31 21 27 24 22 18 13 16 24 28 22 21 24 16 26 26
[126] 20 22 32 22 24 20 26 22 32 35 18 27 29 26 28 25 31 27 25 25 27 25 28 19 22
[151] 29 28 24 31 25 23 21 28 26 29 19 28 28 21 22 18 21 19 28 24 15 27 22 18 19
[176] 29 21 25 27 33 28 26 30 11 27 23 24 31 24 24 20 22 22 21 27 20 27 18 21 25
[201] 29 30 28 25 21 27 20 31 29 25 27 22 26 26 20 27 30 23 21 28 22 26 28 30 22
[226] 21 26 13 27 23 30 20 26 16 26 21 27 25 24 20 20 21 22 26 26 36 20 27 12 21
[251] 29 23 19 25 22 32 19 20 32 34 22 21 23 18 25 27 29 26 13 21 26 20 18 26 25
[276] 28 26 25 16 19 26 34 22 28 26 23 23 18 29 27 24 22 21 18 22 31 22 23 30 21
[301] 15 33 27 25 32 27 22 23 32 26 26 27 20 23 21 35 21 24 31 28 19 28 35 18 28
[326] 22 25 29 21 25 27 25 23 17 32 29 28 23 21 31 24 25 27 18 19 25 23 23 26 24
[351] 21 21 19 21 21 30 27 20 28 26 16 26 35 28 28 33 26 27 24 29 30 26 41 21 27
[376] 29 21 24 19 31 23 19 22 26 24 21 23 19 29 25 22 25 25 28 30 24 24 22 20 24
[401] 23 23 23 23 26 17 30 31 25 22 30 23 31 22 27 41 40 23 25 28 13 27 19 28 23
[426] 26 27 28 17 15 26 26 27 22 24 25 27 22 22 15 34 22 16 34 27 22 29 26 23 34
[451] 23 24 33 32 34 39 37 27 21 18 28 17 24 24 21 19 20 21 19 28 17 26 18 23 24
[476] 25 23 24 23 36 16 22 27 25 16 23 17 21 35 29 21 32 22 23 29 21 28 25 32 20
[501] 17 25 33 26 19 25 17 20 28 23 24 26 19 16 23 20 29 33 23 21 28 18 24 25 18
[526] 21 23 25 26 21 28 29 22 19 27 19 28 24 20 22 24 20 35 23 26 32 25 17 33 17
[551] 24 22 19 22 22 17 31 30 24 24 20 20 18 30 30 28 29 25 23 16 22 17 23 19 25
[576] 29 26 35 34 23 26 17 24 34 23 23 23 19 21 16 28 18 20 23 22 28 22 19 29 33
[601] 27 27 28 19 23 24 27 25 25 14 26 23 21 20 23 29 26 32 24 16 21 29 20 20 35
[626] 35 20 33 19 24 25 27 18 28 28 27 36 29 33 22 21 20 28 28 28 21 23 23 26 20
[651] 24 29 22 28 28 29 27 29 28 21 26 31 23 30 16 20 21 27 26 29 29 16 12 20 21
[676] 28 20 25 21 28 24 28 20 14 23 27 20 28 21 28 20 26 27 22 31 24 24 21 28 18
[701] 23 26 34 28 24 28 32 24 25 26 22 24 27 21 27 24 26 32 32 22

arrival_table <- data.frame(
  Hour = 1:720,
  Patients = simulated_arrival
)

head(arrival_table, 10)

   Hour Patients
1     1       26
2     2       27
3     3       23
4     4       21
5     5       35
6     6       26
7     7       22
8     8       19
9     9       21
10   10       25

Estimate expected daily demand.

# Daily arrivals
daily_arrivals <- numeric(30)

for (day in 1:30) {

  start_hour <- (day - 1) * 24 + 1
  end_hour <- day * 24

  day_arrival <- sum(simulated_arrival[start_hour:end_hour])

  daily_arrivals[day] <- day_arrival
}

daily_arrivals

 [1] 565 612 571 555 585 604 586 578 594 575 570 570 605 606 571 631 581 612 636
[20] 563 585 548 583 563 586 569 631 605 558 618

the_expected_daily_demand <- mean(daily_arrivals)
the_expected_daily_demand

[1] 587.2

# Theoretical expected daily demand
expected_daily <- lambda * 24

# Simulated average daily demand
simulated_daily <- mean(daily_arrivals)

expected_daily

[1] 578.4

simulated_daily

[1] 587.2

##Interpretation:

Based on the Poisson model, the expected daily patient demand is approximately 578 patients. The simulation estimates the average daily demand using 30 days of simulated arrivals. The simulated average should be close to the theoretical expected value, although small differences may occur due to random variation.

Recommend staffing improvements.

# Maximum daily arrivals
max(daily_arrivals)

[1] 636

# Minimum daily arrivals
min(daily_arrivals)

[1] 548

# Standard deviation of daily arrivals
sd(daily_arrivals)

[1] 23.89113

# Estimate daily treatment workload in hours
daily_treatment_hours <- daily_arrivals * treatment_mean / 60

# Descriptive statistics
summary(daily_treatment_hours)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  329.4   342.6   351.0   353.0   363.7   382.3

# Average daily treatment workload
mean(daily_treatment_hours)

[1] 352.9724

# Maximum daily treatment workload
max(daily_treatment_hours)

[1] 382.3067

##Interpretation:

Based on the simulation results, the hospital is expected to receive approximately 578 patients per day. However, daily arrivals vary, with a minimum of 548 patients and a maximum of 636 patients. The standard deviation of approximately 23.89 patients shows that demand changes from day to day.

Based on these results, the hospital could improve staffing by:

Planning staffing around the expected daily demand while considering possible fluctuations.
Maintaining additional staff capacity during high-demand periods to reduce patient waiting times.
Considering average treatment duration when estimating staffing needs, rather than relying only on patient arrival counts.
Monitoring hourly arrival patterns to identify peak periods and adjust staff schedules accordingly.

The estimated daily treatment workload can also help the hospital understand how much treatment time may be required. However, additional information about staff availability, patient severity, and treatment resources would be needed to determine the exact number of staff required.

Other Formats

Homework 3