function [phi, accept] = codem(n, x, phi, sigma2, accept)

% step 1
theta = log(1 + phi) - log(1 -phi);
thetacan = normrnd(theta, sqrt(sigma2));

% step 2
phican = (exp(thetacan)-1)/(exp(thetacan)+1);

% step 3
sumx = sum(x);
logcan = sumx*log(phican) + (n -sumx)*log(1-phican) + thetacan -2*log(1 +exp(thetacan));
logold =  sumx*log(phi) + (n -sumx)*log(1-phi)+ theta-2*log(1 +exp(theta));
logf = logcan - logold;

% step 4
u = rand;

% step 5
if log(u)<logf
    phi = phican;
    accept  = accept  + 1;
end
end