function [prec1,cprec1,beta1,zprec1,p1,kk1,startp] = clustering(y,z,prec,cprec,beta,zprec,p,kk,alpha,K,start)

persistent ptune ctune ztune N J t1 bcovtune vdat;

if start==0;
    ptune=50;
    ctune=20;
    ztune=2;
    sizey=size(y);
    N=sizey(1); J=sizey(2);
    t1=sum(sum(y.*y,1));
    bcovtune=diag(ones(1,K))*0.2;
    vdat=y'*y;
end;
startp=1;

nk=sum(kk,1);
dkk=y*kk;  %NxK

%find prec1 using M-H;
pcand=gamrnd(ptune,prec/ptune);
a=gampdf(pcand,ptune,prec/ptune)/gampdf(prec,ptune,pcand/ptune);
rc=cprec/pcand;
r0=cprec/prec;
tc=sum((dkk.^2)*(1./(rc+nk)'));
t0=sum((dkk.^2)*(1./(r0+nk)'));
hc=sum(((beta.^2)./(rc+nk))/pcand)+1/zprec;
h0=sum(((beta.^2)./(r0+nk))/prec)+1/zprec;
uc=(0.5/hc)*sum(((dkk*(beta./(rc+nk))')-z).^2); %final sum is across patients
u0=(0.5/h0)*sum(((dkk*(beta./(r0+nk))')-z).^2);
lr=exp(prec*(t1-t0)/2 - pcand*(t1-tc)/2 - uc + u0) / (((pcand/prec)^(-J))*((rc/r0)^(-K))*(hc/h0)*prod(rc+nk)/prod(r0+nk))^(N/2);
prior=(1/pcand)^(3/2)/(1/prec)^(3/2);
if rand<(lr*prior)/a
    prec1=pcand;
else
    prec1=prec;
end;

%find cprec1 using M-H;
ccand=gamrnd(ctune,cprec/ctune);
a=gampdf(ccand,ctune,cprec/ctune)/gampdf(cprec,ctune,ccand/ctune);
rc=ccand/prec1;
r0=cprec/prec1;
tc=sum((dkk.^2)*(1./(rc+nk)'));
t0=sum((dkk.^2)*(1./(r0+nk)'));
hc=sum((beta.^2)./(prec1*(rc+nk)))+1/zprec;
h0=sum((beta.^2)./(prec1*(r0+nk)))+1/zprec;
uc=(0.5/hc)*sum(((dkk*(beta./(rc+nk))')-z).^2); %final sum is across patients
u0=(0.5/h0)*sum(((dkk*(beta./(r0+nk))')-z).^2);
lr=exp(prec1*(tc-t0)/2 - uc + u0) / (((rc/r0)^(-K))*(hc/h0)*prod(rc+nk)/prod(r0+nk))^(N/2);
prior=(1/ccand)^(3/2)/(1/cprec)^(3/2);
if rand<(lr*prior)/a
    cprec1=ccand;
else
    cprec1=cprec;
end;

r=cprec1/prec1;
betap=beta;

for j=randperm(J);         %for gene j
    cur=kk(j,:)*(1:K)';    %original cluster of gene j (ie. cluster 3);
    ddelt=zeros(1,K);
    
    vco=kk(:,cur);         %vector of original cluster before move (ie. 1 1 0 0 0 0)'
    vcn=vco;
    vcn(j)=0;              %vector of original cluster after move
    betac=beta(cur);       %beta for original cluster
    mco=nk(cur);           %size of original cluster before move
    mcn=mco-1;             %size of original cluster after move
    ho=sum((beta.^2)./(prec1*(r+nk)))+1/zprec;      %h before move
    dkk=y*kk; 

    for k=1:K;             %move gene j to cluster k
        if k~=cur
            vno=kk(:,k);   %vector of cluster k before move
            vnn=vno;
            vnn(j)=1;      %vector of cluster k after move
            betan=beta(k); %beta for cluster k
            mno=nk(k);     %size of cluster k before move
            mnn=mno+1;     %size of cluster k after move
            hc=ho-(betac^2)/(prec1*(r+mco))-(betan^2)/(prec1*(r+mno))+(betac^2)/(prec1*(r+mcn))+(betan^2)/(prec1*(r+mnn));
            ku=kk; ku(j,cur)=0; ku(j,k)=1;
            dkkmove=y*ku;
            nkmove=nk; nkmove(cur)=nk(cur)-1; nkmove(k)=nk(k)+1;

            ddelt(k)=p(k)*exp(0.5*prec1*(vcn'*vdat*vcn/(r+mcn) + vnn'*vdat*vnn/(r+mnn)) ...
                - (0.5/hc)*sum(((dkkmove*(beta./(r+nkmove))')-z).^2) ...
                - 0.5*prec1*(vco'*vdat*vco/(r+mco) + vno'*vdat*vno/(r+mno)) ...
                + (0.5/ho)*sum(((dkk*(beta./(r+nk))')-z).^2) ...            
                - (N/2)*(log(hc) - log(ho) + log(r+mcn) + log(r+mnn) - log(r+mco) - log(r+mno)));
            
        else
            ddelt(k)=p(k);
        end       
    end
    
    newcl=sum((cumsum([0 ddelt/sum(ddelt)])<rand));
    
    betaprop=betap;
    if newcl~=cur
        kk(j,newcl)=1;
        kk(j,cur)=0;
        nk(cur)=nk(cur)-1;
        nk(newcl)=nk(newcl)+1;
        
        betaavg=(beta(cur)+beta(newcl))/2; %update beta proposal mean
        betaprop(cur)=betaavg; betaprop(newcl)=betaavg;
        betap=betaprop;
    end
end
kk1=kk;
nk1=sum(kk1,1);
dkk1=y*kk1;  %NxK

% finds beta1 using M-H;
bcand=mvnrnd(betaprop,bcovtune);  %beta is a row vector
a=mvnpdf(bcand,betaprop,bcovtune)/mvnpdf(betaprop,bcand,bcovtune);
r=cprec1/prec1;
hc=sum((bcand.^2)./(prec1*(r+nk1)))+1/zprec;
h0=sum((beta.^2)./(prec1*(r+nk1)))+1/zprec;
uc=(0.5/hc)*sum(((dkk1*(bcand./(r+nk1))')-z).^2); %final sum is across patients
u0=(0.5/h0)*sum(((dkk1*(beta./(r+nk1))')-z).^2);
lr=exp(u0 - uc) / (hc/h0)^(N/2);

if rand<lr/a
    beta1=bcand;
else
    beta1=beta;
end;


%find zprec1 using M-H;
zcand=gamrnd(ztune,zprec/ztune);
a=gampdf(zcand,ztune,zprec/ztune)/gampdf(zprec,ztune,zcand/ztune);
r=cprec1/prec1;
hc=sum((beta1.^2)./(prec1*(r+nk1)))+1/zcand;
h0=sum((beta1.^2)./(prec1*(r+nk1)))+1/zprec;
uc=(0.5/hc)*sum(((dkk1*(beta1./(r+nk1))')-z).^2);
u0=(0.5/h0)*sum(((dkk1*(beta1./(r+nk1))')-z).^2);
lr=exp(u0 - uc) / (hc/h0)^(N/2);
prior=(1/zcand)/(1/zprec);
if rand<(lr*prior)/a
    zprec1=zcand;
else
    zprec1=zprec;
end;
   

%finds p1;
ps=gamrnd(nk+alpha,ones(1,K));
p1=ps/sum(ps);

return;

end