// Chi-square, df=1 to 3
twoway (function y=chi2den(1,x), range(0 4)) ///
	(function y=chi2den(2,x), range(0 4)) ///
	(function y=chi2den(3,x), range(0 4)), title("Chisquare") ///
	legend(order(1 "df=1" 2 "df=2" 3 "df=3"))

clear
set obs 1000
forvalues i=1/3 {
	generate x`i' = rchi2(`i')
	}
// Mean = df; SD = sqrt(2*df), so not all equal
summarize x*


// Chi-square, non-integer df
twoway (function y=chi2den(3,x), range(0 4)) ///
	(function y=chi2den(3.5,x), range(0 4)) ///
	(function y=chi2den(4,x), range(0 4)), title("Chisquare, non-integer") ///
	legend(order(1 "df=3" 2 "df=3.5" 3 "df=4")) name(chi2)

// aka Gamma, non-integer shape and scale
twoway (function y=gammaden(3/2,2,0,x), range(0 4)) ///
	(function y=gammaden(3.5/2,2,0,x), range(0 4)) ///
	(function y=gammaden(4/2,2,0,x), range(0 4)), title("Gamma, non-integer") ///
	legend(order(1 "df=3" 2 "df=3.5" 3 "df=4")) name(gamma)

clear
set obs 1000
forvalues i=1/3 {
	local j = 3+(`i'-1)/2
	generate x`i' = rchi2(`j')
	}
// Mean = df; SD = sqrt(2*df), so not all equal
summarize x*

// Chi-square, df=3, non-central
twoway (function y=chi2den(3,x), range(0 4)) ///
	(function y=nchi2den(3,3*0.5^2,x), range(0 4)) ///
	(function y=nchi2den(3,3,x), range(0 4)), title("Chisquare, non-central") ///
	legend(order(1 "ncp=0" 2 "ncp=0.75" 3 "ncp=3")) name(nchi1)

clear
set obs 1000
forvalues i=1/3 {
	local j = (`i'-1)/2
	forvalues k=1/3 {
		generate y`i'`k' = rnormal(`j', 1)^2
		}
	egen x`i' = rowtotal(y`i'*)
	}
// Non-cenrality parameter, ncp, is sum(`j'^2)
// Mean = df+ncp; SD = sqrt(2*(df+2*ncp)), so not all equal
summarize x*

// Chi-square, df=3, non-central
twoway (function y=chi2den(3,x), range(0 4)) ///
	(function y=nchi2den(3,0.5,x), range(0 4)) ///
	(function y=nchi2den(3,1,x), range(0 4)), title("Chisquare, non-central") ///
	legend(order(1 "ncp=0" 2 "ncp=0.5" 3 "ncp=1")) name(nchi2)
	
clear
set obs 1000
forvalues i=1/3 {
	local j = sqrt((`i'-1)/(2*3))
	forvalues k=1/3 {
		generate y`i'`k' = rnormal(`j', 1)^2
		}
	egen x`i' = rowtotal(y`i'*)
	}
// Non-cenrality parameter, ncp, is sum(`j'^2)
// Mean = df+ncp; SD = sqrt(2*(df+2*ncp)), so not all equal
summarize x*