简单易学的机器学习算法——马尔可夫链蒙特卡罗方法MCMC

import random
from scipy.stats import norm
import matplotlib.pyplot as plt

def cauchy(theta):
    y = 1.0/(1.0 + theta ** 2)
    return y

T = 50 #循环的次数
sigma = 1
thetamin = -30
thetamax = 30
theta = [0.0] * (T + 1)
theta[0] = random.uniform(thetamin,thetamax)

t = 0
while t < T:
    t += 1
    theta_star = norm.rvs(loc=theta[t - 1], scale=sigma, size=1, random_state=None) #生成数据[] 一个

    alpha = min(1,cauchy(theta_star[0])/cauchy(theta[t - 1])) #调用上边的函数

    u = random.uniform(0,1)
    if u <= alpha: #随机值小于alpha则接受
        theta[t] = theta_star[0]

    else:
        theta[t] = theta[t-1]
print(theta)