function samples = DPsample_sigma_theta(data,param)
% Gibbs sampler for Dirichlet process model parametet sigma and theta.

% Ole Winther Feb 2009

% calculate count stats from data
n =  data.obs * data.mul' ;     % total number of counts
k = sum( data.mul ) ;           % total number of species observed     

% burn-in and collect samples 
samples.sigma=zeros(param.samples,1);
samples.theta=zeros(param.samples,1);

theta = param.theta(ceil(length(param.theta)*rand)); % select random theta

for t=1:param.burnin+param.samples

    inddef = (theta + k * param.sigma) > 0 ; % select values that give well-defined  probabilities
    sigma = sigmaSample(param.sigma(inddef),theta,data.obs,data.mul,param.sigmaprior(inddef));

    inddef = (param.theta + k * sigma) > 0 ; % select values that give well-defined  probabilities
    theta = thetaSample(sigma,param.theta(inddef),n,k,param.thetaprior(inddef));
    
    if t>param.burnin
        samples.sigma(t-param.burnin)=sigma;
        samples.theta(t-param.burnin)=theta;
    end

end
 
%************** Sample from {theta | sigma} **************/
 

function W=thetaWeights(sigma,theta,n,k,thetap)
ts = theta./abs(sigma) ;
num = max(0,ts + k * (sigma>0)) ; 
dem = max(0,ts + 1 - k * (sigma<=0));

W = log(thetap) + gammaln(num)-gammaln(dem) - logpochhammer(theta+1,n-1) ;
Wmax = max(W);
W = exp(W-Wmax);
W = W / sum(W);
W = cumulate(W);

function theta=thetaSample(sigma,theta,n,k,thetap);

W=thetaWeights(sigma,theta,n,k,thetap);
c1 = find(rand<W,1);
theta = theta(c1);

function Wc = cumulate(W)

Wc = W;
for i=2:length(W)
    Wc(i) = Wc(i) + Wc(i-1) ;
end

%*************** Sample from {sigma | theta} ***************/

function W=sigmaWeights(sigma,theta,obs,mul,sigmap)
n =  obs * mul' ;
k = sum( mul ) ;

ts = theta./abs(sigma) ;
num = max(0,ts + k * (sigma>0)) ; 
dem = max(0,ts + 1 - k * (sigma<=0));

W=gammaln(num)-gammaln(dem)+(k-1)*log(abs(sigma))+log(sigmap);
for i=1:length(mul)
    W=W+mul(i)*logpochhammer(1-sigma,obs(i)-1);
end
Wmax = max(W);
W = exp(W-Wmax);
W = W / sum(W);
W = cumulate(W);

function sigma=sigmaSample(sigma,theta,obs,mul,sigmap);

W=sigmaWeights(sigma,theta,obs,mul,sigmap);
c1 = find(rand<W,1);
sigma = sigma(c1);

% /* ------------------------- Help functions -------------------------- */

function lpchh=logpochhammer(a,n)
lpchh = gammaln(a+n)-gammaln(a) ;
