// Normal

set seed 6834808

twoway (function y=normalden(x,0,1), range(-3 4)) ///
	(function y=normalden(x,0.5,1), range(-3 4)) ///
	(function y=normalden(x,1,1), range(-3 4)), ///
	title("Normal, d = 0.5") ytitle("density") ///
	legend(order(1 "mu=0" 2 "mu=0.5" 3 "mu=1"))

postfile kwpower test d T dfm dfr using kwpowermc.dta, replace
quietly forvalues d=0(0.1)1 {
	noisily display "difference " `d'
	forvalues i=1/2500 {
		clear
		set obs 60 // group size = 20
		generate group = mod(_n, 3)
		generate y = rnormal(group*`d', 1)
		anova y group
		post kwpower (1) (`d') (e(F)) (e(df_m)) (e(df_r))
		kwallis y, by(group)
		post kwpower (2) (`d') (r(chi2_adj)) (r(df)) (.)
	}
}
postclose kwpower

use kwpowermc.dta, clear

generate sig = Ftail(dfm, dfr, T)<0.05 if test==1
replace sig = chi2tail(dfm, T)<0.05 if test==2
bysort test: tabulate d sig, row

collapse (mean) sig, by(test d)
separate sig, by(test)
label variable sig1 "Anova"
label variable sig2 "Kruskal-Wallis"
save kwpowermcnorm.dta, replace
line sig1 sig2 d, title("Normal distributions") ytitle("Power") ///
	xtitle("group differences")