/* Chi-square Random variable*/
	/* generazione di 1000 osservazioni (c=100 samples di size n=10) 
      da una normale di media 3 e varianza 4*/
	data chivar;
	 keep x sample c scale;
	  scale=2; mx=3; n=10; c=100; nsim=n*c; *n=sample size c=#samples;
	   do i=1 to nsim;
	      x=rannor(1)*scale+mx;
		  sample=int((i-1)/n)+1; 	/* la funzione int(x) restituisce la parte intera di x */
	      output;
	   end;
	 run;

	 	/* Calcolo della devianza campionaria */
	proc means data=chivar noprint;
	 var x scale;
	 by sample;
	   output out=chivarmean css=cssx mean=meanx scale; 
	run; 
	proc print data=chivarmean; run;

	/* calcolo della statistica Q=(n-1)(S^2)/sigma^2 */
	data x2;
	 set chivarmean;
	 keep sample q;
	 q=cssx/(scale**2);  
	run;

	proc means data=x2 mean var; var q; run;  /* E(X^2)=g d.f. Var(X^2)=2g */

		/* Plot della distribuzione chi-quadrato teorica ed empirica
	   chi-quadrato con g=(n-1) df è una Gamma(g/2, 2)*/

	proc univariate data=x2 noprint;
	   histogram q / gamma(percents=2.5 20 40 60 80 95 97.5 midpercents alpha=4.5 sigma=2 ) endpoints = 0 to 30 by 5; *alpha=0.5 sigma=sqrt(18);
	   inset n='# di campioni' mean='mean=d.f.' (5.2) var='var=2d.f.' gamma(ksdpval);  * ksdpval: kolmogorov-smirnov edf test pvalue;
	    qqplot q  / gamma(alpha=4.5 sigma=2 ) square ctext=blue;
	run;

	proc format;
	value classiq 0-4='<4'
	             4-12='4-12'
				 12-high='>12'; 
	run;
	proc freq data=x2;
	format q classiq.;
	tables q;
	run;