outer_sum   = 0.0
outer_sum2  = 0.0

for outeri in range(0, nouter):
	inner_sum = 0.0
	for inneri in range(0, ninner):
		X = random.uniform(0, 1)
		inner_sum   += X
	inner_avg   = inner_sum / ninner
	outer_sum  += inner_avg
	outer_sum2 += inner_avg**2

# See Flegal and Jones 2009 paper.
outer_avg = outer_sum / nouter
outer_var = ninner * (outer_sum2 - outer_sum**2 / nouter) / (nouter - 1);
outer_stddev = sqrt(outer_var)