splstr <- function(mu=20,A=c(4,0,-4),B=c(2,-2),C=c(3,-3),
	AB=matrix(0,3,2),AC=matrix(c(2,-2,0,-2,2,0),3,2),BC=matrix(c(1,-1,-1,1),2,2),
	ABC=array(0,c(3,2,2)),
	n=c(3,2,2,2),sigma=c(3,2,1,1))
{
	y <- array(mu,n)
	varcalc <- function(n,sd,mean=0)
	{
		array(rnorm(n=prod(n),mean=mean,sd=sd),n)
	}
	varS <- varcalc(n[1:2],sigma[1])
	varR <- varcalc(n[1:3],sigma[2])
	varC <- varcalc(n[c(1:2,4)],sigma[3])
	var <- varcalc(n,sigma[4])
	for (i in 1:n[1])
	for (j in 1:n[2])
	for (k in 1:n[3])
	for (m in 1:n[4])
		y[i,j,k,m] <- mu + A[i] + varS[i,j] +
			B[k] + AB[i,k] + varR[i,j,k] +
			C[m] + AC[i,m] + varC[i,j,m] +
			BC[k,m] + ABC[i,k,m] + var[i,j,k,m]
	cbind(A=array(1:n[1],prod(n)),
		S=array(1:prod(n[1:2]),prod(n)),
		B=array(matrix(1:n[3],prod(n[1:2]),n[3],byrow=T),prod(n)),
		C=c(matrix(1:n[4],prod(n[1:3]),n[4],byrow=T)),
		Y=y)
}