% code by Ata Kaban 2003, edited by Ella Bingham 2005 (based on bin_plsa_new.m)

% The Aspect Bernoulli method is described in:
% Ata Kaban, Ella Bingham and Teemu Hirsimaki: 
% "Learning to read between the lines: The aspect Bernoulli model", 
% Proceedings of the 4th SIAM International Conference on Data Mining, 
% April 22-24, 2004, Lake Buena Vista, Florida, pp. 462-466.
% and
% Ella Bingham, Ata Kaban and Mikael Fortelius:
% "The Aspect Bernoulli model: Multiple causes of presences and absences". 
% Pattern Analysis and Applications. To appear.

function[A,S,l,allpk]=aspect_bernoulli(x,K,iter,posterior);

% Inputs: 
% x: binary data matrix; each column is an observation
% K: number of aspects (components) to estimate
% iter: maximum number of iterations to make
% posterior: 1/0 to compute the posterior of components or not. Slow; use only if needed. Default=0.

% Outputs:
% A: matrix of parameters, A(t,k) = prob(1|attribute t, aspect k)
% S: matrix of parameters, S(k,n) = prob(aspect k| observation n)
% l: vector of log likelihood values in the iterations
% allpk: cell array of posteriors p(k|t,n,x_{tn}) of all components k 

if nargin < 4; posterior=0; allpk = []; end
[T,N]=size(x);
S=rand(K,N);S=S./repmat(sum(S),K,1);
A=rand(T,K);
A1=1-A;
l=[];

for i=1:iter
   S=S.*(A'*(x./max(eps,A*S))) + S.*((1-A)'*((1-x)./max(eps,1-A*S))); 
   S=S./repmat(sum(S),K,1);  
   A=A.*(( x ./ max(eps,A*S) )*S');
   A1=A1.*(( (1-x) ./ max(eps,A1*S) )*S');
   A=A./(A+A1);
   A1=1-A;
   l=[l, sum(sum(x.*log(max(eps,A*S))+(1-x).*log(max(eps,1-A*S))))/(N*T)];
   if i>1 & abs(l(i)-l(i-1))<1e-5; break; end
end   

% Posterior p(k|n,t,x_{nt}) of component k at observation n and attribute t
if posterior
[T,N]=size(x);
allpk=cell(1,K);
pk = zeros(K,1);
for n=1:N, for t=1:T
    for k=1:K
      pk(k) = S(k,n)*A(t,k)^x(t,n)*(1-A(t,k))^(1-x(t,n));
    end;
    pk = pk/sum(pk);
    for k=1:K
      allpk{k}(t,n)=pk(k);
    end;
end; end;
else allpk = [];
end