function rjmcmc
clear all
clc
N=100000;                                    
m=10;                                       
k=10;                                       
sigma0=.2;
sigmap=.3;
sigmar=.2;

%initial conditions
n0=ceil(rand*m);                      
mub=ones(1,n0)*2;                     
b=normrnd(2,sigmap,n0,1);             
X=5*randn(k,m);
y=X(:,1:n0)*b+sigma0*randn(k,1);
tic 
for i=1:N
nstar=ceil(rand*m);                    
mustar=ones(1,nstar)*2;                 
        if (nstar<n0)                       
            u1=normrnd(0,sigma0,nstar,1);   
            b1=b(1:length(u1));             
            b2=b(length(u1)+1:end);         
            bnstar=b1+u1;                   
            u=[u1;b2];                      
            Unorm = sum(u(:).^2);
            h2 = ((1/sqrt(2*pi*sigmar^2))^n0)*exp((-1/(2*sigmar^2))*(Unorm));
            U1norm = sum(u1(:).^2);
            h1 = ((1/sqrt(2*pi*sigmar^2))^nstar)*exp((-1/(2*sigmar^2))*(U1norm));
   
        else                                
            u=normrnd(0,sigma0,nstar,1);    
            u1=u(1:n0);                     
            u2=u(n0+1:end);                
            bn=[b;zeros(length(u2),1)];     
            bnstar=bn+u;                    
            Unorm = sum(u(:).^2);
            h1 = ((1/sqrt(2*pi*sigmar^2))^nstar)*exp((-1/(2*sigmar^2))*(Unorm));
            U1norm = sum(u1(:).^2);
            h2 = ((1/sqrt(2*pi*sigmar^2))^n0)*exp((-1/(2*sigmar^2))*(U1norm));
        end
    Xstar=X(1:nstar,:);
    Xn=X(1:n0,:);
    y1=y'-bnstar'*Xstar;
    y2=y'-b'*Xn;
    y1 = sum(y1(:).^2);
    y2 = sum(y2(:).^2);
    E1=y1/(2*sigma0^2);
    E2=y2/(2*sigma0^2);
    E3 = bnstar-mustar';
    E3 = sum(E3(:).^2);
    E4 = b-mub';
    E4 = sum(E4(:).^2)';
    rho1=(exp(-(E1-E2)))*((2*pi*sigmap^2)^((n0-nstar)/2))*(exp((-1/(2*sigmap^2))*(E3-E4)))*(h2/h1);
    rho=min(rho1,1);
 
    U=rand;
        if(U<rho)
            n(i+1)=nstar;
        else
            n(i+1)=n0;
        end
end
toc
n;
figure(1)
hist(n)                         
h=hist(n);                      
idx=find(h==max(h))           
b

              

        
    