n.trials = 20000

n.hits.in.4=0
n.hits.in.24=0

for (trial in 1:n.trials) {
	
	# 4 dice
	
	y1to4 = sample(1:6,4,replace=TRUE)
	hitIn4 = any(y1to4 == 6)
	roll4 = paste(y1to4,collapse="")
	
	n.hits.in.4 = n.hits.in.4 + hitIn4
	
	# 24 double dice
	
	y1to24 = matrix(sample(1:6,48,replace=TRUE),24,2)
	y1to24 = apply(y1to24,1, function(x) x[1]*10+x[2] )
	hitIn24 = any(y1to24 == 66)
	
	n.hits.in.24 = n.hits.in.24 + hitIn24
	
	roll24pairs=paste( y1to24, collapse=" ")
	
	print( noquote(
	  paste(trial,
	        " ", roll4, " (",n.hits.in.4,")  ",
	             roll24pairs, " (", n.hits.in.24,")",sep=""
	        )
	  ) 
	)
}


# no. of trials needed to measure expected difference
# in proportion of hits with 

# (95%) ME of .2 of a percentage point  = 0.002

# SE of difference in 2 proportions both close to .5 

# sqrt( .5*.5 / n  + .5*.5 / n) = sqrt( 1/(2*n))

# 95% ME = 2 * sqrt( 1/(2*n))

# ME^2 = 4 / (2*n) = 2/n

# n per game = 2 / ME^2   = 2 / (0.002^2) = 2 * 500^2 = 2 * 250000 = 500,000

# if ME = 1 percentage point = 0.01

#  n per game = 2 / ME^2   = 2 / (0.01^2) = 2 * 100^2 = 2 * 10000 = 20,000