Chisquared Distribution Plot (xfig example)

The Chi²-distribution and its according cumulative distribution function for a large number of degrees of freedom for illustrative purposes. Confidence levels are indicated in the CDF. Feel free to use it in your thesis, there are some customization options in the header.

#!/usr/bin/env isis-script

require("isisscripts");

%---------options-----------------------------
variable CHISQRMIN = 0.0003, 
	 CHISQRMAX = 150,
	 KVALUES = [1,2,4,8,15,30,60,100],
	 COLOR = ["blue","crimson"],
	 LOG = 1,
	 dx=0.001;


%---------we-need-the-gamma-function------------------
%  defined for integer and half-integer values
define gamma_function(n){
	variable gamma = 0;

	if(n mod 1 == 0){
		gamma = factorial(n-1);
	}else if(n mod 1 == 0.5){
		n -= 0.5;
		gamma = sqrt(PI)*factorial(2*n)/(4^n*factorial(n));
	}
	
	return gamma;
}


%-------initializing-panels-----------------------
variable panel = Struct_Type[2];

panel[0] = xfig_plot_new(12, 9);
panel[0].world( CHISQRMIN+LOG*0.1, 0.999*CHISQRMAX, 0, 0.5; xlog=LOG);
panel[0].ylabel("$f_k (\chi^2)$"R);

panel[1] = xfig_plot_new(12, 6);
panel[1].world( CHISQRMIN+LOG*0.1, 0.999*CHISQRMAX, 0, 1.1; xlog=LOG);
panel[1].ylabel("$F_k (\chi^2)$"R);

variable chisqr = [CHISQRMIN:CHISQRMAX:dx];
variable prob_density = Double_Type[length(chisqr)];
variable cumulative_distr = Double_Type[length(chisqr)];


%-------calculate-and-plot-distribution-for-k-values------------------------------
% (k larger than 100 not possible as factorial breaks down)
variable i, k, labely = 0.9, labelx = 0.8, ystep = 0.06, sigma_flag;

foreach k (KVALUES){
	variable kcolor = xfig_mix_colors(COLOR[0], COLOR[1], log(k-1)/5 );
	sigma_flag = 0;

	% calculating probability densities and integrating them in the same loop
	for(i=0; i<length(chisqr)-1; i++){
		prob_density[i] = ( chisqr[i]^(k/2.-1.)*E^(-chisqr[i]/2.) )/( 2.^(k/2.)*gamma_function(k/2.) );
		cumulative_distr[i] = cumulative_distr[i-1]+prob_density[i]*dx;

		% plot a line at edge of confidence level
		if(cumulative_distr[i] >=0.9 && sigma_flag == 0){
			panel[0].plot([chisqr[i],chisqr[i]], [0,prob_density[i]]; line=2, color = kcolor, depth=100-k);
			panel[1].plot([chisqr[i],chisqr[i]], [0.9,1.1]; line=2, color = kcolor, depth=100-k);
			sigma_flag = 1;
		}
	}
	
	% plot and label the probability density for k
	panel[0].plot(chisqr, prob_density; line=0, color = kcolor, depth=100-k);
	panel[0].xylabel(labelx, labely , "k$\,=\,$"R + string(k), -0.5, 0; world0, color = kcolor);
	labely -= ystep;

	% plot according cumulative distribution function
	panel[1].plot(chisqr, cumulative_distr; line=0, color = kcolor, depth=100-k);
}

% dashed line in second panel
panel[1].plot([0.1,200], [1,1]; line=1, color = "gray", depth=100);
panel[1].plot([0.1,200], [0.9,0.9]; line=1, color = "gray", depth=100);
panel[1].xylabel(0.05, 0.77, "90\%-level"R, -0.5, 0; world0, color = "gray");


%-------render-all---------------------------------------
variable multi = xfig_multiplot(panel; xlabel=`$\chi^2$`, cols=1);
multi.render("chisqr.pdf");

quit;