Fitting power-laws to data with errors in both x and y

Modeling
1

Hi, I am new to Stan/PyStan so this might be pretty dumb of me, so apologies in advance.
I am trying to fit a power law m*x^c to data generated in python with x and y both having errors.

#define our model, a line
def function(x, m, c, **kwargs):
    y = (m * x)**c
    return y

#make a function to create and plot our data
def make_data(points, m , c, xerr, yerr, seed):
    np.random.seed(seed)
    xtrue = np.linspace(1,10,points)
    ytrue = function(x = xtrue, m = m, c = c)

    xerr = xerr * np.random.randn(points)
    yerr = yerr * np.random.randn(points)
    xobs = xtrue + (xerr)
    yobs = ytrue + (yerr)
    
    plt.errorbar(xobs, yobs, xerr = xerr, yerr = yerr, fmt = 'x')
    plt.plot(xtrue, ytrue)
    plt.show()
    
    data = {'xtrue': xtrue, 'ytrue':ytrue, 'xobs':xobs, 'yobs':yobs, 'xerr':np.sqrt(abs(xerr)), 'yerr':np.sqrt(abs(yerr))}
    
    return data

data = make_data(points = 30, m = 5, c = 5, xerr = 5, yerr = 5, seed = 123)

and the Stan code is,

data {
    int<lower=1> N;
    vector[N] x;
    vector[N] y;
    vector<lower=0>[N] x_err;
    vector<lower=0>[N] y_err;
}

parameters {
    real m;
    real c;
    vector[N] x_t; // true x values
    real<lower=0> x_err_int;
    real<lower=0> y_err_int;
}

model {
    x_err ~ normal(0, x_err_int);
    y_err ~ normal(0, y_err_int);
    x ~ normal(x_t, x_err);
    for (i in 1:N)
    y[i] ~ normal((m*x_t[i])^c, y_err);
}

I get the following error

Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: normal_lpdf: Location parameter is nan, but must be finite! (in ‘unknown file name’ at line 22)

Initialization between (-2, 2) failed after 100 attempts.
Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

I googled and found that suggested setting an init_r, so I set init_r = -1 when sampling. - that told me every parameter is nan!

Any help appreciated. The same code worked exactly how I wanted it to when I fit a straight line.

2

Is the power positive? If so you need to declare it as real<lower=0> c;. And it’s probably a very good idea to give it a reasonable proper prior in the model block.

If it’s not necessarily positive you’re going to have problems if x is not positive since Stan doesn’t support complex variables.

3

Hi Michael,

Thanks, x is strictly positive so perhaps I should declare that too? I did try in one attempt to give c a prior between 0 and 15 and that did not work and I got similar errors.