(N=2)个数值10和8.5395244,在中间用线性插值插入N1=10个值,结果如下:

    10.0000000  9.7344590  9.6016885  9.4689179  9.3361474  9.2033769
     9.0706064  8.9378359  8.8050654  8.6722949  8.5395244  8.5395244

我想用 Metropolis-Hasting  算法更新插入的N1个值
假设   theta<-c(2,2,2)
Yk代表如上数据的第k 个
目标分布
  dnorm(Y[i];Y[i-1]+mu[i-1]*Dt,beta[i-1]*Dt)*dnorm(Y[i+1];Y[i]+mu[Y[i]*Dt,beta[i]*Dt])分布随机数
  建议分布
  dnorm(.;1/2*(Y[i-1]+Y[i+1]),1/2*Dt*beta[Y[i-1]])
其中 mu[i]=Y[i]+(theta[1]+theta[2]*Y[i]),beta[i]=theta[3]^2*Y[i]

请问R程序如何写?

我的程序在计算比值出现NaN,怎么改？


  使用随机游动metropolis算法产生  missing data
  目标分布
  dnorm(Y[i];Y[i-1]+mu[i-1]*Dt,beta[i-1]*Dt)*dnorm(Y[i+1];Y[i]+mu[Y[i]*Dt,beta[i]*Dt])分布随机数
  建议分布
  dnorm(.;1/2*(Y[i-1]+Y[i+1]),1/2*Dt*beta[Y[i-1]])

####################################################
# alpha<-funcion(new,Y1,x)
# 计算alpha=目标分布f（新）*建议分布g（旧）/目标分布f（旧）*建议分布g（新）

###############################################################

   alpha<-function(i,new,x,theta)
    {
       fnew<-dnorm(new,x[i-1]+(theta[1]+theta[2]*x[i-1])*Dt1,theta[3]^2*x[i-1]*Dt1)*dnorm(x[i+1],new+(theta[1]+theta[2]*new)*Dt1,theta[3]^2*new*Dt1)


        gold<-dnorm(x[i],(x[i+1]+x[i-1])/2,1/2*Dt1*theta[3]^2*x[i-1])
      
        fold<-dnorm(x[i],x[i-1]+(theta[1]+theta[2]*x[i-1])*Dt1,theta[3]^2*x[i-1]*Dt1)*dnorm(x[i+1],x[i]+(theta[1]+theta[2]*x[i]*Dt1),theta[3]^2*x[i]*Dt1)


       gnew<-dnorm(new,(x[i+1]+x[i-1])/2,1/2*Dt1*theta[3]^2*x[i-1])
  
     # return(fnew)
       #return(gold)
      # return(fold)
       #return(gnew)
      return( fnew*(gold)/((fold)*(gnew)))
}
rw.Metropolis<-function(Y1,N1,N,Dt1)

             {
   
    #Y1:update the chain
    x<-Y1                # 新更新的在x中，未更新的在Y1中
   # set i in 0:N-1;j in 1:N1,则x中不需要更新的元素序号：j+i*（11）+1

    u<-runif(N1*(N-1)+N)  # 提前取好均匀分布随机数
    k<-0
    theta<-c(2,2,2)       # 赋初值
     sigma<-theta[3]^2*(x[1])*Dt1
      for(i in 0:N-1)
        for(j in 2:（N1）)
       {
         w<-rnorm(1,(x[j+i*(N1+1)-1]+Y1[j+i*(N1+1)+1])/2,sigma)#从建议分布抽取一个数
         if(u[j+i*(N1+1)]<alpha(j+i*(N1+1），w,x,theta)) x[j+i*(N1+1)]<-w
         # else x[j+i*(N1+1）]<-x[j+i*(N1+1）]
        }
        return（x)
        }