Python | R | MATLAB群体消息和遗传病筛选多元统计模型-CSDN博客

本文链接：https://blog.csdn.net/jiyotin/article/details/139273442

🎯要点

🎯结构化概率模型：Python问题决策影响图结构化概率模型 | 🎯多元统计模型：群体消息传播、遗传病筛选、知识工程、医疗处理决策。

🍇Python概率图图像去噪

估计贝叶斯网络条件概率分布表中的数字只是计算该事件在我们的训练数据中发生的次数。也就是说，为了估计 p(SAT=s1 | Intelligence = i1)，我们只需计算 SAT=s1 且 Intelligence = i1 的数据点在 Intelligence = i1 的总数据点中所占的比例。虽然这种方法可能看起来是临时的，但事实证明，如此获得的参数最大化了观察到的数据的可能性。

但是，对于马尔可夫网络，上述计数方法没有统计依据（因此将导致参数次优）。因此，我们需要使用更复杂的技术。大多数这些技术背后的基本思想是梯度下降 - 我们定义描述概率分布的参数，然后使用梯度下降来找到这些参数的值，以最大化观察数据的可能性。我们有了模型的参数，我们就可以在新数据上使用它们来执行推理。

推理是我们提出整个框架的原因——能够根据我们已知的信息做出预测。
推理在计算上是困难的！在某些特定类型的图上，我们可以相当有效地执行推理，但在一般图上，这是很棘手的。因此我们需要使用近似算法来权衡准确性和效率。

有几个问题我们可以通过推理来回答：

边际推理：寻找特定变量的概率分布。例如，给定一个包含变量 A、B、C 和 D 的图，其中 A 取值 1、2 和 3，求 p(A=1)、p(A=2) 和 p(A=3)。
后验推理：给定一些取值为e的观察变量v_E（E表示证据），找到一些隐藏变量v_H的后验分布p(v_H | v_E=e)。
最大后验推理：给定一些取值为 e 的观察变量 v_E，找到具有最高概率的其他变量 v_H 的设置。

接下来，我们将研究一些用于回答这些问题的流行算法，包括精确的和近似的。所有这些算法都适用于贝叶斯网络和马尔可夫网络。

使用条件概率的定义，我们可以将后验分布写为：
$p\left(v_H \mid v_E=e\right)=\frac{p\left(v_H, v_E=e\right)}{p\left(v_E=e\right)}$
让我们用一个简单的例子来看看如何计算上面的分子和分母。考虑一个具有三个变量的网络，联合分布定义如下：
$\begin{array}{|l|l|l|l|} \hline \text { A } & \text { B } & \text { C } & p ( A , B , C ) \\ \hline 0 & 0 & 0 & 0.05 \\ \hline 0 & 0 & 1 & 0.01 \\ \hline 0 & 1 & 0 & 0.10 \\ \hline 0 & 1 & 1 & 0.05 \\ \hline 1 & 0 & 0 & 0.30 \\ \hline 1 & 0 & 1 & 0.09 \\ \hline 1 & 1 & 0 & 0.25 \\ \hline 1 & 1 & 1 & 0.15 \\ \hline \end{array}$
假设我们要计算 p(A | B=1)。请注意，这意味着我们要计算值 p(A=0 | B=1) 和 p(A=1 | B=1)，它们的总和应为 1。使用上面的等式，我们可以写出
$\mid B=1)=\frac{p(A=0, B=1)}{p(B=1)}$
分子是 A = 0 且 B = 1 的概率。我们不关心 C 的值。因此我们将对 C 的所有值求和。（这来自基本概率 - p(A=0, B =1, C=0) 和 p(A=0, B=1, C=1) 是互斥事件，因此它们的并集 p(A = 0, B=1) 只是各个概率的总和。）

因此，我们将第 3 行和第 4 行相加，得到 p(A=0, B=1) = 0.15。同样，将第 7 行和第 8 行相加得出 usp(A=1, B=1) = 0.40。此外，我们可以通过对包含 B=1 的所有行（即第 3、4、7 和 8 行）求和来计算分母，得到 p(B=1) = 0.55。这给我们带来了以下结果：
$\begin{aligned} & p ( A =0 \mid B =1)=0.15 / 0.55=0.27 \\ & p ( A =1 \mid B =1)=0.40 / 0.55=0.73 \end{aligned}$
如果仔细观察上面的计算，您会注意到我们做了一些重复的计算 - 将第 3 行和第 4 行以及第 7 行和第 8 行添加两次。计算 p(B=1) 的更有效方法是简单地将值 p(A=0, B=1) 和 p(A=1, B=1) 相加。这就是变量消除的基本思想。

💦图像去噪

我们的映射推理问题可以用数学方式写成如下：
$\begin{aligned} Y^* & =\arg \max _Y P(Y \mid X) \\ & =\arg \max _Y \log P(Y \mid X) \\ & =\arg \max _Y[\log P(X, Y)-\log P(X)] \\ & =\arg \max _Y \log P(X, Y) \end{aligned}$
在这里，我们使用了最大对数似然计算中常见的一些标准简化技术。我们将使用 X 和 Y（不带下标）分别表示所有 X_ij 和 Y_ij 值的集合。

因此，我们的联合分布由下式给出：
$Y)=\frac{1}{Z} \prod_{i j} \phi\left(X_{i j}, Y_{i j}\right) \prod_{(i j, k l)} \phi\left(Y_{i j}, Y_{k l}\right)$

代码实现：

import numpy as np
from PIL import Image
import sys

def compute_log_prob_helper(Y, i, j):
  try:
    return Y[i][j]
  except IndexError:
    return 0

def compute_log_prob(X, Y, i, j, w_e, w_s, y_val):
  result = w_e * X[i][j] * y_val

  result += w_s * y_val * compute_log_prob_helper(Y, i-1, j)
  result += w_s * y_val * compute_log_prob_helper(Y, i+1, j)
  result += w_s * y_val * compute_log_prob_helper(Y, i, j-1)
  result += w_s * y_val * compute_log_prob_helper(Y, i, j+1)
  return result

def denoise_image(X, w_e, w_s):
  m, n = np.shape(X)
  Y = np.copy(X)
  max_iter = 10*m*n
  for iter in range(max_iter):

    i = np.random.randint(m)
    j = np.random.randint(n)
    log_p_neg = compute_log_prob(X, Y, i, j, w_e, w_s, -1)
    log_p_pos = compute_log_prob(X, Y, i, j, w_e, w_s, 1)

    if log_p_neg > log_p_pos:
      Y[i][j] = -1
    else:
      Y[i][j] = 1
    if iter % 100000 == 0:
      print ('Completed', iter, 'iterations out of', max_iter)
  return Y

def add_noise(orig):
  A = np.copy(orig)
  for i in range(np.size(A, 0)):
    for j in range(np.size(A, 1)):
      r = np.random.rand()
      if r < 0.1:
        A[i][j] = -A[i][j]
  return A

def convert_from_matrix_and_save(M, filename, display=False):
  M[M==-1] = 0
  M[M==1] = 255
  im = Image.fromarray(np.uint8(M))
  if display:
    im.show()
  im.save(filename)

def get_mismatched_percentage(orig_image, denoised_image):
  diff = abs(orig_image - denoised_image) / 2
  return (100.0 * np.sum(diff)) / np.size(orig_image)

def main():
  # read input and arguments
  orig_image = read_image_and_binarize(sys.argv[1])

  if len(sys.argv) > 2:
    try:
      w_e = eval(sys.argv[2])
      w_s = eval(sys.argv[3])
    except:
      print ('Run as: \npython denoise.py <input_image>\npython denoise.py <input_image> <w_e> <w_s>')
      sys.exit()
  else:
    w_e = 8
    w_s = 10


  noisy_image = add_noise(orig_image)
  denoised_image = denoise_image(noisy_image, w_e, w_s)
  print ('Percentage of mismatched pixels: ', get_mismatched_percentage(orig_image, denoised_image))

  convert_from_matrix_and_save(orig_image, 'orig_image.png', display=False)
  convert_from_matrix_and_save(noisy_image, 'noisy_image.png', display=False)
  convert_from_matrix_and_save(denoised_image, 'denoised_image.png', display=False)

if __name__ == '__main__':
  main()