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#初始值