# Chinese Restaurant Process: Expected cardinality (number of customers) of each block (table)?

Short version of the question: The Chinese Restaurant Process defines a distribution over partitions of $$[T] := \{1, ...., T\}$$. What is the expected cardinality of the $$t$$th block, where $$t \in \{1, ..., T\}$$?

Long version of the question: Suppose we have a Chinese Restaurant Process with parameter $$\alpha$$. As background, the $$CRP(\alpha)$$ specifies a conditional distribution of the form:

$$p(z_t = k|z_1, ..., z_{t-1}) = \begin{cases} \frac{\alpha}{\alpha + t - 1} & k = K_{t-1} + 1\\ \frac{N_{t-1,k}}{\alpha + t - 1} & k \leq K_{t-1} \end{cases}$$

where $$N_{t-1, k}$$ is the integer number of customers at the $$k$$th table after the $$t-1$$ customer has been seated, and $$K_{t-1}$$ is the integer number of non-empty tables.

For each partition of $$T$$ customers, I can write a $$T$$-dimensional vector $$N_T$$ of how many customers are at each table. For instance, if $$T=3$$, the sample space of $$N_t$$ is $$\mathcal{N}_T := \{(3, 0, 0), (2, 1, 0), (1, 2, 0), (1, 1, 1) \}$$. Any of these partitions has some probability mass, given by the CRP:

$$P(N_T) = \frac{\alpha^{T-1} \prod_{k=1}^K N_{T, k}!}{(\alpha+1)_{T -1 \uparrow 1}}$$

where the denominator is the so-called rising factorial, defined as: $$x_{M \uparrow a} := \prod_{m=0}^{M-1}(x + ma)$$

My question is: what is the expected value of this random $$T$$-dimensional vector?

$$\mathbb{E}_{N_T}[N_T] = \sum_{N_T \in \mathcal{N}_T} N_T P(N_T)$$

There is some ambiguity as to what is table $$k$$, but I am going to assume that table $$k$$ is the $$k$$th table encountered in a sequential sampling as you have mentioned. That is, $$N_k = \sum_{t=1}^T 1\{z_t = k\} \implies \mathbb E [N_k] = \sum_{t=1}^T \mathbb P(z_t = k).$$ We also know that $$\mathbb P(z_t = k \mid z_{t-1},\dots,z_1) = \frac{\sum_{j=1}^{t-1} 1\{z_j = k\}}{\alpha + t - 1} 1\{k \le K_{t-1}\} + \frac{\alpha}{\alpha+t-1} 1\{ k = K_{t-1}+1 \}.$$ where $$K_{t-1} = \max\{z_1,z_2,\dots,z_{t-1}\}$$. Taking the expectation of both sides $$(\alpha+t-1) \mathbb P(z_t = k) = \sum_{j=1}^{t-1} \mathbb P(z_j = k \le K_{t-1}) + \alpha \mathbb P(K_{t-1} = k-1)$$ For $$j \le t-1$$, \begin{align*} \mathbb P(z_j = k) &= \mathbb P(z_j = k \le K_{t-1}) + \mathbb P(z_t = k, K_{t-1} < k ) \\ &= \mathbb P(z_j = k \le K_{t-1}) \end{align*} seeing that if $$K_{t-1} < k$$ by definition $$z_j < k$$ for $$j \le t-1$$.

Letting $$\beta_{tk} := \mathbb P(z_t = k)$$ and $$\gamma_{tk} = \mathbb P(K_t = k)$$, we have shown $$(\alpha + t - 1) \beta_{tk} = \sum_{j=1}^{t-1} \beta_{jk} + \alpha \gamma_{t-1,k-1}$$ for $$t \ge 2$$. By subtracting the case of $$t+1$$ from that of $$t$$, we obtain $$\beta_{t+1,k} - \beta_{tk} = \frac{\alpha}{\alpha + t}(\gamma_{t,k-1} - \gamma_{t-1,k-1})$$ hence \begin{align*} \beta_{r+1,k} &= \beta_{2k} + \sum_{t=2}^r \frac{\alpha}{\alpha + t}(\gamma_{t,k-1} - \gamma_{t-1,k-1}), \quad r \ge 2, \\ \beta_{2k} &= \frac{1}{\alpha+1} 1\{k=1\} + \frac{\alpha}{\alpha+1} 1\{k=2\}. \\ \beta_{1k}&= 1\{k=1\}. \end{align*}

The problem reduces to figuring out $$\gamma_{tk}$$. But the distribution of $$K_t$$ is known as the Chinese restaurant table distribution. Its probability mass function is given as $$\gamma_{tk} = \frac{\Gamma(\alpha)}{\Gamma(t+\alpha)} |s(t,k)| \alpha^k \cdot 1\{k \le t\}.$$ where $$s(t,k)$$ are the Stirling numbers of the first kind. Putting the pieces together \begin{align*} \mathbb E(N_k) &= \sum_{r=1}^T \beta_{rk} \\ &= \beta_{1k} + \beta_{2k} + \sum_{r=2}^{T-1} \beta_{r+1,k} \\ &= \beta_{1k} + (T-2) \beta_{2k} + \sum_{r=2}^{T-1} \sum_{t=1}^r \frac{\alpha}{\alpha + t}(\gamma_{t,k-1} - \gamma_{t-1,k-1}). \end{align*}

Here is a piece of R code verifying this calculation:

# library(gmp) # Needed for the Stirling1 numbers
set.seed(123)
sample_crp = function(n, a) {
z = vector("numeric")
z[1] = 1
for (t in 2:n) {
zmax = max(z)
freq = tabulate(z, zmax)
freq[zmax+1] = a

z[t] = sample(zmax+1, 1, T, prob = freq)
}
z
}

gam_f = function(t, k, a) {
if (k > t) return(0)
return( (gamma(a) / gamma(a + t)) * as.numeric(abs(gmp::Stirling1(t, k))) * a^k )
}
bet_f = function(r, k, a) {
if (r == 1) {
return(if (k==1) 1 else 0)
}
if (r == 2) {
if (k==1) return(1/(a+1))
if (k==2) return(a/(a+1))
return(0)
}
temp = bet_f(2, k, a)
for (t in 2:(r-1)) {
temp = temp + a*(gam_f(t, k-1, a) - gam_f(t-1, k-1, a))/(a + t)
}
temp
}
N_f = function(n, k, a) sum(sapply(1:n, function(r) bet_f(r, k, a)))

nrep = 2000
n = 20
a = 0.6
X = matrix(nrow = nrep, ncol = n)
for (rep in 1:nrep) {
X[rep, ] = sample_crp(n, a)
}

tab = apply(X, 1, FUN = function(z) tabulate(z, n))
rbind(estim=rowMeans(tab)[1:7], true=sapply(1:7, function(k) N_f(n, k, a)))


The output is

         [,1]     [,2]     [,3]      [,4]      [,5]      [,6]        [,7]
estim 12.8615 4.770500 1.694000 0.5070000 0.1355000 0.0290000 0.002500000
true  12.8750 4.772263 1.670688 0.5159738 0.1329717 0.0277614 0.004646026


Some details about the above argument:

• What happened when you took the expectation of both sides? In general, the conditional expectation $$\mathbb E(X\mid Y)$$ is a random variable (a function of $$Y$$). As a random variable, $$\mathbb E(X\mid Y)$$ itself has an expectation and we have the general "smoothing" or "tower" property: $$\mathbb E[\mathbb E(X\mid Y)] = \mathbb E(X).$$ Now take $$X = 1_{A}$$ to be the indicator of a set, i.e., $$X(\omega) = 1_A(\omega) = 1\{\omega \in A\}$$. Then, applying the above and noting that the expectation of $$1_A$$ is the probability of $$A$$, we have $$\mathbb E [ \mathbb P(A \mid Y) ] = \mathbb P(A).$$ Note that $$\mathbb P(A \mid Y)$$ is a random variable and not a deterministic quantity. This type of argument is what is going on in that step.
• Wonderful! Thank you! Yes, you correctly deduced what I mean by table $k$. – Rylan Schaeffer Jan 31 at 21:15
• When you write, "Taking the expectation of both sides," your left-hand side is $\mathbb{P}(z_t)$. Could you clarify what happened there? It seems to me you switched from the conditional distribution $\mathbb{P}(z_t | z_{<t})$ to the marginal distribution $\mathbb{P}(z_t)$, but how is that accomplished by taking an expectation? – Rylan Schaeffer Jan 31 at 21:32
• I'm also confused by something right after "We also know that". Your left hand side, the conditional distribution $\mathbb{P}(z_t∣z_{<t})$ is entirely deterministic, right? That is, given $z_{<t}$, we know exactly how much mass to place on each outcome for $z_t$. But you then say to take the expectation of both sides. How does one take the expectation of a deterministic function? Your next step suggests to me that you're treating the indicators as random variables - is this correct? If so, the left hand side is deterministic but the right hand side is random - how can this be? – Rylan Schaeffer Jan 31 at 21:47
• @RylanSchaeffer, I have added some comments to end of the answer. Please have a look and let me know if that makes it clear. $\mathbb P (z_t = k \mid z_{<t})$ is not deterministic (in general) but itself a random variable. It is a function of $z_{< t}$ as you can also see from the RHS of the equation. ($z_{<t}$ is random.) – passerby51 Jan 31 at 23:47
• One quick clarification question: in your expression for $\beta_{r+1, k}$, why is $\beta_{1, k}$ not included in the sum? – Rylan Schaeffer Feb 2 at 16:00