%
%       FBST - MONTE CARLO
%    SEPARE ESTE ARQUIVO EM 4, CONFORME INSTRUCOES ABAIXO:


%%%%%%%%%%%%%%%%%%%%%
%
%   Arquivo 1/4: driverhomo.m
%

% TESTE DE HOMOGENEIDADE
% Problema:
% dada uma tabela de contingencia
%   X = [x11, x12;
%        x21, x22]
% onde x11 = numero de sucessos na amostra do 1o grupo
%      x12 = numero de fracassos na amostra do 1o grupo
%      x21 = numero de sucessos na amostra do 2o grupo
%      x22 = numero de fracassos na amostra do 2o grupo
%
% Testar a hipotese
%       H: p1 = p2
% onde p1 e p2 são as proporções de sucessos no 1o e 2o grupo, respectivamente.
%
X = [150 250; 
   70 130];
omega = X+1;
%omega = [1,1; 1,1];
sx = sum(X,1);
p0 = sx(1) / sum(sx);
theta0 = [p0, p0];

f0 = fxhomo(theta0,X);

resp= [1,0,0,0,0,0,0,0,0,3,50000,0.005,0,0,0];

resp= mtcar('fxhomo', 'gsamphomo', omega, f0, resp, X);

% v: evidencia a favor da hipotese H.
v = 1-resp(13),


%%%%%%%%%%%%%%%%%%%%%
%
%   Arquivo 2/4: fxhomo.m
%
function f = fxhomo(theta,X)

f = exp( log(theta(1))*X(1,1) + log(1-theta(1))*X(1,2) + ...
         log(theta(2))*X(2,1) + log(1-theta(2))*X(2,2) );
      
      
      
%%%%%%%%%%%%%%%%%%%%%
%
%   Arquivo 3/4: gsamphomo.m
%
function [theta, g, omega] = gsamphomo(omega,X)

% omega = hiperparametros da beta
theta = [betarnd(omega(1,1),omega(1,2)), betarnd(omega(2,1),omega(2,2))] ;
   
g = exp( log(theta(1))*omega(1,1) + log(1-theta(1))*omega(1,2) + ...
         log(theta(2))*omega(2,1) + log(1-theta(2))*omega(2,2));
      
      
      
%%%%%%%%%%%%%%%%%%%%%
%
%   Arquivo 4/4: mtcar.m
%
      
 function resp= mtcar(fdens, gsamp, omega, f0, resp, X)
 % fdens(theta|X) - original pdf/likelihood
 % prototype:
 %    f = fdens(X,theta)
 %			X: data; theta: parameter vector
 %
 % gsamp(theta|X, omega) - importance sampling generating function
 % prototype:
 %    [theta, g, omega] = gsamp(X,omega,j,nst)
 %			X: data; omega: hyperparameters; nst: number of strata; j:current stratum; 
 %       theta: parameter vector; g: G(theta)
 
 % resp= [xid,k,ah,aw,awh,awc,aw2,aw2h,aw2c,nst,nmc,targ,v,delta,chi2] 
 % arg      1 2  3  4   5   6   7    8    9  10  11   12 13   14   15 
 % Accumul, Mean, Std, Hit (f>f0), Complement, eVidence 
 % v= eVidence --against-- the hypothesis  
 % [mx,nx]= size(Xs);  
 % mx= Poissom pairs per X , 
 % nx= X samples in Xs for quantile computation   
 % X= Xs(1:mx,[2*xid-1,2*xid]);  
 % nmc= monte carlo points (t1,t2) 
 % nst= stratified points (t3) per (t1,t2) 
 % mx= Poissom pairs per X , 
 % nx= X samples in Xs for quantile computation   

 % Fieller and Finney confidence interval:  
 % sv= Does not exist! N(0,1)/N(0,1) is a Cauchy! 
 % v= mwh/mw = mwh/(mwh+mwc) , wh*wc=0 ,   
 % Cov(A*X) = A*Cov(X)*A' , Cov(X,Y) = E(X*Y)-E(X)*E(Y)    
 % Cov(wh,wc) = E(wh*wc) -E(wh)*E(wc) = -mwh*mwc 
 % Cov(wh,w) = Cov(wh,(wh+wc)) = swh^2  +Cov(wh,wc)    
 % sw^2 = Var(wh+wc) = swh^2 +swc^2   -2*mwh*mwc 
 % Cov(wh,w) = sw^2 -swc^2 +mwh*mwc 
 % v_n= mwh/mw  , u_n= wh -v*w , mu_n= mwh -v*mw 
 % E(u) = 0 , (1/n)*Var(wh -v*w) =  Var(mu) = 
 % =(1/n)*(sw1^2 -2*v*(sw^2-swc^2+mwh*mwc) +v^2*sw^2)  
 % n*(mu)^2/Var(mu)  ~ Chi2(r=1)  as n -> \infty  
 % => A*v^2 -2*B*v +C = 0 
 % A= mw^2 -sw^2*chi2/n  
 % B= mw*mwh -chi2*(sw^2 -swc^2 +mwh*mwc)/n  
 % C= mwh^2 -swh^2*chi2/n  
 % v= (2*B +-sqrt(4*B^2 -4*A*C))/(2*A) 
 % So, a confidence interval: v +-delta  
 % delta= sqrt(B^2 -A*C)/A  ~= sqrt(chi2*D/n) 
 % where D= rwh^2*(1-v)^2 +rwc^2*v^2 
 % on the coeficients (of vaRiation)     
 % rw= sw/mw , rwh= swh/mw , rwc= swc/mw   
 % Solving for n, 
 % n= chi2*D/d^2 , chi2= Chi2(df=1, 1- alpha)    
 % define chi2 (1-alpha)-quantile of Chi2(d.f.=1)
 % alpha     0.005  0.01   0.025  0.05  0.1  
 % (1-a)-q   7.88   6.63   5.02   3.84  2.71 
 
 EPS = 1e-5;
 
 %theta = zeros(1,dim);
 
 xid=  resp(1);   mx=   size(X,1);  
 %alternat input: X= Xs(1:mx,(2*xid-1):(2*xid));   
 k=    resp(2);   ah=   resp(3);   %ac= k -ah, 
 aw=   resp(4);   awh=  resp(5);   awc=  resp(6); 
 aw2=  resp(7);   aw2h= resp(8);   aw2c= resp(9); 
 nst=  resp(10);  nmc=  resp(11);  targ= resp(12); 
 v= resp(13);     delta= resp(14); chi2= resp(15); 
 if( chi2==0 )  chi2= 5.02;   end%chi2  6.7;
 if( k==0 )     mink= 1500;  else   mink=0;  end%mink  


 %t0= mom(X);   t0, 
 %MLE and moments the same for H0 
 %no need for optimization step 
 %usualy f0 is input argument 
 %trivial case: 
 %if( t0(3)==0 )                % trivial case, X over hyp   
 %  resp= [xid,1,zeros(1,13)];  % evidence against hyp = 0 
 %  return;                     % goto is beautiful  
 %end                   
 %t= [t0(1),t0(2),0]; 
 %f0= holrec(t,X,Glx); 
 
 
 % Importance sampling: g a prob density:  
 % Sh(y|x)dF = Sh(y|x)f(x)dx = 
 % Sh(y|x)(f(x)/g(x))g(x)dx = Sh(y|x)(f/g)dG 
 % MC: draw X_j ~g(x)  w_j= f(x_j)/g(x_j) 
 % S ~= sum{j=1:n}(w_j*h(y|x))/sum{j=1:n}(w_j) 
 
 
 %posterior/likelihood = import. sampl. parameters  
 
 pack= 500; 
 flag= 1; 
 wjant = 0;
 while( flag==1 ) 
   k= k+1;   
   
   for j=1:nst    %nst strat. samp. for theta
      
      linsamp = ['[theta, gj, omega] = ', gsamp, '(omega,X);'];
      eval(linsamp);
      lindens = ['fj = ', fdens, '(theta,X);'];
      eval(lindens);
      
      wj=  fj/gj;
      
      
      %if (wj~=wjant)
      %   [wj, wjant, fj, gj],
      %   pause(0.1)
      %   wjant=wj;
      %end
      
      %[fj, gj, wj],
      %pause(0.3);
      
      wj2= wj*wj;     
      hj= ( fj >= f0 );  cj= (1-hj);    
      ah=   ah  +hj; 
      aw=   aw  +wj;  awh=  awh +hj*wj;   awc=  awc +cj*wj; 
      aw2=  aw2 +wj2; aw2h= aw2h +hj*wj2; aw2c= aw2c +cj*wj2; 
      
   end%for_j 
   if( mod(k,pack)==0 )     
     n=    k*nst;  sqrtn= sqrt(n); 
     mh=   ah/n; 
     mw=   aw/n;   mwh=  awh/n;   mwc=  awc/n; 
     mw2=  aw2/n;  mw2h= aw2h/n;  mw2c= aw2c/n; 
     
     sh=   sqrt(max(mh  -mh^2,0)); 
     sw=   sqrt(max(mw2 -mw^2,0)); 
     swh=  sqrt(max(mw2h -mwh^2,0));  
     swc=  sqrt(max(mw2c -mwc^2,0));  
     vh=   mh; 
     v=    mwh/mw; 
     rw=   sw/mw;  rwh= swh/mw;  rwc= swc/mw;    
     A= mw^2 -sw^2*chi2/n;   
     B= mw*mwh -chi2*(sw^2 -swc^2 +mwh*mwc)/n;  
     C= mwh^2 -swh^2*chi2/n;   
     delta= sqrt(max(B^2 -A*C,0))/A; 
     D= rwh^2*(1-v)^2 +rwc^2*v^2;  
     delta2= sqrt(chi2*D/n); 
     [xid,k],
     [v,delta,delta2], %[rw,rwh,rwc], %monitor*** 
     %if (delta<0)
     %   [wj, wj2, ah, aw, awh, awc, aw2, aw2c, sh, sw, swh, swc, v, A, B, C, D],
     %   [A, B, C],
     %   pause(2);
     %end
          
     %[mw2, mw, mw2], 
     if( (k >= mink) & (k >= nmc | (delta>=0 & delta<targ) ) ) 
       flag= 0; 
     end%if_flag  
   end%if_mod_100 
 end%while 
 resp= [xid,k,ah,aw,awh,awc,aw2,aw2h,aw2c,nst,nmc,targ,v,delta,chi2];  
 %arg     1 2  3  4   5   6   7    8    9  10  11   12 13   14   15  
            
 %******************************************************************* 
