MCMCchain: create a MCMC chain using R.

In Hope0/StatComp18033: MLE MCMC and quicksort

Description

Markov chain Monte Carlo (MCMC) methods comprise a class of algorithms for sampling from a probability distribution.

Usage

MCMCchain(N, ntheta, datas, LM, ssigma = 1, xstart = 0, rg, tdf,
  is.log = FALSE)

Arguments

N	the length of the Markov chain
ntheta	the length of parameter
datas	the data and the prior parameter of Bayesian, also the fixed parameters of LM
LM	the likelihood function, the input should be a sample of markov chain, and data.
ssigma	the fixed parameters of transition function, such as variance of normal transition function, the acquiescent is 1
xstart	the initial value of the chain, the length is equal to ntheta, the acquiescent is 0
rg	the function of transition function. The input should be a sample of markov chain and ssigma, the output is another, one example is one dimensional normal-based on rnorm
tdf	the function of density function or probability of transition function. The input should be a sample of markov chain, ssigma, the output of function rg, one example is one dimensional normal-based on dnorm.
is.log	logical; if TRUE, the return values of LM and tdf p are given as log

Value

chain:a markov chain and accept:acceptance times

Examples

## Not run: 
X <- rnorm(10,1,2)
LMn1 <- function(theta,data){
  t <- 1
  for (i in 1:length(data)) t=t*dnorm(data[i],theta[1],2)
  return(t)
}
rnorm1 <- function(mean,sigma) return(rnorm(1,mean,sigma))
dnorm1 <- function(mchain,r,ssigma) return(dnorm(mchain,mean=r,sd=ssigma))
cha <- MCMCchain(1000,1,X,LMn1,ssigma=1,xstart=0,rg=rnorm1,tdf=dnorm1)

LMn2 <- function(theta,data){
  t <- 0
  for (i in 1:length(data)) t=t+dnorm(data[i],theta[1],theta[2],log=T)
  return(t)
}
t1 <- function(i,ssigma){
  j=i
  j[1] <- rnorm(1,i[1],ssigma[1])
  j[2] <- rexp(1,1/i[2])
  return(j)
}
d1<- function(i,j,ssigma){
  p1 <- dnorm(j[1],i[1],ssigma[1])
  p2 <- dexp(j[2],1/i[2])
  return(log(p1*p2))
}
cha2 <- MCMCchain(5000,2,X,LMn2,ssigma=c(1,2),xstart=c(0,1),rg=t1,tdf=d1,is.log=TRUE)
mean(cha2$c[1,][1000:5000])
mean(X)
mean(cha2$c[2,][1000:5000])
var(X)

## End(Not run)

Hope0/StatComp18033 documentation built on May 5, 2019, 11:06 p.m.