mcmc: Maximum likelihood Models

In PEIP: Geophysical Inverse Theory and Optimization

Maximum likelihood Models

Description

Maximum likelihood Models

Usage

mcmc(alogprior, aloglikelihood, agenerate, alogproposal, m0, niter)

Arguments

alogprior	Name of a function that computes the log of the prior distribution.
aloglikelihood	Name of a function the computes the log of the likelihood.
agenerate	Name of a function that generates a random model from the current model using the
alogproposal	Name of a function that computes the log of the proposal distribution r(x,y).
m0	Initial model
niter	Number of iterations to perform

Value

mout	MCMC samples
mMAP	Best model found in the MCMC simulation.
accrate	Acceptance rate

Author(s)

Jonathan M. Lees<[e-mail]>

Examples

fun <-function(m,x)
{
  y=m[1]*exp(m[2]*x)+m[3]*x*exp(m[4]*x)
  return(y)
}

generate <-function( x) {  
  y=x+step*rnorm(4)
  return(y)
}

logprior <-function(m)
{
  if( (m[1]>=0) & (m[1]<=2) &
     (m[2]>=-0.9) & (m[2]<=0) &
     (m[3]>=0) & (m[3]<=2) &
     (m[4]>=-0.9) & (m[4]<=0)  )
    {
      lp=0
    }
  else
    {
      lp= -Inf
    }

  return(lp)
}
loglikelihood <-function(m)
{ 
  fvec=(y-fun(m,x))/sigma
  L=(-1/2)*sum(fvec^2)
  return(L)
}
logproposal <-function(x,y)
  {  
    LR=(-1/2)*sum((x-y)^2/step^2)
    return(LR)
  }

###  Generate the data set.
x=seq(from=1, by=0.25, to=7.0)

mtrue=c(1.0, -0.5, 1.0, -0.75)

ytrue=fun(mtrue,x)

sigma=0.01*rep(1, times= length(ytrue) )

y=ytrue+sigma*rnorm(length(ytrue) )

### set the MCMC parameters
### number of skips to reduce autocorrelation of models
skip=100
### burn-in steps
BURNIN=1000
### number of posterior distribution samples
N=4100
### MVN step size
step = 0.005*rep(1,times=4)

###  We assume flat priors here
m0 = c(0.9003,
   -0.5377,
    0.2137,
   -0.0280)

alogprior='logprior'
aloglikelihood='loglikelihood'
agenerate='generate'
alogproposal='logproposal'

### ###  initialize model at a random point on [-1,1]

###  m0=(runif(4)-0.5)*2
###  this is the matlab initialization:
m0 = c(0.9003,
   -0.5377,
    0.2137,
   -0.0280)

MM =  mcmc('logprior','loglikelihood','generate','logproposal',m0,N)

mout = MM[[1]]
mMAP= MM[[2]]
pacc= MM[[3]]

PEIP documentation built on Aug. 21, 2023, 9:10 a.m.