Fitting power function to data

Question

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;

  }


  }
  }

Answer 1

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.