Fitting an exponential model

I want to fit the following model.
y = e^(b0 + (b1 * X))

Using python and curve_fit from scipy, I get this result wich is correct.

import numpy as np
from scipy.optimize import curve_fit

x = np.array([1,2])
y = np.array([6,7])

def f(t,b0,b1):
    return np.exp(b0 + (b1 * t))
p, cov = curve_fit(f, x, y)
#parameters
print(p) 

I tried to implement this with pystan, but the values that I get are very different from the solution found with scipy.
(beta must be equal to something like -0.00180411)

model_code = """
data {
  int<lower=0> N;
  vector[N] y;
  vector[N] x;
}
parameters {
  real beta0;
  real beta1;
  real<lower=0> sigma;
}
transformed parameters {
  real m;
  real b0;
  real b1;
  b0 = exp(beta0);
  b1 = exp(beta1);
  for (i in 1:N) 
    m = b0+ (b1 * x[i]);
}
model {
  // priors
  beta0 ~ normal(0,10); 
  beta1 ~ normal(0, 10);     
  sigma ~ cauchy(0,5); 
  // likelihood
  for (n in 1:N)
   y[n] ~ normal(m, sigma);
}
"""
dat = {'N': 2,
           'y': y,
           'x': x }

sm = pystan.StanModel(model_code=model_code)
fit = sm.sampling(data=dat, iter=50000, chains=1)
fit

How can I obtain the solution of scipy with stan ?
Thank you.

1 Like

Without having looked too long at it (so maybe there is more)… but you declare m as real such that it is scalar. Thus m will be set to the Nth value at the end of that loop. You probably meant to declare m as vector and assign values in the loop to ith element.

However, you can also just do (after declaring m as vector of length N)

m = b0 + b1 * x;

and then you should also do instead of the for loop in the model block again a vectorized statement like

y ~ normal(m, sigma)

this should bring you much closer to what you are expecting.

How to predict forecasted values for this exponential model?

1 Like