# Author: Adolfo J. Rumbos
# Date: September 5, 2008

# Randomization Test for Westvato data (continued)
#
# Defining pHat()
#
# This script is to be run after the Simulation2roundHourlyWorkers.R has been 
# run.  It assumes that the Xbar vector has been obtained.  
# This program defines a function, pHat(), which 
# computes the proportion, p_Hat, of samples in the simulation which yield a 
# mean greater than or equal to an observed value, obs. This yields 
# an estimate of the probability that, under the assumption that Westvaco's 
# choice of the three hourly workers to be laid off was entirely random 
# (i.e., each worker, regardless of her or his age, had an equal likelihood
# of being chosen).
# 
# The function takes on the two arguments: X and obs, where 
# X is a vector that contains values for a test statistic obtined from
# many repetions of a sampling, and obs is the  value of the statistic
# observed in the data.
# The result is the proportion of samples that yield a test statistic 
# equal to, or larger than, the observed value
#
#
pHat <- function(X,obs)   # Sets up to define the function Phat
{

L<-length(X) 		# Computes length of X
NYes <- 0			# Sets up counter of "yeses"

for (i in 1:L)		# Sets up loop through vector X
{
NYes <- NYes + (X[i]>=obs) # adds 1 to counter if Xbar[i]>= obs
}
result <- NYes/L 		# Computes proportion of NYes in L
return(result)		# Returns estimate of p_hat