Demo entry 6782356

Central Limit Theorem

   

Submitted by Peter Saeta on Jan 13, 2019 at 19:55
Language: Igor. Code size: 2.1 kB.

// An example function that takes a variable in the range (-1, 1)
// and returns a real value. It will be called with random numbers
// uniformly distributed in the range (-1, 1).
Function nop(x)
	Variable x									// parameters are declared in Igor
	return x
End

// SampleN takes a function of a single variable, such as nop, and
// a number of times to sample the function and add the results. It
// does this one-hundred-thousand times, putting the results in the
// wave "accumulator". It then prepares a histogram of the values in
// accumulator, and plots the values, along with uncertainties
// given by the square root of the number in each bin.
Function SampleN(f, ntimes)
	FUNCREF nop f
	Variable ntimes
	
	Make/O/N=100000 accumulator				// create accumulator with 1e5 points
	accumulator = 0								// set all values to zero
	Variable n
	for (n = ntimes ; n > 0; n -= 1)		// loop ntimes times
		accumulator += f(enoise(1))			// add a new value to each element
	endfor
	
	// Prepare the histogram
	
	Make/N=51/O myhisto						// create myhisto to hold histogram
	// Compute the histogram from accumulator, putting the results
	// in myhisto, centering the x-coordinate of each bin (/C flag)
	// and computing a wave W_SqrtN holding the square root of the
	// number of counts in each bin (/N flag). The width of the bins
	// is set using an algorithm described in the documentation (/B=4).
	Histogram /B=4 /C/DEST=myhisto/N accumulator
	WAVE W_SqrtN									// declare wave created by Histogram
	DoWindow/F MyGraph							// see if we have already made a graph
	if (V_flag == 0)							// if not,
		Display/K=1 myhisto					// graph myhisto and then
		DoWindow/C MyGraph						// rename the window so we can find
		HMC#FixGraph()							// fix the graph, and add error bars
		ErrorBars myhisto Y,wave=(W_SqrtN,W_SqrtN)
		SetAxis left 0,*						// make sure y axis starts at 0
	endif
	
	Execute/P "CurveFit/G/H=\"1000\"/TBOX=784 gauss myhisto /W=W_SqrtN /I=1 /D /R"
	Execute/P "HMC#AddChiSqInfoToPlot()"
End

This snippet took 0.00 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).