4
$\begingroup$

Context

I looked literally everywhere but I couldn't find a complete implementation of the Random Walk Metropolis-Hastings algorithm using the log scale. By log scale I mean that we are working with the logarithm of the target distribution (which would usually be a posterior). For simplicity here I have a minimal working example in R but I would be happy also with a python implementation.

Problems

I do not understand the following things:

  1. If I use the log scale, should I also use a log-proposal?
  2. If I use the log scale, does it matter if the function that evaluates the log of my target distribution is just proportional to the log target?
  3. How should I implement it?

Minimal Working Example Context

Bayesian Logistic Regression (with 0-1 labels) log-posterior is given by

logistic regression

I want to sample from this posterior using Random-Walk Metropolis-Hastings algorithm. According to posts such as this and this and this blog post it is better to use the log posterior.

Minimal Working Example - Random-Walk Metropolis-Hastings function in R

library(MASS)
rwmh_log <- function(start, niter, logtarget){
    # Set current z to the initial point and calculate its log target to save computations
    z  <- start    # It's a column vector
    pz <- logtarget(start)
    # Generate matrix containing the samples. Initialize first sample with the starting value
    samples <- matrix(0, nrow=niter, ncol=nrow(start))
    # Generate uniform random numbers in advance, to save computation. Take logarithm?
    log_u <- log(runif(niter))
    # Proposal is a multivariate standard normal distribution. Generate samples and
    # later on use linearity property of Gaussian distribution
    normal_shift <- mvrnorm(n=niter, mu=c(0,0,0), Sigma=diag(nrow(start)))
    for (i in 2:niter){
        # Sample a candidate
        candidate <- z + normal_shift[i, ]
        # calculate log target of candidate and store it in case it gets accepted
        p_candidate <- logtarget(candidate)
        # use decision rule explained in blog posts
        if (log_u[i] <= pz - p_candidate){
            # Accept!
            z  <- candidate
            pz <- p_candidate
        }
        # Finally add the sample to our matrix of samples
        samples[i, ] <- z
    }
    return(samples)
}

Logistic Regression Data Generation to run the minimal working example

set.seed(123)
# Number of observations, mean and variance-covariance matrix for class 0
n1 <- 100
m1 <- c(6, 6)
s1 <- matrix(c(1, 0, 0, 10), nrow=2, ncol=2)
# Number of observations, mean and variance-covariance matrix for class 1
n2 <- 100
m2 <- c(-1, 1)
s2 <- matrix(c(1, 0, 0, 10), nrow=2, ncol=2)
# Generate explanatory data by sampling bivariate normal distributions
class1 <- mvrnorm(n1, m1, s1)
class2 <- mvrnorm(n2, m2, s2)
# Generate class labels. First n1 are of class 0, last n2 are of class 1
y  <- c(rep(0, n1), rep(1, n2))
X  <- rbind(class1, class2)

Plot of the generated data set (just for fun)

library(ggplot2)
data <- data.frame(X, y)
ggplot(data=data, aes(x=X1, y=X2, color=as_factor(y))) + 
    geom_point() + 
    theme(plot.title=element_text(hjust=0.5, size=20)) + 
    labs(color="Class", title="Linearly Separable Dataset")

dataset

Logistic regression log posterior to run the minimal working example

# NOTICE THAT X NEEDS A COLUMN OF 1s FOR THE BIAS
X <- cbind(1, X)
log_posterior_unnormalized <- function(beta){
    log_prior      <- -0.5*sum(beta^2)
    log_likelihood <- -sum(log(1 + exp((1 - 2*y) * (X %*% beta))))
    return(log_prior + log_likelihood)
}

Running the Minimal Working Example

start    <- matrix(0, nrow=3, ncol=1) # 3 because of the bias
niter    <- 2000
samples  <- rwmh_log(start=start, niter=niter, logtarget=log_posterior_unnormalized)

Results of the Minimal Working Example

It errors, saying 

Error in if (log_u[i] <= pz - p_candidate) { : 
  missing value where TRUE/FALSE needed

which basically means that both pz and p_candidate become -Inf. What is happening??

Running this algorithm on the command line we can see that it does something for the first 102 iterations and breaks at the 103th. The samples generated up to that point are pasted below. You can see how they blow up.

> samples[1:102, ]
              [,1]       [,2]       [,3]
  [1,]  0.00000000  0.0000000  0.0000000
  [2,]  1.41003898  0.7678739 -1.1686514
  [3,]  1.41003898  0.7678739 -1.1686514
  [4,]  1.41003898  0.7678739 -1.1686514
  [5,]  1.12738511  1.3976872 -0.4979555
  [6,]  1.12738511  1.3976872 -0.4979555
  [7,]  1.12738511  1.3976872 -0.4979555
  [8,]  1.68208878  2.7622693  0.2584510
  [9,]  1.68208878  2.7622693  0.2584510
 [10,]  1.68208878  2.7622693  0.2584510
 [11,]  1.68208878  2.7622693  0.2584510
 [12,]  1.48830009  3.6458290  0.5262860
 [13,]  3.13261741  4.2619841  1.1795437
 [14,]  3.13261741  4.2619841  1.1795437
 [15,]  2.02763486  5.2880473  0.7658672
 [16,]  2.02763486  5.2880473  0.7658672
 [17,]  1.96977976  6.9251877  0.6729261
 [18,] -0.68979341  7.6298268  1.1032108
 [19,]  0.24805648  7.7550307  1.6386097
 [20,]  0.40805008  8.9963608  1.0833313
 [21,]  0.10764974  9.1156086  2.8628342
 [22,] -0.38268725 10.3344574  3.1492587
 [23,] -0.38268725 10.3344574  3.1492587
 [24,]  0.85171606 10.1059391  4.4215254
 [25,]  0.30011080 11.1485574  3.7030592
 [26,]  0.61490763 12.6665140  3.2527206
 [27,] -0.89593039 12.3662036  5.6501731
 [28,] -0.53032703 13.1415187  5.6613023
 [29,] -2.97562829 13.3948230  7.2948707
 [30,] -2.97562829 13.3948230  7.2948707
 [31,] -2.67296962 14.8767398  7.1043539
 [32,] -2.34678595 15.0634645  7.4827778
 [33,] -2.28767265 15.1808296  7.7828163
 [34,] -2.28767265 15.1808296  7.7828163
 [35,] -2.22941911 15.9571394  7.8020756
 [36,] -2.22941911 15.9571394  7.8020756
 [37,] -2.22941911 15.9571394  7.8020756
 [38,] -3.61326092 17.1933051  8.8868507
 [39,] -3.61326092 17.1933051  8.8868507
 [40,] -3.61326092 17.1933051  8.8868507
 [41,] -3.61326092 17.1933051  8.8868507
 [42,] -1.81956128 17.3817482  9.3416200
 [43,] -3.16750568 17.4750150 10.0015226
 [44,] -0.79642112 17.6399588  9.8016328
 [45,] -1.32371016 18.3821248  9.1565188
 [46,] -1.32371016 18.3821248  9.1565188
 [47,] -2.34696374 19.0425533  9.5953375
 [48,] -2.51424008 18.6796395 10.4786403
 [49,] -2.51424008 18.6796395 10.4786403
 [50,] -2.51424008 18.6796395 10.4786403
 [51,] -3.36106702 20.5262499 11.9090427
 [52,] -2.63987051 20.3130838 12.9556715
 [53,] -2.63987051 20.3130838 12.9556715
 [54,] -0.43401179 19.8467996 13.6708499
 [55,] -1.07079158 19.0419739 14.5880249
 [56,] -1.07079158 19.0419739 14.5880249
 [57,] -1.07079158 19.0419739 14.5880249
 [58,] -1.65936235 21.3411496 14.1030373
 [59,] -0.98841697 21.9482180 14.3336541
 [60,] -0.09823962 23.8173624 14.0384963
 [61,]  0.02893740 24.2862548 14.9104612
 [62,]  0.02893740 24.2862548 14.9104612
 [63,] -0.72302884 24.7495870 15.4289650
 [64,]  0.13822989 26.9019761 15.0382800
 [65,]  0.13822989 26.9019761 15.0382800
 [66,]  0.04357883 26.8745540 16.2482905
 [67,] -0.88960758 27.9014297 16.9891905
 [68,] -0.68093651 29.0824771 18.7134528
 [69,] -0.43058184 29.8667399 18.7786067
 [70,] -1.29463115 31.0491849 19.9036095
 [71,] -1.24440787 29.5035870 21.8790285
 [72,] -1.24440787 29.5035870 21.8790285
 [73,] -1.55557884 31.5517166 20.5560774
 [74,] -1.99585038 32.3840898 20.3167258
 [75,] -1.99585038 32.3840898 20.3167258
 [76,] -1.99585038 32.3840898 20.3167258
 [77,] -1.89981602 33.6591331 22.0290308
 [78,] -1.89981602 33.6591331 22.0290308
 [79,] -2.44680063 34.6652981 22.4020355
 [80,] -2.44680063 34.6652981 22.4020355
 [81,] -1.60876820 35.8310586 22.4224862
 [82,] -1.45422863 36.5512324 22.7365438
 [83,] -1.18233421 35.5138059 24.0647585
 [84,] -0.24570214 37.0486802 24.1860769
 [85,]  0.03168349 37.4168439 24.8989192
 [86,]  0.03168349 37.4168439 24.8989192
 [87,]  0.03168349 37.4168439 24.8989192
 [88,]  0.03168349 37.4168439 24.8989192
 [89,]  0.24215122 37.5586010 26.5258005
 [90,]  0.24215122 37.5586010 26.5258005
 [91,]  0.24215122 37.5586010 26.5258005
 [92,]  0.58862001 39.3607620 27.9298507
 [93,] -0.27926896 40.1595364 29.2239346
 [94,] -0.27926896 40.1595364 29.2239346
 [95,] -0.27926896 40.1595364 29.2239346
 [96,] -0.27926896 40.1595364 29.2239346
 [97,] -0.17760335 40.9974928 29.4057818
 [98,]  0.43759118 42.1254915 29.5706227
 [99,] -0.13618441 42.4222319 29.9347374
[100,]  0.93734292 42.7479366 30.4868951
[101,]  0.93734292 42.7479366 30.4868951
[102,]  0.93734292 42.7479366 30.4868951

Looking at the value of the log_posterior_unnormalized on such samples we can see that it does indeed blow up to infinity.

> apply(samples[1:107, ], 1, log_posterior_unnormalized)
  [1]   -138.6294   -310.8512   -310.8512   -310.8512   -814.2805   -814.2805   -814.2805  -2196.6386  -2196.6386
 [10]  -2196.6386  -2196.6386  -2969.5265  -3877.5393  -3877.5393  -4299.3816  -4299.3816  -5418.8160  -6098.4219
 [19]  -6520.9842  -7102.3063  -8222.7219  -9253.5215  -9253.5215  -9933.2233 -10196.3848 -11024.7336 -12248.1868
 [28] -12816.9682 -13964.8656 -13964.8656 -14897.0547 -15285.9802 -15566.7016 -15566.7016 -16129.8982 -16129.8982
 [37] -16129.8982 -17662.7135 -17662.7135 -17662.7135 -17662.7135 -18160.6234 -18610.1692 -18687.7953 -18765.9209
 [46] -18765.9209 -19482.5643 -19800.5161 -19800.5161 -19800.5161 -22021.5857 -22595.2153 -22595.2153 -22838.3542
 [55] -22862.6686 -22862.6686 -22862.6686 -24129.0525 -24740.6719 -25905.4325 -26827.0686 -26827.0686 -27470.8678
 [64] -28779.8978 -28779.8978 -29565.7594 -30755.0321 -32768.3554 -33383.2801 -34951.9200 -35192.1701 -35192.1701
 [73] -35742.4804 -36158.4290 -36158.4290 -36158.4290 -38238.1653 -38238.1653 -39191.6603 -39191.6603 -40079.6649
 [82] -40818.7965 -40987.1734 -42216.0336 -42980.1015 -42980.1015 -42980.1015 -42980.1015 -44203.8974 -44203.8974
 [91] -44203.8974 -46483.0324 -47911.4547 -47911.4547 -47911.4547 -47911.4547 -48647.6199 -49606.0710 -50046.5853
[100]        -Inf        -Inf        -Inf

This makes me think that maybe the decision rule is not working properly? Maybe it's accepting samples that it should not accept?

Strange Result: Should I sample the NEGATIVE log posterior?

Somehow if I use the Random Walk Metropolis-Hastings algorithm to sample the negative log posterior I get somehow sensible results? How is that possible? Here is my code to see what's going on.

start    <- matrix(0, nrow=3, ncol=1) # 3 because of the bias
niter    <- 100000
samples  <- rwmh_log(start=start, niter=niter, logtarget=log_posterior_unnormalized)

If then we use the following plotting code

samplesdf <- data.frame(samples) %>% mutate(rn=row_number())
trace1 <- ggplot(data=samplesdf, aes(x=rn, y=X1)) + geom_line()
trace2 <- ggplot(data=samplesdf, aes(x=rn, y=X2)) + geom_line()
trace3 <- ggplot(data=samplesdf, aes(x=rn, y=X3)) + geom_line()
grid.arrange(trace1, trace2, trace3, ncol=1)

we get

traceplot

and with the following plotting code we get the histograms for each parameter coordinate

hist1 <- ggplot(data=samplesdf, aes(x=X1, stat(density))) + 
  geom_histogram(binwidth=0.05, alpha=0.5, fill="turquoise1", color="turquoise4")
hist2 <- ggplot(data=samplesdf, aes(x=X2, stat(density))) + 
  geom_histogram(binwidth=0.05, alpha=0.5, fill="turquoise1", color="turquoise4")
hist3 <- ggplot(data=samplesdf, aes(x=X3, stat(density))) + 
  geom_histogram(binwidth=0.05, alpha=0.5, fill="turquoise1", color="turquoise4")
grid.arrange(hist1, hist2, hist3, ncol=1)

giving

histograms

These results can be made even better by starting at the map estimate and using the inverse of the approximate hessian matrix as the variance-covariance matrix for the normal proposal.

set.seed(123)
# Number of observations, mean and variance-covariance matrix for class 0
n1 <- 100
m1 <- c(6, 6)
s1 <- matrix(c(1, 0, 0, 10), nrow=2, ncol=2)
# Number of observations, mean and variance-covariance matrix for class 1
n2 <- 100
m2 <- c(-1, 1)
s2 <- matrix(c(1, 0, 0, 10), nrow=2, ncol=2)
# Generate explanatory data by sampling bivariate normal distributions
class1 <- mvrnorm(n1, m1, s1)
class2 <- mvrnorm(n2, m2, s2)
# Generate class labels. First n1 are of class 0, last n2 are of class 1
y  <- c(rep(0, n1), rep(1, n2))
X  <- rbind(class1, class2)
# Add colum of 1s for the bias
X <- cbind(1, X)
# log posterior
log_posterior_unnormalized <- function(beta){
  log_prior      <- -0.5*sum(beta^2)
  log_likelihood <- -sum(log(1 + exp((1 - 2*y) * (X %*% beta))))
  return(log_prior + log_likelihood)
}
# define negative log posterior to be minimized by optim function
logtarget <- function(x) -log_posterior_unnormalized(x)
# optimize it to find hessian and starting point
r <- optim(c(0,0,0), logtarget, method="BFGS", hessian=TRUE)
# Start at the MAP estimate
niter <- 100000
start <- r$par
z <- start
pz <- logtarget(z)
samples <- matrix(0, nrow=niter, ncol=3)
# Use the inverse of the approximate hessian matrix for our normal proposal
vcov <- solve(r$hessian) 
normal_shift <- mvrnorm(n=niter, mu=c(0,0,0), Sigma=vcov)
samples[1, ] <- start
log_u <- log(runif(niter))
for (i in 2:niter){
  # Sample a candidate
  candidate <- z + normal_shift[i, ]
  # calculate log target of candidate and store it in case it gets accepted
  p_candidate <- logtarget(candidate)
  # use decision rule explained in blog posts
  if (log_u[i] <= pz - p_candidate){
    # Accept!
    z  <- candidate
    pz <- p_candidate
  }
  # Finally add the sample to our matrix of samples
  samples[i, ] <- z
}


samplesdf <- data.frame(samples) %>% mutate(rn=row_number())
trace1 <- ggplot(data=samplesdf, aes(x=rn, y=X1)) + geom_line()
trace2 <- ggplot(data=samplesdf, aes(x=rn, y=X2)) + geom_line()
trace3 <- ggplot(data=samplesdf, aes(x=rn, y=X3)) + geom_line()
grid.arrange(trace1, trace2, trace3, ncol=1)

hist1 <- ggplot(data=samplesdf, aes(x=X1, stat(density))) + 
  geom_histogram(binwidth=0.05, alpha=0.5, fill="turquoise1", color="turquoise4")
hist2 <- ggplot(data=samplesdf, aes(x=X2, stat(density))) + 
  geom_histogram(binwidth=0.05, alpha=0.5, fill="turquoise1", color="turquoise4")
hist3 <- ggplot(data=samplesdf, aes(x=X3, stat(density))) + 
  geom_histogram(binwidth=0.05, alpha=0.5, fill="turquoise1", color="turquoise4")
grid.arrange(hist1, hist2, hist3, ncol=1)

giving the following two plots

traceplot2

histograms2

Improve this question
$\endgroup$
2
  • $\begingroup$ @Xi’an Thank you! I don’t think it does, the problem is not so much that the probability values are too small, but that the chain is choosing to go to areas of the parameter space having progressively smaller and smaller mass, thus spending all of its time in areas with negligible probability $\endgroup$
    – Euler_Salter
    Jan 13, 2020 at 10:00
  • $\begingroup$ This other question also wonders about the log-transform. $\endgroup$
    – Xi'an
    Jan 13, 2020 at 13:58

1 Answer 1

Reset to default
5
$\begingroup$

Short answer: this is a code error.

  1. If I use the log scale, should I also use a log-proposal?

The probability of acceptance is always $$\min\left\{1,\dfrac{\pi(x^\text{new})}{\pi(x^\text{old})}\times\dfrac{q(x^\text{old}|x^\text{new})}{q(x^\text{new}|x^\text{old})} \right\}$$ which can also be written as $$\min\left\{1,\exp[\ln\pi(x^\text{new})-\ln\pi(x^\text{old})]\times\dfrac{q(x^\text{old}|x^\text{new})}{q(x^\text{new}|x^\text{old})} \right\}$$ or as $$\min\left\{1,\exp[\ln\pi(x^\text{new})-\ln\pi(x^\text{old})+\ln q(x^\text{old}|x^\text{new})-\ln q(x^\text{new}|x^\text{old})] \right\}$$ or yet as $$\min\left\{1,\dfrac{\pi(x^\text{new})}{\pi(x^\text{old})}\times\exp[\ln q(x^\text{old}|x^\text{new})-\ln q(x^\text{new}|x^\text{old})] \right\}$$ In short one can use any bijective transform of the target or proposal, provided the ratio remains the same.

  1. If I use the log scale, does it matter if the function that evaluates the log of my target distribution is just proportional to the log target?

Proportionality to the target, $\pi(x)\propto\tilde\pi(x)$ say, does not pass to the log-target: $\ln\pi(x)\,{\!\not\mathrel{\propto\,}}\ln\tilde\pi(x)$

Should I sample the NEGATIVE log posterior?

No but the R code takes the reverse of $$\exp[\ln\pi(x^\text{new})-\ln\pi(x^\text{old})]$$

    candidate <- z + normal_shift[i, ]
    # calculate log target of candidate and store it in case it gets accepted
    p_candidate <- logtarget(candidate)
    # use decision rule explained in blog posts
    if (log_u[i] <= pz - p_candidate){

which explains for the diverging behaviour.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.