function [r,bestC] = DPgibbs(X,alpha,G0alpha,maxIter,initFull,Y)
  % X        data set
  % X(n,w) = # of times word w appears in document n
  % alpha    concentration parameter for DP
  % G0alpha  concentration parameter for G0 prior
  % maxIter  maximum gibbs iterations
    
[N,W] = size(X);

% initialize
c = ones(N,1);
K = 1;

if initFull,
  c = 1:N;
  K = N;
end;

bestC = c;
bestV = -Inf;

r = zeros(2,maxIter);
maxK = floor(10*max(10,log(N)));
%logStir = computeLogStirlingTable(N,maxK);
%for n=1:N,
%  logStir(n,:) = logStir(n,:) - addLog(logStir(n,:));
%end;
tic;

subTime = 0;
% iterate
for iter=1:maxIter,
  order=permuted([1:N]');
  for nn=1:N,
    n = order(nn);
    
    lp_new = log(alpha) + log_marginal_posterior_normal(X(n,:), G0alpha, []);
    lp_old = ones(K,1) * -Inf;
    for k=1:K,
      d = c; d(n) = -1;
      ck = find(d==k);
      if length(ck) > 0,
        lp_old(k) = log(length(ck)) + ...
                    log_marginal_posterior_normal(X(n,:), G0alpha, X(ck,:));
      end;
    end;

    sm = addLog([lp_new ; lp_old]);
    p  = random_position(exp([lp_new ; lp_old] - sm)) - 1;
    if p == 0,
      K = K+1;
      c(n) = K;
    else
      c(n) = p;
    end;
  end;
  
  tm = toc;
  lp = data_likelihood_normal(X,G0alpha,c) + ...
       log_DP_prior_count(c,alpha,N);
%       logStir(N,min(maxK,length(unique(c))));
  %clf;
  %plot_colors(X,c);
  %[fs pr re] = compute_f(Y, c);
  fs = 0; pr = 0; re = 0;
  if mod(iter,1) == 0,
    disp(sprintf('%d  \t%g  \t\t%g  \t%g', length(unique(c)), lp, tm-subTime,fs));
  end;
    r(1,iter) = tm-subTime; r(2,iter) = lp; r(3,iter) = fs; r(4,iter) = pr; ...
        r(5,iter) = fs;
  if lp > bestV,
    bestV = lp;
    bestC = c;
  end;

  endTime = toc;
  subTime = subTime + (endTime-tm);
end;


function p = random_position(a)
  r = rand();
  p = 0;
  for i=1:length(a);
    if r < a(i),
      p = i;
      break;
    else
      r = r - a(i);
    end;
  end;
  
  
