function ll = DPlikelihood(sigma,theta,obs,mul)

n =  obs * mul' ;
k = sum( mul ) ;

sigma = max(sigma,eps-theta/k); % to ensure non-negativeness of probabilities

ts = theta./sigma ;

% a bit of bookkeeping to be able to handle negative sigma
num = abs(ts) + k * (sigma>0) ; 
dem = abs(ts) + 1 - k * (sigma<=0);


if length(sigma(:)) == 1
    ll=gammaln(num)-gammaln(dem)+(k-1)*log(abs(sigma)) ...
        - logpochhammer(theta+1,n-1) + mul*logpochhammer(1-sigma,obs-1)';
else
    ll=gammaln(num)-gammaln(dem)+(k-1)*log(abs(sigma)) ...
        - logpochhammer(theta+1,n-1);
    for i=1:length(mul)
        ll = ll + mul(i)*logpochhammer(1-sigma,obs(i)-1);
    end
end
    
    
function lpchh=logpochhammer(a,n)
lpchh = gammaln(a+n)-gammaln(a) ;