Metropolis Hastings для модели линейной регрессии

Я пытаюсь реализовать алгоритм Метрополиса-Гастингса для простой линейной регрессии на C (без использования других библиотек (boost, Eigen и т.д.) и без двумерных массивов)*. Для лучшего тестирования кода/оценки графиков трассировки я переписал код для R (см. ниже), сохранив как можно больше C-кода.

К сожалению, цепочки не сходятся. мне интересно, если

  1. ошибка в самой реализации?
  2. "просто" неудачный выбор дистрибутивов предложений?

Предполагая последнее, я думаю о том, как найти хорошие параметры распределения предложений (в настоящее время я выбрал произвольные значения), чтобы алгоритм работал. Даже с тремя параметрами, как в данном случае, найти подходящие параметры довольно сложно. Как обычно решают эту проблему, если, скажем, выборка Гиббса не является альтернативой?

*Я хочу использовать этот код для Cuda

#### posterior distribution
logPostDensity <- function(x, y, a, b, s2, N)
{
sumSqError = 0.0
for(i in 1:N)
{
  sumSqError = sumSqError + (y[i] - (a + b*x[i]))^2
}
return(((-(N/2)+1) * log(s2)) + ((-0.5/s2) * sumSqError))

}

# x = x values
# y = actual datapoints
# N = sample size
# m = length of chain
# sigmaProp = uniform proposal for sigma squared
# paramAProp = uniform proposal for intercept
# paramBProp = uniform proposal for slope

mcmcSampling <- function(x,y,N,m,sigmaProp,paramAProp,paramBProp)
{

  paramsA = vector("numeric",length=m) # intercept
  paramsB = vector("numeric",length=m) # slope
  s2 = vector("numeric",length=m) # sigma squared

  paramsA[1] = 0
  paramsB[1] = 0
  s2[1] = 1

  for(i in 2:m)
  {

    paramsA[i] = paramsA[i-1] + runif(1,-paramAProp,paramAProp)

    if((logPostDensity(x,y,paramsA[i],paramsB[i],s2[i-1],N)
        - logPostDensity(x,y,paramsA[i-1],paramsB[i-1],s2[i-1],N))
       < log(runif(1)))
    {
      paramsA[i] = paramsA[i-1]
    }

    paramsB[i] = paramsB[i-1] + runif(1,-paramBProp,paramBProp)

    if((logPostDensity(x,y,paramsA[i],paramsB[i],s2[i-1],N)
        - logPostDensity(x,y,paramsA[i-1],paramsB[i-1],s2[i-1],N))
       < log(runif(1)))
    {
      paramsB[i] = paramsB[i-1]
    }

    s2[i] = s2[i-1] + runif(1,-sigmaProp,sigmaProp)

    if((s2[i] < 0) || (logPostDensity(x,y,paramsA[i],paramsB[i],s2[i],N)
                       - logPostDensity(x,y,paramsA[i],paramsB[i],s2[i-1],N))
       < log(runif(1)))
    {
      s2[i] = s2[i-1]
    }


  }

  res = data.frame(paramsA,paramsB,s2)
  return(res)
}


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

set.seed(321)
x <- runif(100)
y <- 2 + 5*x + rnorm(100)

summary(lm(y~x))


df <- mcmcSampling(x,y,10,5000,0.05,0.05,0.05)


par(mfrow=c(3,1))
plot(df$paramsA[3000:5000],type="l",main="intercept")
plot(df$paramsB[3000:5000],type="l",main="slope")
plot(df$s2[3000:5000],type="l",main="sigma")

c r linear-regression mcmc
person beginneR    schedule 17.11.2015    source источник
comment
Во-первых, для формального байесовского анализа вам нужны априорные распределения по «a», «b» и «s2», которых у вас сейчас нет. Но если вы предполагаете неинформативный априор, это не должно иметь большого значения. Я бы попробовал адаптивный Metropolis Hastings, где вы динамически настраиваете распределение вашего предложения (Jump) (то есть значения переменных prop) на основе текущей скорости принятия (которую вам придется отслеживать с помощью трех счетчиков). Я думаю, что причина, по которой вы не получаете конвергенцию, заключается в том, что вы не исследуете полное пространство параметров.   -  person Alexey Shiklomanov    schedule 17.11.2015
comment
Спасибо. Как вы предложили, я предполагаю неинформативные априоры. Я слышал об адаптивной выборке MH, но не знаю, как ее реализовать.   -  person beginneR    schedule 18.11.2015
comment
Это довольно просто. По сути, просто отслеживайте свой уровень принятия по блоку образцов, возьмите соотношение между этим и целевым уровнем принятия (обычно это 0,44), а затем умножьте ширину распространения вашего предложения на это соотношение. Например, для заданного блока из 100 образцов, если вы принимаете 88 из них, ваш фактический коэффициент принятия составляет 0,88, поэтому соотношение составляет 0,88/0,44 = 2, поэтому вы должны умножить значение Prop для этого параметра на 2. Вот пример реализации.   -  person Alexey Shiklomanov    schedule 18.11.2015
comment
Спасибо, я посмотрю, что я могу сделать.   -  person beginneR    schedule 18.11.2015

Ответы (1)


arrow_upward
0
arrow_downward

В секции перехвата была одна ошибка (paramsA). Все остальное было в порядке. Я реализовал то, что предложил Алексей в своих комментариях. Вот решение:

pow <- function(x,y)
{
  return(x^y)
}


#### posterior distribution
posteriorDistribution <- function(x, y, a, b,s2,N)
{
sumSqError <- 0.0
for(i in 1:N)
{
  sumSqError <- sumSqError + pow(y[i] - (a + b*x[i]),2)
}
return((-((N/2)+1) * log(s2)) + ((-0.5/s2) * sumSqError))

}

# x <- x values
# y <- actual datapoints
# N <- sample size
# m <- length of chain
# sigmaProposalWidth <- width of uniform proposal dist for sigma squared
# paramAProposalWidth <- width of uniform proposal dist for intercept
# paramBProposalWidth <- width of uniform proposal dist for slope

mcmcSampling <- function(x,y,N,m,sigmaProposalWidth,paramAProposalWidth,paramBProposalWidth)
{

  desiredAcc <- 0.44

  paramsA <- vector("numeric",length=m) # intercept
  paramsB <- vector("numeric",length=m) # slope
  s2 <- vector("numeric",length=m) # sigma squared

  paramsA[1] <- 0
  paramsB[1] <- 0
  s2[1] <- 1

  accATot <- 0
  accBTot <- 0
  accS2Tot <- 0

  for(i in 2:m)
  {
    paramsA[i] <- paramsA[i-1] + runif(1,-paramAProposalWidth,paramAProposalWidth)
    accA <- 1
    if((posteriorDistribution(x,y,paramsA[i],paramsB[i-1],s2[i-1],N) - 
        posteriorDistribution(x,y,paramsA[i-1],paramsB[i-1],s2[i-1],N)) < log(runif(1)))
    {
      paramsA[i] <- paramsA[i-1]
      accA <- 0
    }


    accATot <- accATot + accA

    paramsB[i] <- paramsB[i-1] + runif(1,-paramBProposalWidth,paramBProposalWidth)
    accB <- 1
    if((posteriorDistribution(x,y,paramsA[i],paramsB[i],s2[i-1],N) - 
        posteriorDistribution(x,y,paramsA[i-1],paramsB[i-1],s2[i-1],N)) < log(runif(1)))
    {
      paramsB[i] <- paramsB[i-1]
      accB <- 0
    }

    accBTot <- accBTot + accB

    s2[i] <- s2[i-1] + runif(1,-sigmaProposalWidth,sigmaProposalWidth)
    accS2 <- 1

    if((s2[i] < 0) || (posteriorDistribution(x,y,paramsA[i],paramsB[i],s2[i],N) - 
                       posteriorDistribution(x,y,paramsA[i],paramsB[i],s2[i-1],N)) < log(runif(1)))
    {
      s2[i] <- s2[i-1]
      accS2 <- 0
    }

    accS2Tot <- accS2Tot + accS2

    if(i%%100==0)
    {

      paramAProposalWidth <- paramAProposalWidth * ((accATot/100)/desiredAcc)
      paramBProposalWidth <- paramBProposalWidth * ((accBTot/100)/desiredAcc)
      sigmaProposalWidth <- sigmaProposalWidth * ((accS2Tot/100)/desiredAcc)

      accATot <-  0
      accBTot <-  0 
      accS2Tot <-  0

    }


  }
    res <- data.frame(paramsA,paramsB,s2)
    return(res)

}
person beginneR    schedule 18.11.2015