stable state, it only slightly fluctuates from one iteration to the next, i.e., there is no systematic and significant increase and decrease in likelihood. We can use this as a part of the stopping criterion. Specifically, we monitor likelihood changes over a number of consecutive iterations. If the average of these changes is less than some threshold, the estimation process terminates. More robust approaches to determining the stable state are discussed elsewhere, e.g. [Ritter and Tanner 1992]. The formula for the likelihood is defined as follows.

$$f({\bf t};\psi,\phi,\theta)=\prod_{i=1:N_{K}}\left(\frac{N_{x_{i},z_{i},t_{i}}+\eta/N_{T}}{N_{x_{i},z_{i}}+\eta}\right)\tag{3}$$

To avoid a numerical precision problem in model implementation, one usually uses log likelihood log(f(t; **ψ, φ, θ**)) instead. Note that we use the strategy mentioned in [Escobar and West 1995] (Section 6) to estimate α, β and η from data.

The sampling results in the stable state are used to estimate model parameters. Again, we define Nr,z as the number of all tuples associated with resource r and topic z, with Nr, Nx,u, Nu, N**x,u,t** and Nx,z defined in a similar way. From Eq. (18) and Eq. (19) in the Appendix, the formulas for computing such parameters are as follows:

$$\phi_{r,z}={\frac{N_{r,z}+\alpha/N_{Z}}{N_{r}+\alpha}}$$
$$(4)$$

Nr + α(4)
$$\psi_{u,x}={\frac{N_{u,x}+\beta/N_{X}}{N_{u}+\beta}}$$
$$\quad({\mathfrak{H}})$$
$$(6)$$
$$\theta_{x,z,t}={\frac{N_{x,z,t}+\eta/N_{T}}{N_{x,z}+\eta}}$$

Parameter estimation via Gibbs sampling is less prone to the local maxima problem than the generic EM algorithm, as argued in [Rosen-Zvi et al. 2004]. In particular, this scheme does not estimate parameters φ, ψ, and θ directly. Rather, they are integrated out, while the hidden variables z and x are iteratively sampled during the training process. The process estimates the "posterior distribution" over possible values of φ, ψ, and θ. At a stable state, z and x are drawn from this distribution and then used to estimate φ, ψ, and θ. Consequently, these parameters are estimated from a combination of "most probable solutions", which are obtained from multiple maxima. This clearly differs from the generic EM with point estimation, which we used in our previous work [Plangprasopchok and Lerman 2007]. Specifically, the point estimation scheme estimates φ, ψ, and θ from single local maximum.

Per training iteration, the computational complexity of Gibbs sampling is more expensive than EM. This is because we need to sample hidden variables (z and x) for each data point (tuple), whereas EM only requires updating parameters. In general, the number of the data points is larger than the dimension of parameters. However, it has been reported in [Griffiths and Steyvers 2004] that to reach the same performance, Gibbs sampling requires fewer floating point operations than the other popular approaches: Variational Bayes and Expectation Propagation [Minka 2001]. Moreover, to our knowledge, there is currently no ACM Journal Name, Vol. x, No. y, zz 2010.