Nonnegative tensor canonical polyadic decomposition (CPD) has found wide-spread applications in various signal processing tasks. However, the implementation of most existing algorithms needs the knowledge of tensor rank, which is difficult to acquire. To address this issue, by interpreting the nonnegative CPD problem using probability density functions (pdfs), the problem is recasted as probabilistic inference with integrated feature of automatic rank determination. Furthermore, to scale the inference algorithm to massive data, its implementation under modern distributed computing architecture is investigated, giving rise to a distributed probabilistic nonnegative tensor CPD algorithm. Numerical studies are presented to show the remarkable performance of the proposed algorithms in terms of accuracy and scalability.