中餐馆两种采样方式:
已知条件概率
算法1:直接从联合分布中采样
N:餐厅的总人数
T:样本总数(采样的次数)
: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:样本总数(采样的次数)
: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
如何验证采样算法的正确性?
通过两方面进行比较:
- 被占桌子个数K的期望(样本均值与理论均值)
- 被占桌子个数P(K=k)的概率(样本概率与理论概率)