# reading the data

ground = c(0,1,2,3,4)
aircraft=c(0,1,2,3,4)

# simple bootstrap to estimate means, se and confidence intervals

ms1=numeric(1000000)
for(i in 1:1000000){
  ms1[i]=mean(sample(ground,replace=TRUE))}
ms1s=sort(ms1)
cat("mean ground = ",mean(ms1), "SE = ",sd(ms1),"95% PU = ",ms1s[25001], ", ", ms1s[975000],"\n")

par(mfrow=c(1,1))
hist(ms1,main=NULL,xlab="ground",ylab="frequency")

ms2=numeric(1000000)
for(i in 1:1000000){
  ms2[i]=mean(sample(aircraft,replace=TRUE))}
ms2s=sort(ms2)
cat("mean aircraft = ",mean(ms2), "SE = ",sd(ms2),"95% PU = ",ms2s[25001], ", ", ms2s[975000],"\n")

hist(ms2,main=NULL,xlab="aircraft",ylab="frequency")

# PERMUTATION TEST BASED ON RESAMPLING WITHOUT REPLACEMENT
# Combine the two datasets into a single dataset
# i.e., under the null hypothesis, there is no difference between the two groups
combined = c(ground,aircraft)

# Observed difference
diff.observed = mean(ground) - mean(aircraft)

number_of_permutations = 100000

diff.random = NULL
for (i in 1 : number_of_permutations) {
  
  # Sample from the combined dataset without replacement
  shuffled = sample (combined, length(combined))
  
  ground.random = shuffled[1 : length(ground)]
  aircraft.random = shuffled[(length(aircraft) + 1) : length(combined)]
  
  # Null (permutated) difference
  diff.random[i] = mean(ground.random) - mean(aircraft.random)
}

# P-value is the fraction of how many times the permuted difference is equal or more extreme than the observed difference

pvalue = sum(abs(diff.random) >= abs(diff.observed)) / number_of_permutations
cat ("P value from the permutation test: ",pvalue)