function [x,y,v]=bayes2(N,rho,alpha,sigma,y)
% [x,y,v]=bayes2(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
% v - probabilities at the end

% mu's unknown, sigma known! (maybe --> sigma unknown)

     mm=0; % mean of noise to add to image
mu=[.1,.9]'; % initial mu
beta=[.5,.5,.5,.5]';
if nargin<=4,
	 % create 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=1; % initial sigma
     end
end
     tol=0.01;
count=1; % start iteration counter
yold=y; 

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

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

while(max(max(abs(v-vold)))>tol) % while our probs have change less than tol
     [count,scoreb((v<=.5),y,mu,sigma,beta)] % crude display
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;

%update mu estimate

mu(1)=alpha*sum(sum(sum(V(:,:,1:g),3).*y./g))/(m*n)+(1-alpha)*mu(1);
mu(2)=alpha*sum(sum(sum((1-V(:,:,1:g)),3).*y./g))/(m*n)+(1-alpha)*mu(2);

count=count+1;
end

x=(v>.5);
scoreb(x,y,mu,sigma,beta)
scoreb(yold,y,mu,sigma,beta)
mu
