% Welch's  ANOVA
clear all
close all force
  disp('Welch''s ANOVA Simulations')
  lw =2.5; 
  set(0, 'DefaultAxesFontSize', 16);
  fs = 16;
  msize = 10;

  s = RandStream('mt19937ar','Seed',0);
 RandStream.setGlobalStream(s);
 %


%
format long
k=4; nk=200;
pvals=[]; polds=[];
M=50000;
for i = 1:M
yy=randn(k, nk);
    y1=3*yy(1,:); y2= yy(2,:); ...
    y3 = yy(3,:); y4= yy(4,:);
yibar=[mean(y1) mean(y2) mean(y3)  mean(y4)]; 
  ni=[nk nk nk nk];
  si2 = [var(y1) var(y2) var(y3) var(y4)];
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) ;
pvals=[pvals pval];
%
pold = anova1([y1 y2 y3 y4], ...
    [ones(1,nk)  2*ones(1,nk)  3*ones(1,nk)  4*ones(1,nk)],'off');
polds=[polds pold];
end
figure;
bihist( pvals,  polds, 50) 
legend('Welch','ANOVA')
print -depsc 'C:\STAT\ANOVA\ANOVAeps\welchpvals.eps'
anovap=sum(polds<0.05)/M
welchp=sum(pvals<0.05)/M
  %anovap=  0.08054
  %welchp = 0.05014
format short