# Fitting power function to data

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.

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;

}

}
}

• Can you say more about your implementation of MH ? – peuhp Mar 4 '16 at 9:54
• Moreover, Is your prior based on a prior belief or simply sometimes you tried among others ? – peuhp Mar 4 '16 at 10:05
• Sorry, MH = Metropolis Hastings. – beginneR Mar 4 '16 at 10:30
• 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 ? – peuhp Mar 4 '16 at 13:26
• 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) – Glen_b Mar 5 '16 at 0:28

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.

• 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. – beginneR Mar 6 '16 at 9:12
• 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. – Kevin S. Van Horn Mar 6 '16 at 21:29
• 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. – beginneR Mar 7 '16 at 8:24
• Yes, I wrote them out to a file and then read in the file. – Kevin S. Van Horn Mar 8 '16 at 15:27
• Here is a fully reproducible example with dataset generation: gist.github.com/jsta/707bfe7b025dfb4a805b2c473a2378e7 – jsta Apr 9 '18 at 13:33