function gibbs_nonlinear
close all
%% Plot target
xv = -1.5:0.01:1.5;
yv = -2:0.01:2;
[xx,yy] = meshgrid(xv,yv);
a = 1;
b = 0.1;
z = logpi([xx(:) yy(:)]', a,b);
z = reshape(z, length(yv), length(xv));
contour(xv,yv,exp(z));
axis xy

%% Run sampler
M = 150;
X1 = zeros(1,M);
X2 = zeros(1,M);
sq = 0.5;
hold on;
plot(X1(1),X2(1),'r.','MarkerSize',12);
xlabel('x_1');
ylabel('x_2');
set(gca,'FontSize',12);
pause(0.1);
ar = 0;
for(m = 2:M)
    X1p = X1(m-1) + sq*randn(1);
    ap = logpi([X1p;X2(m-1)],a,b) - logpi([X1(m-1);X2(m-1)],a,b);
    if(rand(1) < exp(ap))
        X1(m) = X1p;
        ar = ar+1;
        
        plot([X1(m-1) X1(m)],[X2(m-1) X2(m-1)],'r.','MarkerSize',12);
    else
        X1(m) = X1(m-1);
        plot(X1p,X2(m-1),'kx','linewidth',1,'Markersize',8);      
    end
    pause(0.1)
    
    X2(m) = X1(m)^3 + sqrt(b)*randn();
    plot([X1(m) X1(m)],[X2(m-1) X2(m)],'r.','MarkerSize',9);
    pause(0.1)
end

end
%----------------------------------
function z = logpi(x,a,b)
z = 1/a*(2*x(1,:)+sin(x(1,:)*2*pi)).^2 + 1./(b) .* (x(2,:) - x(1,:).^3).^2;
z = -1/2*z;
end