%Welch ANOVA Simulations
%Example from the original paper Welch 1951 page 335
clear all

k=3;
%
yibar=[27.8 24.1 22.2];
si2=[60.1 6.3 15.4]
%ni =[20 10 10] %Original in Welch (1951).
ni=[28  10   7]
wei =ni./si2 
nwei=wei/sum(wei) 
yhat= sum(nwei.*yibar)
Q=sum( (1-nwei).^2./(ni-1))
W=sum(wei.*(yibar - yhat).^2)/(k-1) 
v2 = W/(1+2*(k-2)*Q/(k^2 -1)) 
nu = (k^2-1)/(3 * Q) 
pval = 1-fcdf(v2, k-1, nu) 
%%
% Standard ANOVA (See hint to Exercise 11.4 (MTHFR)
% pn how the fond F from \bar Y_i  and s_i^2's)
ybar=sum(ni.*yibar)/sum(ni)
MSTr= 1/(k-1) * sum( ni .* (yibar - ybar).^2 )
MSE = 1/(sum(ni)-k) * sum( (ni-1).*si2)
pval=1-fcdf(MSTr/MSE, k-1, sum(ni)-k)
 
