function [x,xx] = blogis (y,iters,x0)
% x = blogis (y,iters,x0)
% Bayesian Logistic regression with a Gibbs sampler.
% $Id: blogis.m,v 1.3 2017/04/29 13:58:43 daichi Exp $
if nargin < 3
  x0 = 0;
  if nargin < 2
    iters = 100;
  end
end
% body
x = x0;
N = length(y);
u = zeros(1,N);
for j = 1:iters
  % sample u
  s = sigmoid (x);
  for i = 1:N
    if y(i) > 0
      r = unif (0,s);
    else
      r = unif (s,1);
    end
    u(i) = log (r / (1 - r));
  end
  % determine border
  if (sum(y==1) > 0)
    a = max (u(y == 1));
  else
    a = -inf;
  end
  if (sum(y==0) > 0)
    b = min (u(y == 0));
  else
    b = inf;
  end
  % sample x
  if a < b
    while (1)
      newx = randn;
      if (a < newx & newx < b)
	x = newx;
	break;
      end
    end
  else
    fprintf(1,'iter %d: max = %g, min = %g\n', j, a, b);
  end
  xx(j) = x;
end
x = mean(xx);
hist(xx,20);
axis([-4 4 0 floor(iters/5)]);
