中餐馆两种采样方式：

已知条件概率

算法1：直接从联合分布中采样

N:餐厅的总人数
T：样本总数（采样的次数）
$\alpha$：Dirichlet参数

代码1：
#算法1：直接从联合分布里采样,根据中餐馆条件概率采样，先得到z1，再根据z1得到z2.。。。。最后得到联合概率样本
#N表示中餐馆有10个人，alpha = 3表示dirichlet参数，T=50表示样本数，即迭代多少次
#verbose表示详细信息，verbose=FALSE，意思就是设置运行的时候不显示详细信息
#给定初始Z[0]=1,其他初始化为0
import numpy as np
def Draw_CRP_Direct_Sample(N = 10, alpha = 3, T = 50, VERBOSE = False):
    Z_table = np.zeros((T,N))#(T,N)用0填充的数组，数值表示每个人坐的桌子的编号
    for t in range(T):
        Z = np.zeros(N,dtype=int)#包含N个人的样本
        for i in range(N):
            if i == 0:
                Z[i] = 1#初始化Z[0] = 1,即第一个人坐的的类别
            else:
                if VERBOSE:
                    print('初始Z=',Z)#一个样本
                unique, counts = np.unique(Z, return_counts=True)
                #unique返回样本Z中去重之后的主题编号（从小到大排序），counts计数(主题权重）返回每个主题的个数
                #对于一维数组或者列表，unique函数去除其中重复的元素，并按元素由小到大返回一个新的无元素重复的元组或者列表
                #return_index=True表示返回新列表元素在旧列表中的位置，并以列表形式储存在s中
                #return_inverse=True 表示返回旧列表元素在新列表中的位置，并以列表形式储存在p中
                #return_counts=True  表示返回新列表元素在旧表中出现的次数，并以列表形式存储在counts中

                # remove the zeros unsigned tables
                if unique[0] == 0:
                    unique = np.delete(unique,0)#删除主题编号为0的元素

                if VERBOSE:
                   print("unique,counts,alpha", unique,counts,alpha)
                # added alpha to the end of the counts (weights) array
                counts = np.append(counts,alpha) #counts=[counts,alpha]
                # also the new table index will be the max of table seen so far
                unique = np.append(unique,max(unique)+1)#添加一个新的主题
                print("np.append(counts,alpha)",counts)#将counts和alpha拼接后的新counts
                print("np.append(unique,max(unique)+1)",unique)

                if VERBOSE:
                  print("sum(counts)=",sum(counts))
                #轮盘赌法 
                #u是随机产生的数，使得不同数据对应不同的数据概率，并且在整体上保留了“区域概率越大，对应数据越多”            
                u = np.random.uniform()*sum(counts)
                #np.cumsum累加求和 例子：counts=[1,2,3,4,5] np.cumsum(counts)  array([ 1,  3,  6, 10, 15], dtype=int32)
                a_counts = np.cumsum(counts)#累加求和列表，这里是对每个主题包含的个数累加求和，也可以算出每个类别的概率进行累加求和，效果是一样的
                if VERBOSE:
                    print("acounts,counts, u, (a_counts > u)",a_counts,counts, u, a_counts > u)

                # first index where accumuated sum is greater than random variable
                index =  np.argmax(a_counts > u) #返回最大值的索引
                print("index：", index)
                Z[i] = unique[index]

            if VERBOSE:
                print("最终Z=",Z)
                print("\n\n") 

        Z_table[t,:] = Z
    return Z_table

算法2:用吉布斯采样从联合样本中采样

N:餐厅的总人数
T：样本总数（采样的次数）
$\alpha$：Dirichlet参数
b:burn-in，采样完丢弃的样本数

代码2：Gibbs sampling
import numpy as np
#算法2：用吉布斯采样从联合分布中采样
#burn_in = 10代表有前10个样本会去掉
def Draw_CRP_Gibbs_Sample(N = 10, alpha = 3, T = 50, burn_in = 10, VERBOSE = False):
    Z = np.ones(N,dtype=int)#一个样本初始化为1
    Z_table = np.zeros((T,N))#所有样本初始化为0 
    for t in range(T+burn_in):
        for i in range(N):
            if VERBOSE:
                print("初始Z0=",Z)
            # remove current table assignment删除当前表赋值，处理到第几项，第几项就赋为0，Z2更新的就是Z1中为0的那一项
            Z[i] = 0

            if VERBOSE:
                print("Z1=",Z)

            unique, counts = np.unique(Z, return_counts=True)

            # remove the zeros in unassigned tables 删除未分配表中的0
            if unique[0] == 0:
                unique = np.delete(unique,0)
                counts = np.delete(counts,0)

            if VERBOSE:
                print("unique,counts,alpha", unique,counts,alpha)

            # added alpha to the end of the counts (weights) array
            counts = np.append(counts,alpha)

            # also the new table index will be the max of table seen so far
            unique = np.append(unique,max(unique)+1)

            print("np.append(counts,alpha)",counts) 
            print("np.append(unique,max(unique)+1)",unique)


            if VERBOSE:
                print("sum(counts)=",sum(counts))
            u = np.random.uniform()*sum(counts)

            a_counts = np.cumsum(counts)

            if VERBOSE:
                print("a_counts,counts, u, a_counts > u",a_counts,counts, u, a_counts > u)

            # first index where accumuated sum is greater than random variable
            index =  np.argmax(a_counts > u)

            print("index:", index)

            Z[i] = unique[index]

            if VERBOSE:
                print("Z2=",Z)
                print("\n") 

        old_table = np.unique(Z)#old_table统计出现的类别，可能不是从1开始
        print("old_table=",old_table)
        new_table = np.array(range(1,len(old_table)+1))#new_table与old_table形状大小一样，目的是将old_table中的类别变为从1开始
        print("new_table=",new_table)

        for k in range(len(old_table)):
            Z[Z == old_table[k]]=new_table[k]#将old_table中旧类别转化为新类别

        if t >= burn_in:  #舍弃前burn-in个样本
            Z_table[t-burn_in,:] = Z

        if VERBOSE:
            print("Z3=",Z)
            print("\n\n\n") 

    if VERBOSE:
        print("Z_table=",Z_table)

    return Z_table

如何验证采样算法的正确性？

通过两方面进行比较：

1. 被占桌子个数K的期望（样本均值与理论均值）
2. 被占桌子个数P(K=k)的概率（样本概率与理论概率）

参考

徐亦达：中国餐馆过程演示