function P_k = negative_binomial_poisson_mixture(k,p_nbin, p_poisson)

n_nbin_max = 100*p_poisson;
n_nbin_values = 0:n_nbin_max;

P_n_nbin = poisspdf(n_nbin_values,p_poisson);

P_k_given_n_nbin = zeros(numel(k), numel(n_nbin_values));

P_k_given_n_nbin(2:end,1) = 0;
P_k_given_n_nbin(1,1) = 1;

for k_i = 1:numel(k)
    P_k_given_n_nbin(k_i,2:end) = nbinpdf(k(k_i),  n_nbin_values(2:end), p_nbin);
end

P_k = P_n_nbin*P_k_given_n_nbin';

end