In this paper, we study probabilistic modeling of heterogeneously attributed multi-dimensional arrays. The model can manage the heterogeneity by employing an individual exponential-family distribution for each attribute of the tensor array. These entries are connected by latent variables and are shared information across the different attributes. Because a Bayesian inference for our model is intractable, we cast the EM algorithm approximated by using the Laplace method and Gaussian process. This approximation enables us to derive a predictive distribution for missing values in a consistent manner. Simulation experiments show that our method outperforms other methods such as PARAFAC and Tucker decomposition in missing-values prediction for crossnational statistics and is also applicable to discover anomalies in heterogeneous office-logging data. Keywords-tensor factorization; Bayesian probabilistic model; Gaussian process; data fusion;