Program 3 - Histogram

Program Code

Off[General::"spell1",General::"spell"]

BeginPackage["Hystogram`"];

Unprotect[Hystogram];

Hystogram::usage = "Hystogram[f, x0, numtokeep, numtodrop, numbins] creates a plot showing the 
density of the state variable x for a function f(x).";

Clear[Hystogram, Density, Fraction, Counts];

Options[Hystogram] = {
	PlotStyle -> Hue[2/3],
	Scale -> Density,
	PlotRange -> All
};

Begin["`Private`"];

RuleQ[r_] := Head[r] === Rule;
RealQ[r_] := NumericQ[r] && Im[r] == 0;

Hystogram::badfunc = "`2` detected in function `1`.";
Hystogram::badopt = "Ignoring unknown `2`option `1` in argument list.";
Hystogram::badpr = "Ignoring invalid PlotRange command `1`.";
Hystogram::badpt = "Initial point `1` is not a real number.";
Hystogram::badit1 = "Number of iterations `1` is not a positive integer.";
Hystogram::badit2 = "Number of iterations to drop `1` is not a non-negative integer.";
Hystogram::badbin = "Number of bins `1` is not a positive integer greater than 1.";

Hystogram [f_, x0_, numiterations_, numtodrop_, numbins_, args___]:=
		Module[{g, var, x, points, maxval, minval, dx1, dx2, plotopts, baropts, type, gropts, pr, checklist},
	Check[f[x0], Return[]];
	Switch[Head[f],
		Function, g = f,
		Symbol, If[MemberQ[Attributes[f], NumericFunction], g = f, g = #&],
		_, var = DeleteCases[Union[Flatten[Apply[List, {f}, {0, Infinity}]]], _?RealQ];
			Switch[Length[var],
				0, Message[Hystogram::badfunc, f, "No variables"]; Return[],
				1, If[Head[var[[1]]] === Symbol,
					g[x_] := f /. var[[1]]->x,
					Message[Hystogram::badfunc, f, "Unknown symbol"]; Return[]
				],
				_Integer, Message[Hystogram::badfunc, f, "Too many variables"]; Return[]
			];
	];
	
	If[!RealQ[x0], Message[Hystogram::badpt, x0]; Return[]];
	If[!IntegerQ[numiterations] || numiterations < 1,
		Message[Hystogram::badit1, numiterations]; Return[]];
	If[!IntegerQ[numtodrop] || numtodrop < 0,
		Message[Hystogram::badit2, numtodrop]; Return[]];
	If[!IntegerQ[numbins] || numbins < 2,
		Message[Hystogram::badbin, numbins]; Return[]];
	
	plotopts = {args};
	
	checklist[l_, test_, mess_, str_] := Module[{i, results = l},
		For[i = 1, i <= Length[results], Null,
			If[!test[results[[i]]],
				Message[mess, results[[i]], str];
				results = Delete[results, i],
				i++
			]
		];
		Return[results]
	];
	SetAttributes[checklist, HoldAll];
	plotopts = checklist[plotopts, RuleQ, Hystogram::badopt, ""];
	
	{pr, baropts, type} = {PlotRange, PlotStyle, Scale} /. plotopts /. Options[Hystogram];
	plotopts = DeleteCases[plotopts, (PlotRange -> _) | (PlotStyle -> _) | (Scale -> _)];
	
	gropts = #[[1]]& /@ Options[Graphics];
	plotopts = checklist[plotopts, MemberQ[gropts, #1[[1]]]&, Hystogram::badopt, "graphics "];
	
	x = Nest[g, x0, numtodrop];
	points = NestList[g, x, numiterations];
	
	If[MemberQ[{All, Automatic}, pr] || VectorQ[pr, RealQ] && Length[pr] == 2,
		pr = {All, pr},
		If[!VectorQ[pr, (MemberQ[{All, Automatic}, #] || VectorQ[#, RealQ] && Length[#] == 2)&] || Length[pr] != 2,
			Message[Hystogram::badpr,pr];
			pr = {All, All}
		];
	];
	
	{minval, maxval} = {Min[points], Max[points]};
	dx1 = (maxval - minval) / numbins;
	If[!MemberQ[{All, Automatic}, pr[[1]]],
		{minval, maxval} = pr[[1]]
	];
	dx2 = (maxval - minval) / numbins;
	
	bins = Compile[{{pts, _Real, 1}, {xmin, _Real}, {w, _Real}, {nb, _Integer}, {ni, _Integer}},
		Module[{results, i, j},
			results = Table[0, {nb}];
			For[i = 1, i <= ni, i++, 
				j = Round[(pts[[i]] - xmin) / w] + 1;
				If[j >= 1 && j <= numbins, results[[j]]++];
			];
			results
		]
	] [points, minval, (maxval - minval) / (numbins - 1), numbins, numiterations];

	Switch[type,
		Counts, Null,
		Fraction, bins = bins / numiterations,
		_, bins = bins / (numiterations * dx1)
	];
	
	bars = Table[Rectangle[{minval + (i-1) dx2, 0}, {minval + i dx2, bins[[i]]}], {i, 1, numbins}];
	
	Return[Show[Graphics[Join[{baropts}, bars]], plotopts, Frame -> True, PlotRange -> pr]];
];

End[];

SetAttributes[Hystogram, ReadProtected];
Protect[Hystogram];

EndPackage[];

On[General::"spell1",General::"spell"];

Examples

Hystogram[2  x^2 - 1, .3, 2000, 100, 100];

Hystogram[3.3 x (1 - x), 0.2, 1000, 0, 100, 
      Scale -> Counts];