Difference between revisions of "Chisquared Distribution Plot (xfig example)"

From Remeis-Wiki
Jump to navigation Jump to search
Line 2: Line 2:
  
 
[[Category:SLxfig]]
 
[[Category:SLxfig]]
[[File:chisq.png|560px|]]
+
[[File:chisq.png]]
  
 
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.
 
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.

Revision as of 13:43, 13 April 2018

File:Chisq.png

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;