function [x,y]=bayes(N,rho,alpha,sigma,y)
% [x,y]=bayes(N,rho,alpha,sigma,y)
% Inputs:
% N - number of samples each iteration
% rho - fraction of best performing samples to take
% alpha - smoothing paramter
% sigma - std. deviation for image noise
% y - image data
% Outputs:
% x - reconstructed image
% y - original image

% mu's known, sigma known!

     mm=0;
mu=[0,1]';
beta=[.5,.5,.5,.5]';
if nargin<=4,
% example image
     y=zeros(20);
     y(6,8:12)=ones(1,5);
     y(7,7:13)=ones(1,7);
     y(8:12,6:14)=ones(5,9);
     y(13,7:13)=ones(1,7);
     y(14,8:12)=ones(1,5);
     if nargin<=3,
     sigma=.01;
     end
     if nargin==0,
     N=400;
     rho=.05;
     alpha=.5;
     end
end

     tol=0.01;

count=1; % start counter
yold=y;

% add noise
     y=y+mm+sigma*randn(size(y));

y(y>1)=1; % truncate image to [0,1]
y(y<0)=0;
%y=round(y);

% initial probs. --> equal prob of being 0 or 1
[m,n]=size(y);
v=ones(m,n)*1/(2);
vold=zeros(m,n);

while(max(max(abs(v-vold)))>tol) % have our probs changed less than tol since last it?
     [count,scoreb((v>=.5),y,mu,sigma,beta)] % display information about current iteration
V=[];
S=[];
for i=1:N,
% generate x's
	V(:,:,i)=(rand(m,n)<v);
% score x's
S(i)=scoreb(V(:,:,i),y,mu,sigma,beta);
end
% update parameters
vold=v;
[S,I]=sort(S);
g=floor((rho)*N);
w=(1/g)*sum(V(:,:,I(1:g)),3);
v=alpha*w+(1-alpha)*v;
count=count+1;
end
x=(v>=.5);
scoreb(x,y,mu,sigma,beta)
scoreb(yold,y,mu,sigma,beta)
