De Blasi 2015 - Species Sampling and Unseen Species Prediction

Summary

Given a sample with $k$ distinct species, Gibbs-type priors give the exact law of the number $K_m^{(n)}$ of new species in $m$ further draws, estimators of rare-species variety and of discovery probabilities. For PY they are explicit; the frequentist Good-Toulmin estimator becomes unstable once $m>n$, the Bayesian PY estimator does not.

Main Content

Data summaries. $K_n$ distinct species, frequencies $N_n=(N_{1,n},\dots ,N_{K_n,n})$, or frequency counts $M_n=(M_{1,n},\dots ,M_{n,n})$ with $M_{i,n}$ = number of species seen exactly $i$ times; ${\sum }_i{N}_i=n$, ${\sum }_i{M}_i=K_n$, ${\sum }_{i}i{M}_i=n$. Targets: overall species variety $K_m^{(n)}=K_{n+m}-K_n$; rare species variety ${\hat{M}}_m^{(n)}(\tau )={\sum }_{i=1}^{\tau }E(M_{i,n+m}\mid X_{1:n})$ (Eq. 28); local discovery probabilities $U_{n+m,i}=P(X_{n+m+1}\in {\Delta }_{i,n+m}\mid X_{1:n})$ (Eq. 29), ${\Delta }_{i,n+m}$ the labels with frequency $i$ in the enlarged sample ($i=0$: new species).

Frequentist benchmarks. Turing: ${\check{U}}_{n,i}=(i+1)M_{i+1,n}/n$ (Eq. 30), so ${\check{U}}_{n,0}=M_{1,n}/n$ and coverage $1-M_{1,n}/n$. Good-Toulmin: ${\check{U}}_{n+m,0}=n^{-1}{\sum }_{i\ge 1}(-\lambda {)}^{i-1}i{M}_{i,n}$, ${\check{K}}_m^{(n)}={\sum }_{i\ge 1}(-1{)}^{i-1}{\lambda }^i{M}_{i,n}$ with $\lambda =m/n$ (Eq. 31); alternating sums give inadmissible values for $m>n$ roughly.

Law of new species under Gibbs-type priors (Eq. 32)

$$P(K_m^{(n)}=j\mid X_{1:n})=\frac{V_{n+m,k+j}}{V_{n,k}}\frac{\mathcal{C}(m,j;\sigma ,-n+k\sigma )}{{\sigma }^j},$$

with the non-central generalised factorial coefficient $\mathcal{C}(m,j;\sigma ,-n+k\sigma )=\frac{1}{j!}{\sum }_{r=0}^j(-1{)}^r(\genfrac{}{}{0ex}{}{j}{r}\left)\right(n-\sigma (r+k){)}_m$. So $K_n$ is sufficient for predicting the number of new species.

Pitman-Yor estimators (Eqs. 34-35, 39, 42-43)

PY$(\sigma ,\theta )$: $P(K_m^{(n)}=j\mid K_n=k)=\frac{(\theta /\sigma +k{)}_j}{(\theta +n{)}_m}\mathcal{C}(m,j;\sigma ,-n+k\sigma )$; squared-loss estimator of overall variety

$${\hat{K}}_m^{(n)}=(k+\frac{\theta }{\sigma })[\frac{(\theta +n+\sigma {)}_m}{(\theta +n{)}_m}-1];$$

Bayesian Turing analogue ${\hat{U}}_{n,i}=\frac{i-\sigma }{\theta +n}{M}_{i,n}$ ($i\ge 1$) and new-species discovery probability at step $n+m+1$

$${\hat{U}}_{n+m,0}=\frac{\theta +k\sigma }{\theta +n}\frac{(\theta +n+\sigma {)}_m}{(\theta +n+1{)}_m}.$$

General Gibbs-type: ${\hat{U}}_{n,0}=V_{n+1,k+1}/V_{n,k}$ and ${\hat{U}}_{n,i}=\frac{V_{n+1,k}}{V_{n,k}}(i-\sigma )M_{i,n}$ (Eq. 39).

Asymptotics (Eq. 36). Conditionally on $X_{1:n}$, $K_m^{(n)}/m^{\sigma }\to Z_{n,k}$ a.s. as $m\to \infty$, with $Z_{n,k}\stackrel{d}{=}{B}_{k+\theta /\sigma ,n/\sigma -k}{U}_{(\theta +n)/\sigma }$ ($B$ beta, $U_q$ with density $f_{U_q}(u)=\frac{\Gamma (q\sigma +1)}{\sigma \Gamma (q+1)}{u}^{q-1-1/\sigma }{f}_{\sigma }(u^{-1/\sigma })$, $f_{\sigma }$ the positive $\sigma$-stable density). Setting $n=k=0$ recovers the $\sigma$-diversity. Credible intervals via simulation of $Z_{n,k}$. Rare species (Sec. 4.2). ${\hat{M}}_m^{(n)}(\tau )={\hat{N}}_m^{(n)}(\tau )+{\hat{O}}_m^{(n)}(\tau )$ (new species with frequency $\le \tau$ plus old species updated to frequency $\le \tau$); explicit sufficient statistic $(K_n,M_{1,n},\dots ,M_{\tau ,n})$ for old species (Eq. 37), $K_n$ for new (Eq. 38). PY: ${\hat{N}}_{i,m}=\left(\genfrac{}{}{0ex}{}{m}{i}\right)(1-\sigma {)}_{i-1}(\theta +k\sigma )\frac{(\theta +n+\sigma {)}_{m-i}}{(\theta +n{)}_m}$; $M_{i,n+m}/m^{\sigma }\to \frac{\sigma (1-\sigma {)}_{i-1}}{i!}{Z}_{n,k}$. For discovery of species with frequency $i$ in the enlarged sample: Eq. 40 and, for PY, Eq. 43.

Application (Table 2, Fig. 5). Naegleria gruberi ESTs from two cDNA libraries (Susko and Roger 2004). Frequency counts $M_i$, $n$ and distinct genes $j$:

Library	$M_1$	$M_2$	$M_3$	$M_4$	$M_5$	$M_6$	$M_7$	$M_8$	$M_9$	$M_{10}$	$M_{11}$
Aerobic	346	57	19	12	9	5	4	2	4	5	4
Anaerobic	491	72	30	9	13	5	3	1	2	0	1

Library	$M_{12}$	$M_{13}$	$M_{14}$	$M_{15}$	$M_{16}$	$M_{17}$	$M_{18}$	$M_{27}$	$M_{55}$	$j$	$n$
Aerobic	1	0	0	0	1	1	1	1	1	473	959
Anaerobic	0	1	3	0	0	0	0	0	0	631	969

Empirical-Bayes $(\sigma ,\theta )=(0.66,155.5)$; compare ${\check{U}}_{n+m,0}$ (Good-Toulmin) with ${\hat{U}}_{n+m,0}$ (PY) for $m=\mathrm{1..2000}$. The Good-Toulmin estimator is erratic for $m$ large relative to $n$; both agree for small $m$; the anaerobic library shows clearly higher gene diversity. One can choose the additional sample size $m$ as the largest value with discovery probability above a cost-based threshold. Footnote: estimators (40) and (43) correct a minor inaccuracy in the earlier discovery-probability estimators of the authors’ reference [19].

Examples

Recomputing the aerobic library (runnable)

$n=959$, $k=473$, $\sigma =0.66$, $\theta =155.5$. Turing gives ${\check{U}}_{n,0}=346/959=0.361$; PY gives ${\hat{U}}_{n,0}=(155.5+0.66\cdot 473)/(155.5+959)=0.420$. As $m$ grows the Good-Toulmin estimator oscillates and diverges while PY decays smoothly:

$m$	GT ${\check{U}}_{n+m,0}$	PY ${\hat{U}}_{n+m,0}$	GT ${\check{K}}_m^{(n)}$	PY ${\hat{K}}_m^{(n)}$
0	0.361	0.420	0	0
500	0.311	0.370	167	196
959	0.366	0.340	305	359
1500	$1.8\times 10^9$	0.314	$4.8\times 10^{10}$	535
2000	$9.9\times 10^{15}$	0.296	$3.6\times 10^{17}$	688

import numpy as np
from scipy.special import gammaln
M = {1:346,2:57,3:19,4:12,5:9,6:5,7:4,8:2,9:4,10:5,11:4,12:1,16:1,17:1,18:1,27:1,55:1}
n = sum(i*c for i, c in M.items()); k = sum(M.values())     # 959, 473
sig, th = 0.66, 155.5
lr = lambda x, m: gammaln(x + m) - gammaln(x)
U_py = lambda m: (th + k*sig)/(th + n) * np.exp(lr(th + n + sig, m) - lr(th + n + 1, m))
K_py = lambda m: (k + th/sig) * (np.exp(lr(th + n + sig, m) - lr(th + n, m)) - 1)
U_gt = lambda m: sum((-m/n)**(i-1) * i * c for i, c in M.items()) / n
K_gt = lambda m: sum((-1)**(i-1) * (m/n)**i * c for i, c in M.items())
for m in (0, 500, 959, 1500, 2000):
    print(m, round(U_gt(m), 3), round(U_py(m), 3), round(K_gt(m), 1), round(K_py(m), 1))

Relevance to marketing measurement

Replace species by SKUs, creative variants, keywords or customers. ${\hat{K}}_m^{(n)}$ answers “how many new distinct items will appear in the next $m$ transactions”, ${\hat{U}}_{n+m,0}$ “what fraction of the next transactions will be a first-seen item”. The exchangeability and the discrete-support assumptions are strong; check with held-out time splits.

Connections

Teh 2010 - Blackwell-MacQueen Urn and Chinese Restaurant Process, De Blasi 2015 - EPPF and Predictive Structure. Related vault area on unseen counts: Pareto-NBD Model.

See Also

De Blasi 2015 - Gibbs-Type Priors - Overview