标签：采样 Metropolis stats mus python random np sigmas
文章目录
马氏链
MH采样
代码
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
# 正态分布
x_=np.linspace(-20,20,100)
y_=stats.norm.pdf(x_,0,5)# 正态分布
# y_=stats.expon(scale=1).pdf(x_)# 指数分布
# 采样数10000
Samp_Num=10000
result=[]
init=1
result.append(init)
# p=lambda r:stats.expon(scale=1).pdf(r)# 指数分布
p=lambda r:stats.norm.pdf(r,0,5) # 正态分布
q=lambda v:stats.norm.rvs(loc = v,scale = 2, size = 1)
for i in range(Samp_Num):
y=q(result[i])# 从分布q(y|x_t)中采样
alpha=min(1,p(y)/p(result[i]))# 接受概率(简化)
u=np.random.rand(1)# 从uniform(0,1)中采样
if u<alpha:
result.append(y[0])# 接受
else:
result.append(result[i])# 拒绝
if i%1000==0:
print(i)
plt.hist(result, 50, density=1, facecolor='blue', alpha=0.5)
plt.plot(x,raw_y)
plt.show()
Gibbs采样
代码
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
def p_x_given_y(y, mus, sigmas):
mu = mus[0] + sigmas[1, 0] / sigmas[0, 0] * (y - mus[1])
sigma = sigmas[0, 0] - sigmas[1, 0] / sigmas[1, 1] * sigmas[1, 0]
return np.random.normal(mu, sigma)
def p_y_given_x(x, mus, sigmas):
mu = mus[1] + sigmas[0, 1] / sigmas[1, 1] * (x - mus[0])
sigma = sigmas[1, 1] - sigmas[0, 1] / sigmas[0, 0] * sigmas[0, 1]
return np.random.normal(mu, sigma)
def gibbs_sampling(mus, sigmas, iter=10000):
samples = np.zeros((iter, 2))
y = np.random.rand() * 10
for i in range(iter):
x = p_x_given_y(y, mus, sigmas)
y = p_y_given_x(x, mus, sigmas)
samples[i, :] = [x, y]
return samples
mus = np.array([5, 5])
sigmas = np.array([[1, .9], [.9, 1]])
x,y = np.random.multivariate_normal(mus, sigmas, int(1e5)).T
sns.jointplot(x,y,kind='kde')
samples = gibbs_sampling(mus, sigmas)
sns.jointplot(samples[:, 0], samples[:, 1])
标签：采样,Metropolis,stats,mus,python,random,np,sigmas 来源： https://blog.csdn.net/NP_hard/article/details/120441853
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。