1

I'm trying to implement a random walk Metropolis-Hastings algorithm in R. I have used the the self-defined functions logit and invlogit to apply and undo the logit function. I have also used the normal distribution to add random noise. Given these two facts, once you use the inverse transformation on the transformed parameter + random noise, the proposal distribution is no longer symmetric, which is why I then apply the corrective term log(yt*(1 - yt)) - log(xt*(1 - xt)) to the acceptance probability.

My problem is that there seems to be two ways to implement this algorithm in R. If both of these ways were equal, then, as I understand it, I should be getting equal values for both when calculating acceptanceRate. However, this is not the case, which leads me to believe that one implementation is flawed (has bugs) and the other doesn't.

However, two other possibilities are that (1) both ways are incorrect or (2) both ways are correct and I'm misunderstanding something. I'm new to R coding, and I still can't understand why these deviations exist in the value of acceptanceRate between the two implementations.

Note: my specific problem of interest is why I'm getting different value for acceptanceRate between the two implementations.

Implementation 1

log.posterior <- function(p) (12+p)*log(p) + (9-p)*log(1-p)

B <- 10000           ## number of realisations we want to have
chain <- rep(0, B+1)  ## vector to hold realisations
chain[1] <- 0.5       ## initial value
num.accept <- 0       ## keep track on how often we accept proposals
for(i in 1:B){
  xt <- chain[i] ## current point
  logit <- function(p) log(p/(1-p))
  invlogit <- function(lo) 1/(1 + exp(-lo))
  yt <- invlogit(rnorm(1, mean = logit(xt), sd = 0.45))      ## proposal
  lapt <- log.posterior(yt) - log.posterior(xt) + log(yt*(1 - yt)) - log(xt*(1 - xt))   ## acceptance probability on the log scale)
  if( runif(1) <= exp(lapt) ){
    chain[i+1] <- yt    ## accept proposal if runif(1) is less or equal to the acceptance probility
    num.accept <- num.accept + 1 ## proposal was accepted
  }else
    chain[i+1] <- xt    ## reject proposal
}

acceptanceRate <- num.accept/B

See how implementation 1 uses yt <- invlogit(rnorm(1, mean = logit(xt), sd = 0.45))? Everything is accumulated and done together.

Implementation 2

log.posterior <- function(p) (12+p)*log(p) + (9-p)*log(1-p)

B <- 10000           ## number of realisations we want to have
chain <- rep(0, B+1)  ## vector to hold realisations
chain[1] <- 0.5       ## initial value
num.accept <- 0       ## keep track on how often we accept proposals
for(i in 1:B){
  xt <- chain[i] ## current point
  logit <- function(p) log(p/(1-p))
  xt <- logit(xt)
  yt <- xt + rnorm(1, mean = 0, sd = 0.45)      ## proposal
  invlogit <- function(lo) 1/(1 + exp(-lo))
  xt <- invlogit(xt)
  yt <- invlogit(yt)
  lapt <- log.posterior(yt) - log.posterior(xt) + log(yt*(1 - yt)) - log(xt*(1 - xt))   ## acceptance probability on the log scale)
  if( runif(1) <= exp(lapt) ){
    chain[i+1] <- yt    ## accept proposal if runif(1) is less or equal to the acceptance probility
    num.accept <- num.accept + 1 ## proposal was accepted
  }else
    chain[i+1] <- xt    ## reject proposal
}

acceptanceRate <- num.accept/B

Notice that implementation 2 breaks everything down into separate pieces are proceeds sequentially.

halfer
  • 19,824
  • 17
  • 99
  • 186
The Pointer
  • 2,226
  • 7
  • 22
  • 50
  • 1
    too broad.......... – Mitch Wheat Apr 19 '18 at 12:27
  • @MitchWheat But it's about a specific value? Did you read the entire thing? – The Pointer Apr 19 '18 at 12:28
  • 1
    This is off topic. Voted to close. – ABCD Apr 19 '18 at 12:41
  • @SmallChess is this better suited towards stats.stackexchange? I thought, since it was a code-specific problem, it would be suitable for here? – The Pointer Apr 19 '18 at 12:42
  • @Suren Oh? It seemed to me like it was. I will take it elsewhere if people think that’s more suitable? – The Pointer Apr 19 '18 at 12:46
  • I take it back. Maybe it is the implementation. – kangaroo_cliff Apr 19 '18 at 12:51
  • @Suren Yes, I was quite sure it’s my inexperience with R and lack of familiarity with how it deals with variables and such that may be causing the problems. – The Pointer Apr 19 '18 at 12:52
  • In implementation 1, you seem to use the `invlogit` before you declare it. These function doesn't seem to be declared inside the `for` loop anyway. – kangaroo_cliff Apr 19 '18 at 12:54
  • You might want to check out https://codereview.stackexchange.com/ – kath Apr 19 '18 at 12:55
  • @Suren Thanks. That was a transcription error on my part. Fixed now. – The Pointer Apr 19 '18 at 12:59
  • @kath Thanks. Will take a look. – The Pointer Apr 19 '18 at 13:00
  • @ThePointer Apologies for hijacking this post, but I noticed you took down your recent rstan question. We were exchanging comments earlier. If still an issue, please undelete, and I can post an answer to clarify. – Maurits Evers Apr 24 '18 at 12:18
  • @MauritsEvers That's fine. I think I managed to fix the issue, but I've encountered new issues. If you have any time, I would, of course, appreciate any assistance you could provide: https://stats.stackexchange.com/questions/342470/modelling-parameter-r-max-limits-i-1-dots-10-p-i-min-limits-i – The Pointer Apr 24 '18 at 13:37
  • 1
    @ThePointer Can I suggest that you re-open the previous question. I'd like to show a few optimisations, that might also help with your CV question. – Maurits Evers Apr 24 '18 at 13:46
  • 1
    @MauritsEvers I have undeleted it. Someone voted to close it, so I thought it best to delete it. Thank you for the assistance. – The Pointer Apr 24 '18 at 14:05
  • @MauritsEvers any chance you could look at my similar quesiton [here](https://stackoverflow.com/questions/59709258/random-walk-metropolis-hastings-implementation-in-r-using-log-scale)? – Euler_Salter Jan 13 '20 at 02:20
  • 1
    @Euler_Salter I had a look at your question. I haven't touched this subject since the time I posted this, so I don't remember anywhere near enough to help you. The best I could do is give your question an upvote. I wish I could be of more help. – The Pointer Jan 13 '20 at 02:43
  • 1
    @ThePointer don't worry! Thank you again :) – Euler_Salter Jan 13 '20 at 02:49

2 Answers2

2

Apparently, OP compared the two functions that depend on random number generators without setting the seed (set.seed).

I don't see what is wrong with it. For a small chain, I get the same outcome.

log.posterior <- function(p) (12+p)*log(p) + (9-p)*log(1-p)
invlogit <- function(lo) 1/(1 + exp(-lo))
logit <- function(p) log(p/(1-p))

set.seed(1)

B <- 100           ## number of realisations we want to have
chain <- rep(0, B+1)  ## vector to hold realisations
chain[1] <- 0.5       ## initial value
num.accept <- 0       ## keep track on how often we accept proposals
for(i in 1:B){
  xt <- chain[i] ## current point

  xt <- logit(xt)
  yt <- xt + rnorm(1, mean = 0, sd = 0.45)      ## proposal

  xt <- invlogit(xt)
  yt <- invlogit(yt)
  lapt <- log.posterior(yt) - log.posterior(xt) + log(yt*(1 - yt)) - log(xt*(1 - xt))   ## acceptance probability on the log scale)
  if( runif(1) <= exp(lapt) ){
    chain[i+1] <- yt    ## accept proposal if runif(1) is less or equal to the acceptance probility
    num.accept <- num.accept + 1 ## proposal was accepted
  }else
    chain[i+1] <- xt    ## reject proposal
}

acceptanceRate <- num.accept/B

# acceptanceRate 
# [1] 0.69

# chain[30:40]  
# [1] 0.7674114 0.6612332 0.5867199 0.5867199 0.5744098 0.6033942 0.5359917  [8] 0.5359917 0.5359917 0.6040635 0.6040635

log.posterior <- function(p) (12+p)*log(p) + (9-p)*log(1-p)
logit <- function(p) log(p/(1-p))
invlogit <- function(lo) 1/(1 + exp(-lo))

set.seed(1)
B <- 100           ## number of realisations we want to have
chain <- rep(0, B+1)  ## vector to hold realisations
chain[1] <- 0.5       ## initial value
num.accept <- 0       ## keep track on how often we accept proposals
for(i in 1:B){
  xt <- chain[i] ## current point

  yt <- invlogit(rnorm(1, mean = logit(xt), sd = 0.45))      ## proposal
  lapt <- log.posterior(yt) - log.posterior(xt) + log(yt*(1 - yt)) - log(xt*(1 - xt))   ## acceptance probability on the log scale)
  if( runif(1) <= exp(lapt) ){
    chain[i+1] <- yt    ## accept proposal if runif(1) is less or equal to the acceptance probility
    num.accept <- num.accept + 1 ## proposal was accepted
  }else
    chain[i+1] <- xt    ## reject proposal
}

acceptanceRate <- num.accept/B

# acceptanceRate 
# [1] 0.69

# chain[30:40]  
# [1] 0.7674114 0.6612332 0.5867199 0.5867199 0.5744098 0.6033942 0.5359917  [8] 0.5359917 0.5359917 0.6040635 0.6040635
kangaroo_cliff
  • 6,067
  • 3
  • 29
  • 42
2

The issue is that you are using random numbers, so to have reproducible results you need to use set.seed before running your algorithms.
I pulled out the definition of the function from the for-loop and used set.seed. I got the same results in both cases:

log.posterior <- function(p) (12+p)*log(p) + (9-p)*log(1-p)
logit <- function(p) log(p/(1-p))
invlogit <- function(lo) 1/(1 + exp(-lo))

1st implementation

set.seed(42)
B <- 10000  ## number of realisations we want to have
chain <- rep(0, B+1)  ## vector to hold realisations
chain[1] <- 0.5       ## initial value
num.accept <- 0       ## keep track on how often we accept proposals
for(i in 1:B){
  xt <- chain[i] ## current point
  yt <- invlogit(rnorm(1, mean = logit(xt), sd = 0.45))      ## proposal
  lapt <- log.posterior(yt) - log.posterior(xt) + log(yt*(1 - yt)) - log(xt*(1 - xt))   ## acceptance probability on the log scale)
  if( runif(1) <= exp(lapt) ){
    chain[i+1] <- yt    ## accept proposal if runif(1) is less or equal to the acceptance probility
    num.accept <- num.accept + 1 ## proposal was accepted
  }else
    chain[i+1] <- xt    ## reject proposal
}

acceptanceRate1 <- num.accept/B

rm(B, chain, num.accept, i, lapt, xt, yt)

2nd implementation

set.seed(42)
B <- 10000           ## number of realisations we want to have
chain <- rep(0, B+1)  ## vector to hold realisations
chain[1] <- 0.5       ## initial value
num.accept <- 0       ## keep track on how often we accept proposals

for(i in 1:B){
  xt <- chain[i] ## current point
  xt <- logit(xt)
  yt <- xt + rnorm(1, mean = 0, sd = 0.45)      ## proposal
  xt <- invlogit(xt)
  yt <- invlogit(yt)
  lapt <- log.posterior(yt) - log.posterior(xt) + log(yt*(1 - yt)) - log(xt*(1 - xt))   ## acceptance probability on the log scale)
  if( runif(1) <= exp(lapt) ){
    chain[i+1] <- yt    ## accept proposal if runif(1) is less or equal to the acceptance probility
    num.accept <- num.accept + 1 ## proposal was accepted
  }else
    chain[i+1] <- xt    ## reject proposal
}

acceptanceRate2 <- num.accept/B

acceptanceRate1
# [1] 0.7029
acceptanceRate2
# [1] 0.7029
kath
  • 7,624
  • 17
  • 32
  • Ahhh, ok, I see what's going on now! Thank you for taking the time to review! – The Pointer Apr 19 '18 at 13:05
  • any chance you could look at my question [here](https://stackoverflow.com/questions/59709258/random-walk-metropolis-hastings-implementation-in-r-using-log-scale) ? – Euler_Salter Jan 13 '20 at 02:21