"timer"<-
function(myfun)
{
	tstart <- seconds()
	eval(parse(text = myfun))
	tend <- seconds()
	return(tend - tstart)
}
"seconds"<-
function()
{
	d <- stamp(print = F)
	hr <- as.numeric(substring(d, 12, 13))
	min <- as.numeric(substring(d, 15, 16))
	sec <- as.numeric(substring(d, 18, 19))
	return(3600 * hr + 60 * min + sec)
}
"ideal.boot"<-
function(x, theta = mean, ...)
{
#######################################################################
# Computes the ideal bootstrap se of statistic theta(x), where
# x is the observed data. Generates partitions of n recursively and 
# uses formulae 6.8, 6.16 & 6.17 of Efron and Tibshirani (1992).
# This algorithm provides a general solution to problem 6.10.
#
# Uses auxiliary functions:
#  generate.partition : recursively generates a partition
#  bootprob           : probability of a partition
#
	n <- length(x)
	j <- numeric(n)
	X1 <- X2 <- 0
	while(!is.null(j <- generate.partition(j))) {
		thetastar <- theta(x[rep(1:n, j)], ...)
		probability <- bootprob(j)
		X1 <- X1 + thetastar * probability
		X2 <- X2 + thetastar^2 * probability
	}
	return(ideal.se = sqrt(X2 - X1^2), ideal.mean = X1)
}
"trimmean"<-
function(x)
{
	mean(x, trim = 0.25)
}
"prob610.data"<-
c(1.2, 3.5, 4.7, 7.3, 8.6, 12.4, 13.800000000000002, 18.1)
"generate.partition"<-
function(j)
{
###################################################################################
#generate all non-negative integer solutions to the equation
# j_1+j_2+...+j_n = n
#Algorithm:
#Step 1. Start with j == (n,0,...,0) = (j_1, ..., j_n)
#Step 2. Find the largest index, i, such that j_i + j_(i+1) + ... + j_n < n.
#Step 3. Increment j_i and then set everything to its left to 0 except for
#        j_1 which is set to the largest possible value.
#Step 4. Repeat Steps 2 & 3 until (0,...,n) is reached. This occurs when
#        there is no value of j_i from Step 2 which can be found.
#
	n <- sum(j)
	if(n == 0)
		return(c(length(j), rep(0, length(j) - 1)))
	ji <- sum(cumsum(rev(j)) < n)
	if(!ji)
		return(NULL)
	else ji <- (n:1)[ji]
	if((ji == 2)) {
		j[ji - 1] <- j[ji - 1] - 1
	}
	else {
		j[1:(ji - 1)] <- 0
		j[1] <- n - 1 - sum(j[ji:n])
	}
	j[ji] <- j[ji] + 1
	return(j)
}
"bootprob"<-
function(j)
{
############################################################################
#computes the probability of a bootstrap sample of a vector x, where there
#are j[1] of x[1], j[2] of x[2], ..., j[n] of x[n], 
#where n=length(x).
#
#Reference: Efron and Tibshirani (1993, formulae 6.8, 6.16, 6.17)
#
#use logs and then exponeniate
#
	N <- sum(j)
	exp(lgamma(1 + N) - sum(lgamma(j + 1)) - N * log(N))
}
prob610time_timer(ideal.boot(prob610.data, theta=trimmean))
print(prob610time)