You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
100 lines
2.7 KiB
100 lines
2.7 KiB
import numpy as np
|
|
|
|
PA = [0.3, 0.3, 0.2, 0.2]
|
|
PB = [0.4, 0.4, 0.1, 0.1]
|
|
PC = [0.2, 0.2, 0.3, 0.3]
|
|
PS = [[[0.2, 0.6, 0.2],
|
|
[0.1, 0.3, 0.6],
|
|
[0.05, 0.2, 0.75],
|
|
[0.01, 0.1, 0.89]],
|
|
[[0.6, 0.3, 0.1],
|
|
[0.2, 0.6, 0.2],
|
|
[0.1, 0.3, 0.6],
|
|
[0.05, 0.2, 0.75]],
|
|
[[0.75, 0.2, 0.05],
|
|
[0.6, 0.3, 0.1],
|
|
[0.2, 0.6, 0.2],
|
|
[0.1, 0.3, 0.6]],
|
|
[[0.89, 0.1, 0.01],
|
|
[0.75, 0.2, 0.05],
|
|
[0.6, 0.3, 0.1],
|
|
[0.2, 0.6, 0.2]]]
|
|
|
|
def normal(x):
|
|
s = sum(x)
|
|
res = []
|
|
for i in x:
|
|
if i == float(0):
|
|
res.append(0.0)
|
|
continue
|
|
res.append(float(i / s))
|
|
return res
|
|
|
|
def direct_cal():
|
|
res = [0.0, 0.0, 0.0]
|
|
for XA in range(4):
|
|
for XB in range(4):
|
|
for XC in range(4):
|
|
for sBC in range(3):
|
|
res[sBC] += PA[XA] * PB[XB] * PC[XC] \
|
|
* PS[XA][XB][0] \
|
|
* PS[XA][XC][1] \
|
|
* PS[XB][XC][sBC]
|
|
return normal(res) # normal(x)-x/sum(x)将计数变成概率
|
|
|
|
# 拒绝采样方法
|
|
def reject_sampling():
|
|
n = 5000
|
|
res = [0, 0, 0]
|
|
for i in range(n):
|
|
XA = np.random.choice(4, p=PA)
|
|
XB = np.random.choice(4, p=PB)
|
|
XC = np.random.choice(4, p=PC)
|
|
sAB = np.random.choice(3, p=PS[XA][XB])
|
|
sAC = np.random.choice(3, p=PS[XA][XC])
|
|
sBC = np.random.choice(3, p=PS[XB][XC])
|
|
if sAB == 0 and sAC == 1:
|
|
res[sBC] += 1
|
|
return normal(res)
|
|
|
|
# 似然加权采样方法
|
|
def likelihood_weighting():
|
|
n = 5000
|
|
res = [0, 0, 0]
|
|
for i in range(n):
|
|
w = 1
|
|
XA = np.random.choice(4, p=PA)
|
|
XB = np.random.choice(4, p=PB)
|
|
XC = np.random.choice(4, p=PC)
|
|
|
|
w = w * PS[XA][XB][0] # sAB加权
|
|
w = w * PS[XA][XC][1] # sAC加权
|
|
|
|
sBC = np.random.choice(3, p=PS[XB][XC])
|
|
res[sBC] += w
|
|
return normal(res)
|
|
|
|
# Gibbs采样方法
|
|
def Gibbs():
|
|
global PA, PB, PC
|
|
n = 4999
|
|
res = [0, 0, 0]
|
|
XA, XB, XC, sAB, sAC, sBC = 0, 0, 1, 0, 1, 1
|
|
for k in range(n):
|
|
PA = normal([PA[i] * PS[i][XB][sAB] * PS[i][XC][sAC] for i in range(4)])
|
|
XA = np.random.choice(4, p=PA)
|
|
|
|
PB = normal([PB[i] * PS[XA][i][sAB] * PS[i][XC][sBC] for i in range(4)])
|
|
XB = np.random.choice(4, p=PB)
|
|
|
|
PC = normal([PC[i] * PS[XA][i][sAC] * PS[XB][i][sBC] for i in range(4)])
|
|
XC = np.random.choice(4, p=PC)
|
|
|
|
sBC = np.random.choice(3, p=PS[XB][XC])
|
|
res[sBC] += 1
|
|
return normal(res)
|
|
|
|
print(direct_cal())
|
|
print(reject_sampling())
|
|
print(likelihood_weighting())
|
|
print(Gibbs()) |