% A small-scale simulation to explore
% the sensitivity and specificity of
% the weighted $t$-test proposed in the
% Bioinformatics article.


%%%%%%%%%%%%%%%%%%%%%%%%
% Set a few constants, allocate space for results
%%%%%%%%%%%%%%%%%%%%%%%%
  
  n = 50000; % the common size of the libraries
  n_sim = 10000; % the number of simulated comparisons
  n_A = 4;   % number of libraries in group 1
  n_B = 4;   % number of libraries in group 2

  p_A = [1 2 5 10 100]/n;  % the mean proportion in group 1
  p_B_mult = [1 2 3 5 10]; % the scaling factor for group 2 (fold change)
  
  k_over = [2 5 10 50]; % factor by which the beta-binomial variance
			% should exceed the binomial variance. 

  resmat = zeros(5,5,4,12);
  
%%%%%%%%%%%%%%%%%%%%%%%%
% Set a few constants
%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%
% Main processing loop
%%%%%%%%%%%%%%%%%%%%%%%%

for(i1 = 1:length(p_A))
  zed = clock;
  [i1 zed(4:6)]
  for(i2 = 1:length(p_B_mult))
    for(i3 = 1:length(k_over))
      A = bioinf_sim_fn(n,p_A(i1),p_B_mult(i2),...
			n_A,n_B,k_over(i3),n_sim);
      resmat(i1,i2,i3,[1 5 9]) = sum(A(:,[2 5 8]) < 0.1); 
      resmat(i1,i2,i3,[2 6 10]) = sum(A(:,[2 5 8]) < 0.05); 
      resmat(i1,i2,i3,[3 7 11]) = sum(A(:,[2 5 8]) < 0.01); 
      resmat(i1,i2,i3,[4 8 12]) = sum(A(:,[2 5 8]) < 0.005);
    end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%
% Main processing loop
%%%%%%%%%%%%%%%%%%%%%%%%