[R] Biased calculations from using rmultinom?

[Sebastian S]

Dear R users,

This is a problem that has puzzled me for quite some time, and I wonder if 
anyone could offer me an insight to what is happening?

I generated a multinomial matrix by using 
A<-rmultinom(100,1,c(0.4,0.3,0.2,0.1), and multiply this matrix by a set of 
values say B<-runif(n,100,50). I then calculate var(colSums(A*B)). I was 
able to work out the theoretical mean of var(colSums(A*B)) by using 
fun.var.match.gen below.

I was however not sure, why is there a persistent over estimate from my 
fun.var.match.gen? I have enclosed the tests I have used to demonstrate this 
fact. The difference is very small but very persistent throughout 
simulations. Is there something I have overlooked here?

All comments very welcome!

Sebastian.

#### Functions used:

fun.split.matrix<-
function(object, index)
{
	result <- lapply(split(object, index), function(x, object)
	matrix(x, ncol = ncol(object)), object)
	return(result)
}

fun.var.match.gen<-function(q,n,mean.cost.m,var.cost.m){
	k<-length(q)
	mean.o<-sum(n*q*mean.cost.m)/k

	result<-sum(n*var.cost.m*q+mean.cost.m^2*(n*q*(1-q+n*q))-2*mean.o*n*q*mean.cost.m+mean.o^2)/(k-1)
}


# Test:

p1<-c(10,10,1,8,8,6,10,9,6,5)
p1<-p1/sum(p1)
e<-5
f<-60000

n1<-4
n.rep<-50000
n.gen<-n1*n.rep

A<-fun.runif(n.gen,e,f)

temp1a<-rmultinom(n.gen,1,p1)*A
# Randomly select items to split
sam<-sample(rep(1:n.rep,each=n1))
temp3<-sapply(fun.split.matrix(temp1a,sam),colSums)

junk.vars<-apply(temp3,2,var)
junk.m<-mean(junk.vars)
junk.theo<-fun.var.match.gen(p1,n1,(e+f)/2,(f-e)^2/12)

plot(junk.vars)
abline(h=junk.theo,col=2)
abline(h=junk.m)


http://adsfac.net/link.asp?cc=FXT002.7542.0