# Find the script "urn.rb" in the current Ruby library search path
# and execute it.  Note that the current working directory (i.e., the
# directory from which you started ruby) is in the Ruby library search
# path.  This makes available the Urn class.
require 'urn'

# Level of significance
alpha = 0.05
theta = 0.50
n_balls  = 10
n_draws = 100

# Draw the sample by calling the "sample" class method.
draws = Urn.sample(n_balls,theta,n_draws,10000)

# Observed test statistic
observed_test_statistic = 0.6

# Find the script "dbd/broccoli.rb" in the current Ruby library search path
# and execute it.  This makes available several modules used below.
require 'dbd/broccoli'

# Get the Monte Carlo integral for "h(X) = X", i.e., get the mean
mean_under_Ho = DBD::Broccoli::Integration::MonteCarlo::integrate(draws) do |x|
  x
end

# Get the Monte Carlo integral for "h(X) = (X-E(X))^2", i.e., get the variance
variance_under_Ho = DBD::Broccoli::Integration::MonteCarlo::integrate(draws) do |x|
  (x-mean_under_Ho)**2
end

# Having to type "DBD::Broccoli::Integration::MonteCarlo" is a pain, so
# let's include this module into the current environment.
include DBD::Broccoli::Integration::MonteCarlo

# Now there is less to type.  Let's do it again:
mean_under_Ho = integrate(draws) do |x|
  x
end

variance_under_Ho = integrate(draws) do |x|
  (x-mean_under_Ho)**2
end

# Blocks be delimited by 'do' and 'end' (as above) or simply by '{' and '}'
mean_under_Ho = integrate(draws) { |x| x }
variance_under_Ho = integrate(draws) { |x| (x-mean_under_Ho)**2 }

# It turns out that the module "DBD::Broccoli::Integration::MonteCarlo" already
# has the mean and variance defined:
mean_under_Ho = mean(draws)
variance_under_Ho = variance(draws)

# Use Monte Carlo integration to estimate the 1-alpha quantile of the null dist.
# Again, we can explicitly refer to the methods of interest.
critical_value = DBD::Broccoli::Integration::MonteCarlo::quantile(1-alpha,draws)
alpha = 1 - DBD::Broccoli::Integration::MonteCarlo::probability(critical_value,draws)
p_value = 1 - DBD::Broccoli::Integration::MonteCarlo::probability(observed_test_statistic,draws)

# Or, we can take advantage of the fact that we included the module earlier.
critical_value = quantile(1-alpha,draws)
alpha = 1 - probability(critical_value,draws)
p_value = 1 - probability(observed_test_statistic,draws)

# Likewise, let's include DBD::Broccoli::Tools, which provides the method "seq"
# (among other things).  By the way, "require" and "include" are typically placed
# at the top of a chuck of code, but you don't have to do this.
include DBD::Broccoli::Tools

theta = seq(0.35,1.00,0.05)
power = Array.new(theta.length)
for i in 0...(theta.length)
  draws = Urn.sample(n_balls,theta[i],n_draws,5000)
  power[i] = 1 - probability(critical_value,draws)
end

# Make R available to use from within Ruby.  Access it via the $R global variable.
require 'dbd/R'

# Convert the Ruby arrays "theta" & "power" into a R vectors having the same name
$R.assign("theta",theta)
$R.assign("power",power)

# Or, more conveniently, we can just do the following:
$R.theta = theta
$R.power = power

# Let's see what version of R we've got running.
# Note that we cannot use the convenience feature above since '$R.R.version.string'
# would not be syntactically correct in Ruby.
R_version_long = $R.pull("R.version.string")

# Just get the version number
R_version_short = R_version_long.split[2]

# Evaluate the R code defined in this "here document".  See page 62 of the pickaxe book.
# Note that the value of the Ruby variance 'alpha' gets substituted because passing
# the commands to R.
$R.eval <<EOF
  plot(theta,power,type="l",ylab="Power of the Test",
       xlab=expression(paste(theta," (True Proportion of Black Balls)")),
       main="Power for Test of Ho: Proportion = 0.5")
  points(0.5,#{alpha},pch=16)
EOF

print "Press <enter> to continue."
$stdout.flush
gets