5
$\begingroup$

I am trying to implement an MH algorithm to fit a power function to my data.

The power function has the following form:

$\hat{y} = a * x^b$

The data are assumed to be normally distributed around the predicted values, so:

$y = \mathcal{N}(a*x^b, \sigma)$

I chose a normal prior for a and for b and an exponential distribution for sigma.

$b = \mathcal{N}(1,0.5)$
$a = \mathcal{N}(50,20)$
$\sigma = EXP(0.01)$

Currently, my implementation doesn't work (chains for a and b are mixing very poorly) even though the code is correct (checked it with a simple linear model).

I thought about reframing the model to make it linear ($\log(y) = \log(a) + b*\log(x)$). This resolves the problem of poor mixing but I cannot use this reparameterized model because it is not possible to transform the posterior for $\sigma$ back to original scale. Also, the parameters of the linear model cannot be directly transformed to those obtained from a non-linear fit due to log-transformation (values can become outliers on log scale).

So, I am wondering,

  • if there is anything I am missing that could be the reason why the implementation doesn't work?

  • what are good (better) choices of priors here? The normal of b is obviously not a good choice since the exponent (b) cannot be negative.

More info
For all parameters, I use normal proposal distributions. This should be fine for a and b, but is it for sigma? Maybe for sigma the proposal should weight values near zero more heavily?

My desired acceptance rate is 45 %. To come close to this value, I adapt the proposal distribution widths after every 100 steps.

Code (in C)
The following code shows my implementation. The code as posted here can not be run as a single program but it should nonetheless be possible to detect any errors in the code. meanPowA, meanPowB, sdPowA and sdPowB are the prior distribution parameter.

Posterior (log)

float LogPost(float *x, float  *y, float  a, float  b, float  sigma, int N)
{

float  sumSqError = 0.0;

float logDensityA = (a - meanPowA )/sdPowA;
logDensityA = (logDensityA*logDensityA)/2;

float logDensityB = (b - meanPowB )/sdPowB;
logDensityB = (logDensityB*logDensiyB)/2;

float logDensitySigma = lambda * sigma;

for(int i = 0; i < N; i++)
{
  float resid = y[i] - (a * powf(x[i],b));
  sumSqError += resid*resid;
}

float s2 = sigma*sigma;

return -N * logf(sigma) + (-0.5/s2) * sumSqError - logDensityA - logDensityB - logDensitySigma;

}

MCMC implementation

// Metropolis Hastings
void MCMC(
float *x,
float *y,
int sampleSize,
int chainLength)
{


// standard deviations of proposals
float sigmaPropsalWidth = 0.5;
 float aProposalWidth = 0.5;
 float bProposalWidth = 0.5;

// desired acceptance rate
float desiredAcc = 0.45;


int accATot = 0;
int accA = 0;
int accBTot = 0;
int accB = 0;
int accSigmaTot = 0;
int accSigma = 0;

int chainLength = 5000;
for (int i = 1; i < chainLength; i++)
{

  a[i] = a[i - 1] + rnorm(0,aProposalWidth);

  accA = 1;

  if ((LogPost(x, y, a[i], b[i - 1], sigma[i - 1], sampleSize) -
  LogPost(x, y, a[i - 1], b[i - 1], sigma[i - 1], sampleSize)) < logf(runif(0, 1)))
  {
    a[i] = a[i - 1];
    accA = 0;
  }

  accATot += accA;


  b[i] = b[i - 1] + rnorm(0,bProposalWidth);

  accB = 1;

  if ((LogPost(x, y, a[i], b[i], sigma[i - 1], sampleSize) -
  LogPost(x, y, a[i], b[i - 1], sigma[i - 1], sampleSize)) < logf(runif(0, 1)))
  {
    b[i] = b[i - 1];
    accB = 0;
  }

  accBTot += accB;


  sigma[i] = sigma[i - 1] + rnorm(0,sigmaPropsalWidth);

  accSigma = 1;

  if (sigma[i] <= 0) {
    sigma[i] = sigma[i - 1];
    accSigma = 0;
  }
  else if ((LogPost(x, y, a[i], b[i], sigma[i],sampleSize) -
  LogPost(x, y, a[i], b[i], sigma[i - 1], sampleSize)) < logf(runif(0, 1)))
  {
    sigma[i] = sigma[i - 1];
    accSigma = 0;
  }


  accSigmaTot += accSigma;

  if ((i % 100) == 0)
  {

    sigmaProposalWidth *= ((accSigmaTot / 100.0) / desiredAcc);
    aProposalWidth *= ((accATot / 100.0) / desiredAcc);
    bProposalWidth *= ((accBTot / 100.0) / desiredAcc);

    accATot = 0;
    accBTot = 0;
    accSigmaTot = 0;

  }


  }
  }
$\endgroup$
  • 1
    $\begingroup$ Can you say more about your implementation of MH ? $\endgroup$ – peuhp Mar 4 '16 at 9:54
  • 1
    $\begingroup$ Moreover, Is your prior based on a prior belief or simply sometimes you tried among others ? $\endgroup$ – peuhp Mar 4 '16 at 10:05
  • $\begingroup$ Sorry, MH = Metropolis Hastings. $\endgroup$ – beginneR Mar 4 '16 at 10:30
  • $\begingroup$ So your proposal is 3d with no correlation between the components ? have you tryied one variable at a time MH ? Moreover, is the mixing bad even when initialising close to the expected solution ? $\endgroup$ – peuhp Mar 4 '16 at 13:26
  • 1
    $\begingroup$ You'd need to say more about your implementation. (I'd suggest also explaining why you choose to use MCMC on this problem as I suspect it could modify the kind of answers you'd get) $\endgroup$ – Glen_b Mar 5 '16 at 0:28
2
$\begingroup$

Don't bother with custom MCMC code; just use Stan. I coded your model up in rstan (the R interface to Stan) with 100 data points and it runs in seconds. Here's the R code:

model_code <- '
data {
  int<lower=0> N;
  vector[N] x;
  vector[N] y;
}
parameters {
  real a;
  real b;
  real<lower=0> sigma;
}
model {
  b ~ normal(1, 0.5);
  a  ~ normal(50, 20);
  sigma ~ exponential(0.01);
  {
    vector[N] mu;
    for (i in 1:N) mu[i] <- a * x[i] ^ b;
    y ~ normal(mu, sigma);
  }
}
'

estimate.model <- function(x,y) {
  N <- length(x)
  data <- list(N = N, x = x, y = y)
  fit <- stan(model_code = model_code, data = data)
  return(fit)
}

This returns a stanfit object containing posterior draws and various diagnostic results.

$\endgroup$
  • $\begingroup$ I should mention that my code is embedded in a larger C-project. Therefore I cannot use Stan. And unfortunately (even though Stan is written in C++), there is no C/C++ interface yet. $\endgroup$ – beginneR Mar 6 '16 at 9:12
  • 2
    $\begingroup$ I have used Stan embedded in a larger C# project. I just spawned off a separate process to run the executable generated by CmdStan. You might want to consider something similar. $\endgroup$ – Kevin S. Van Horn Mar 6 '16 at 21:29
  • $\begingroup$ Ok, interesting. And did you then import the posterior samples from Stan into C#? The "problem" in my case is that I repeatedly have to update the posterior by new incoming data. So I would have to generate a new executable each time. This would probably take some time. $\endgroup$ – beginneR Mar 7 '16 at 8:24
  • $\begingroup$ Yes, I wrote them out to a file and then read in the file. $\endgroup$ – Kevin S. Van Horn Mar 8 '16 at 15:27
  • $\begingroup$ Here is a fully reproducible example with dataset generation: gist.github.com/jsta/707bfe7b025dfb4a805b2c473a2378e7 $\endgroup$ – jsta Apr 9 '18 at 13:33

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.