function(n, x, idcol, ntrial)
# n is the sample size, x is the size vector (it doesn't have to sum
# to one, idcol is the vector sampled--usually a set of labels--and
# ntrial is the number of simulations.  A list is output--it contains
# idcol as its first element, the vector of pi_i's as its second
# element and a matrix of pi_ij's (with the pi_i's along its diagonal)
# as its third element--look at accompanying text to see how to use
# the output to compute estimators and standard errors.
{
	px <- x/sum(x)
	tally.piij <- matrix(rep(0, length(x) * length(x)), ncol = length(x))
	id <- seq(1, length(x))
	for(i in 1:ntrial) {
		samp.pi <- sample(id, size = n, prob = px)
		tally.piij[samp.pi, samp.pi] <- tally.piij[samp.pi, samp.pi] + 
			1
	}
	tally.piij <- tally.piij/ntrial
	tally.pi <- diag(tally.piij)
	out <- list(idcol, pihat = tally.pi, piijhat = tally.piij)
	out
}