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), a novel centralized inference algorithm is developed with an 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 using synthetic data and real-world data are presented to show the remarkable performance of the proposed algorithms in terms of accuracy and scalability.