利用MCMC算法得到一般正态分布样本 R和Python

R:

snorm<-function(x){
sigma=2 
y=(1/sqrt(2*pi*sigma*sigma))*exp((-1/(2*sigma*sigma))*(x-4)^2)
#y=dnorm(x,4,2)
}

T=50
Theta<-c(NA,n=T)
Theta[1]=rnorm(1,0,1)
t=1

 

while(t<T){
t=t+1
theta_star<-rnorm(1,Theta[t-1],1)
p<-min(1,snorm(theta_star)/snorm(Theta[t-1])) 
u<-runif(1)  
if(u<p){
Theta[t]=theta_star
 }   
else{     
Theta[t]=Theta[t-1]
    }
}

Python:

import numpy as np
import math
import random
from scipy.stats import norm

def normin(x):
    sigma=2
    y=(1/np.sqrt(2*np.pi*sigma*sigma))*np.exp((-1/(2*sigma*sigma))*pow(x-4,2))
#y=norm.pdf(theta,loc=4,scale=2)
    return y

T=5000
t=0
theta1=[0.0]*(T+1)
theta1[0]=10#初始值


T=5000
t=0
theta1=[0.0]*(T+1)
theta1[0]=10#初始值