function[phat,vhat,v_min,alphahat,betahat] = ...
    fit_ab2(n_vec,x_vec)
% FIT_AB2(N_VEC,X_VEC) -- Given the vectors of library sizes
% (n_vec) and associated counts for a particular tag (x_vec),
% this fits the alpha and beta values of the underlying beta
% distribution using the method of moments. It returns 
% PHAT - the estimated proportion
% VHAT - the associated variance estimate
% VMIN - the minimum possible variance
% ALPHAHAT - the first of the beta distribution parameters
% BETAHAT - the second of the beta distribution parameters
%
% n_vec and x_vec are assumed to be 1 by n.

% Get the associated proportions

p_vec = x_vec ./ n_vec;

% Get the minimum variance

v_min = sum(p_vec .* (1- p_vec) ./ n_vec);

% assume a starting weight vector of all equal values

n = length(n_vec);
w_vec = ones(size(n_vec))/n;
%w_vec = sum(n_vec)./n_vec;
w_vec = n_vec ./ sum(n_vec);
w_vec = w_vec / sum(w_vec);

phat = w_vec * p_vec';
v_min = phat*(1-phat)/sum(n_vec);

% Iterate to convergence (this is quite rapid)

n_rep = 15;

resmat1 = zeros(n_rep,4);
resmat2 = zeros(n_rep,n);
k_vec = zeros(2,1);
k_vec2 = zeros(4,1);
counter = 1;

for(i = 1:n_rep)

  phat = w_vec * p_vec';
  w2 = w_vec*w_vec';
  vhat = (((w_vec.^2) * (p_vec.^2)') - w2*(phat^2)) / (1 - w2);
  
  if(vhat > v_min)
    
    % Fit beta and alpha
    
    w2n = sum(w_vec.^2 ./ n_vec);
    betahat = (vhat - phat*(1-phat)*w2) / (phat*w2n - (vhat/(1-phat)));
    alphahat = (phat/(1-phat))*betahat;
    
    % refit the weight vector
  
    w_vec = ((alphahat + betahat)*n_vec) ./ ...
	    ((alphahat + betahat + n_vec));
    w_vec = w_vec ./ sum(w_vec);
    
  else
    
    alphahat = 0;
    betahat = 0;
    vhat = v_min;
    break;
    
  end
  
%  resmat1(i,:) = [phat vhat betahat alphahat];
%  resmat2(i,:) = w_vec;

  k_vec(mod(i+1,2)+1) = alphahat + betahat;
  if(~(mod(i,2)))
    k_temp = mean(k_vec);
    w_vec = (k_temp*n_vec) ./ (k_temp + n_vec);
    w_vec = w_vec ./ sum(w_vec);
  end

%  k_vec2(counter) = alphahat + betahat;
%  if(~(mod(counter,4)))
%    delta = (k_vec2(1)-k_vec2(2)) / ...
%	    ((k_vec2(1)-k_vec2(2)) - (k_vec2(3)-k_vec2(4)));
%    k_temp = k_vec2(1) + delta*(k_vec2(3)-k_vec2(1))
%    w_vec = (k_temp*n_vec) ./ (k_temp + n_vec);
%    w_vec = w_vec ./ sum(w_vec);
%    counter = 0;
%  end
%  counter = counter + 1;

  
end

%resmat1(:,1)'
%resmat1(:,2)'
%resmat1(:,3)'
%resmat1(:,4)'
%resmat2