Central Limit Theorem Code

This code computes the sample mean, Ȳ, of two independent identical distributions, Y_i

This is done repeatedly, using the Ȳ of the previous iteration as Y_i in the next one, to efficiently compute Ȳ for very large numbers of the original Y_i

(The graphics drawing functions are omited for the sake of brevity, but the comments should be descriptive enough.)

function getstats(distr) {
	var scale = 0;
	var mean = 0, meansq = 0;
	for (var c = 0; c < distr.length; c++) {
		var y = c / (distr.length-1); // y goes from 0 to 1.
		mean += y*distr[c]; // y * p(y)
		meansq += y*y*distr[c]; // y² * p(y)
		scale += distr[c];
	}
	mean /= scale;   // E[Y]
	meansq /= scale; // E[Y²]
	var stdev = Math.sqrt(meansq - mean*mean); // sqrt(E[Y²] - E[Y]²)
	return "&mu; = " + +mean.toPrecision(3) + ", &sigma; = " + +stdev.toPrecision(3);
}

function computeYbars(Yi) {
	output_box.innerHTML = ""; // Reset the output area on the page.
	stats_box.innerHTML = getstats(Yi); // Calculate the mean and standard deviation of the original Yi distribution.
	
	var n = 2;
	for (var c = 0; c <= 10; c++, n <<= 1) { // Double n 10 times, to generate Y bar in powers of 2, up to n = 2048.
		var Ybar = [];
		
		var half = Yi.length >> 1;
		
		// Average the overlap of the two distributions to get Y bar,
		// and average the values into the original discrete bins
		// to keep the number of discrete values from exploding.
		var max = 0;
		for (var y = 0; y < Yi.length; y++) { // Loop through every possible y value.
			var lower = y < half ? 0 : Yi.length - (Yi.length - y)*2 + 1; // Two cases.
			var upper = y < half ? y*2 : Yi.length-1;
			
			Ybar[y] = 0;
			for (var k = lower; k <= upper; k++) // Add the probability of every possible sum of (Y1 + Y2)/2 that has this y value. 
				Ybar[y] += Yi[k]*Yi[(upper - k) + lower];
			
			if (Ybar[y] > max) max = Ybar[y];
		}
		
		// Scale to fit the screen.
		var scale = (main_canvas.height * 0.9) / max;
		for (var y = 0; y < Ybar.length; y++)
			Ybar[y] *= scale;
		
		// Draw it, and compute the mean and standard deviation.
		drawNewDistribution(n, Ybar, getstats(Ybar));
		
		Yi = Ybar; // Use this Ybar as Yi for the next round.
	}
}