mcmc logistic

Author: Aragorn

code_in_notes/MCMC_independent.py

@@ -12,7 +12,7 @@
 domain=lambda x:(x[0]>=-1)*(x[1]>=-1)*(x[0]<=1)*(x[1]<=1)
 q=lambda y: 1/(2*np.pi)*np.exp(-1*y[0]**2/2-y[1]**2/2)
 h=lambda x: domain(x)
-h2=lambda x: np.sin(x[0])**2+np.log(abs(1+x[1]))
+h2=lambda x: np.sin(10*x[0])**2+np.log(abs(1+x[1]*10))
 #从均匀分布中采样
 def sample_q():
     return (nprd.normal(),nprd.normal())

code_in_notes/MCMC_logistic.py

@@ -0,0 +1,82 @@
+#!/usr/bin/python3
+## file: MonteCarlo.py
+
+import numpy as np
+from numpy import random as nprd
+
+##设定参数
+M =20000
+#真值
+beta_0=1
+beta_1=1
+beta_2=-1
+#样本量
+N=200
+#Logistic函数
+Logistic=lambda x: 1.0/(1+np.exp(-1*x))
+
+##产生数据
+def gen_logit(N):
+    Data=[]
+    for n in range(N):
+        x1=nprd.normal()*1.414+1
+        x2=nprd.chisquare(2)
+        d_star=beta_0+beta_1*x1+beta_2*x2
+        p_star=Logistic(d_star)
+        d=(1 if nprd.uniform()<p_star else 0)
+        Data.append((d,x1,x2))
+    return Data
+
+#计算接受率
+def rho(beta_x,beta_y,Data):
+    log_post_pai_x=(-1*beta_x[0]**2-beta_x[1]**2-beta_x[2]**2)/2
+    log_post_pai_y=(-1*beta_y[0]**2-beta_y[1]**2-beta_y[2]**2)/2
+    log_ratio=log_post_pai_y-log_post_pai_x
+    for data in Data:
+        w_b_x=beta_x[0]+data[1]*beta_x[1]+data[2]*beta_x[2]
+        w_b_y=beta_y[0]+data[1]*beta_y[1]+data[2]*beta_y[2]
+        F_b_x=Logistic(w_b_x)
+        F_b_y=Logistic(w_b_y)
+        log_post_pai_x=np.log((F_b_x if data[0]==1 else 1-F_b_x))
+        log_post_pai_y=np.log((F_b_y if data[0]==1 else 1-F_b_y))
+        log_ratio+=(log_post_pai_y-log_post_pai_x)
+    return min(1,np.exp(log_ratio))
+
+#随机游走
+def q_sampler(beta):
+    return [b+nprd.normal(0,0.1) for b in beta]
+
+##独立的MCMC算法，输入：
+##    N_samples  : 抽样次数
+##      rho(x,y,Data)   : 计算接受率
+##   q_sampler(x): 给定x，从q中抽样的函数
+##      x0       : 初始值
+##     data      : 数据
+def MH_RW(N_samples, rho, q_sampler, x0, data):
+    X=[]
+    x=x0
+    for i in range(N_samples):
+        y=q_sampler(x)
+        if nprd.uniform()<=rho(x,y,data):
+            X.append(y)
+            x=y
+        else:
+            X.append(x)
+    return X
+
+## 从后验抽样：
+data=gen_logit(N)
+beta_post=MH_RW(M, rho, q_sampler, [0,0,0], data)
+beta0_post=[b[0] for b in beta_post]
+sub_beta0=beta0_post[int(M*0.2):]
+beta1_post=[b[1] for b in beta_post]
+sub_beta1=beta1_post[int(M*0.2):]
+beta2_post=[b[2] for b in beta_post]
+sub_beta2=beta2_post[int(M*0.2):]
+#后验均值
+mean_beta0=np.mean(sub_beta0)
+mean_beta1=np.mean(sub_beta1)
+mean_beta2=np.mean(sub_beta2)
+print("Mean beta0=",mean_beta0)
+print("Mean beta1=",mean_beta1)
+print("Mean beta2=",mean_beta2)

code_in_notes/MCMC_random_walk.py

@@ -6,15 +6,13 @@
 
 ##设定参数
 M =10000
-pai_cons=90*np.sqrt(2)/(2*np.pi)
-pai=lambda x: pai_cons* \
-    np.exp(-90*(x[0]-0.5)**2-45*(x[1]+0.1)**2)
+pai=lambda x: np.exp(-90*(x[0]-0.5)**2-45*(x[1]+0.1)**2)
 domain=lambda x:(x[0]>=-1)*(x[1]>=-1)*(x[0]<=1)*(x[1]<=1)
 h=lambda x: domain(x)
-h2=lambda x: np.sin(x[0])**2+np.log(abs(1+x[1]))
+h2=lambda x: np.sin(10*x[0])**2+np.log(abs(1+x[1]*10))
 #从均匀分布中采样
 def sample_q(x):
-    return (x[0]+0.3*nprd.normal(),x[1]+0.3*nprd.normal())
+    return (x[0]+0.2*nprd.normal(),x[1]+0.2*nprd.normal())
 
 
 ##独立的MCMC算法，输入：
@@ -48,3 +46,4 @@ def MH_RW(N_samples, pai, q_sampler, x0):
 subH2=H2[int(M*0.2):]
 integral2=np.pi/(90*np.sqrt(2))*np.mean(subH2)
 print("Intgral2=",integral2)
+