function [mu, s2, effekt, effektint, s2p, f] = faktorForsok(Y, alpha, printF, plotF)
%function [mu, s2, effekt, effektint] = faktorForsok(Y, alpha, printF, plotF)
%
% Skattar effekter och konfidensintervall i ett 2^k faktorförsök med n replikat
% för k = 2 eller 3 och n >= 2. För k=2 illustreras Även effekterna.
%
% Inparametrar:
%   Y - en matris med observationer där kolonnerna är observationer för de
%       olika faktorerna i standardordning. (1),(a),(b),(ab) för ett
%       2^2-försök. För ett 2^3-försök är ordningen
%       (1),(a),(b),(ab),(c),(ac),(bc),(abc). 
%   alpha - 1-konfidensgrad. alpha=0.05 om den inte anges.
%
%   printF - 1 för att skriva ut resultat/0 ingen utskrift
%    plotF - 1 för att illustrera effekterna i ett k=2 försök, 
%            0 (eller k=3) ingen illustration.
%
% Utparametrar:
%    mu - Medelvärden för var och och en av faktorkombinationerna.
%    s2 - Motsvarande variansskattningar.
%    effekt - Skattade effekter.
%    effektint - Motsvarande konfidensintervall.
%    s2p - poolad variansskattning
%      f - antal frihetsgrader

% 2015-09-11 Johan Lindström

% Delvis baserad på kod från facdesign.m av:
% 96-01-31, Orginalversion av okänd författare.
% 2003-10-03, Joakim Lübeck. Helt omarbetad och till stor del
%             korrigerad version.

if nargin<2 || isempty(alpha), alpha=0.05; end
if nargin<3 || isempty(printF), printF=1; end
if nargin<4 || isempty(plotF), plotF=1; end

%check input sizes
if size(Y,2)~=4 && size(Y,2)~=8
  error('Y måste ha 4 (k=2) eller 8 (k=3) kolonner')
end
if size(Y,1)<2
  error('Y måste ha minst 2 rader (observationer)')
end
n_fact = size(Y,2);
k = log2(n_fact);
n = size(Y,1);

%beräkna väntevärde och varians för varje faktor-komb
mu = mean(Y);
s2 = var(Y);
%varians
s2p = mean(s2);
%frihetsgrader
f = 4*(n-1);

%beräkna faktor-kombinationerna som en regression
X = faktorDesign(Y);
[effekt, effektint] = regress(Y(:), X, alpha);

kombstr = cell(1,n_fact);
kombstr(1:4) = {'(1)','a','b','ab'};
if k==3
  kombstr(5:end) = {'c','ac','bc','abc'};
end

if printF
  fprintf('\n 2^%u-Faktorförsök med %u observationer per faktor\n',k,n);
  fprintf('\n Kombination    Medel        s^2\n');
  for i = 1:n_fact
    fprintf(1, '%8s %12.3f %10.3f \n', kombstr{i}, mu(i), s2(i));
  end
   
  fprintf(1, '\n    Faktor     Effekt        %.3g%% konf.int\n', ...
          100*(1-alpha));
  for i = 2:n_fact
    fprintf(1, '%8s %12.3f %10.3f %10.3f\n', ...
     upper(kombstr{i}), effekt(i), effektint(i,:));
  end
  
  fprintf(1,'\n Varians och frihetsgrader\n');
  fprintf(1,'           n: %10u\n', n);
  fprintf(1,'           f: %10u\n', f);
  fprintf(1,'         s2p: %10.3f\n', s2p);
  fprintf(1,'   d(effekt): %10.3f\n', sqrt(s2p/(n*2^k)));
end

if plotF && k==2
  mu_all = effekt(1);
  A = effekt(2);
  B = effekt(3);
  AB= effekt(4);

  subplot(231)
  [a,b] = meshgrid([-1 0 1]);
  surf([-1 0 1],[-1 0 1],mu_all + A*a + B*b + AB*a.*b, 'facealpha', 0.5, ...
       'facecolor', 'interp', 'edgecolor', 'interp')  
  hold on
  for k=1:2
    for l=1:2
      plot3(2*k-3, 2*l-3, Y(:,k+2*(l-1)),'*k')
    end
  end
  hold off
  axis([-1.1 1.1 -1.1 1.1 min(Y(:)) max(Y(:))])
  view(-35,25)
  set(gca,'xtick',-1:1)
  set(gca,'xticklabel',{'-','','+'})
  set(gca,'ytick',-1:1)
  set(gca,'yticklabel',{'-','','+'})
  xlabel('Faktor A')
  ylabel('Faktor B')
  zlabel('Respons')

  subplot(234)
  plot([mu(1:2)' mu(3:4)'],'.-')
  text(2.1, mu(2), 'B-')
  text(2.1, mu(4), 'B+')
  axis([.75 2.25 ...
        min(mu(:))-.1*range(mu(:)) ...
        max(mu(:))+.1*range(mu(:))])
  caxis([min(mu) max(mu)])
  xlabel('A')
  ylabel('Respons')
  set(gca,'xtick',[1 2])
  set(gca,'xticklabel',{'-','+'})
  
  for i=1:4
    subplot(2,3,i+1+floor(i/3))
    surf([-1 0 1],[-1 0 1],mu_all + A*a + B*b + AB*a.*b, 'facealpha', 0.5, ...
         'facecolor', 'interp', 'edgecolor', 'interp')%, 'cdata', ones(3))
    hold on
    switch i
      case 1, Atmp = 0; Btmp=0; ABtmp=0; title_text = 'Medel';
      case 2, Atmp = A; Btmp=0; ABtmp=0; title_text = 'Effekt A';
      case 3, Atmp = 0; Btmp=B; ABtmp=0; title_text = 'Effekt B';
      case 4, Atmp = A; Btmp=B; ABtmp=0; title_text = 'Effekt A+B (utan samspel)';
    end
    surf([-1 0 1],[-1 0 1],mu_all + Atmp*a + Btmp*b + ABtmp*a.*b, ...
         'facealpha', 0.5, 'facecolor', 'interp', ...
         'edgecolor', .5*ones(1,3), 'cdata', .6*ones(3,3,3))  
    for k=1:2
      for l=1:2
        plot3(2*k-3, 2*l-3, Y(:,k+2*(l-1)),'*k')
      end
    end
    hold off
    axis([-1.1 1.1 -1.1 1.1 min(Y(:)) max(Y(:))])
    caxis([min(mu) max(mu)])
    view(-38,28)
    set(gca,'xtick',-1:1)
    set(gca,'xticklabel',{'-','','+'})
    set(gca,'ytick',-1:1)
    set(gca,'yticklabel',{'-','','+'})
    xlabel('Faktor A')
    ylabel('Faktor B')
    zlabel('Respons')
    title(title_text)
  end
  
  colormap( jet(128) )
end
