Assistance with an exponential model with weights

Modeling
1

I’m trying to build a Stan model for an exponential fit with the formula y = \alpha e^{(-\beta x)}
I also want to add weights to the model.

What I’ve got so far is the following model, however, it’s not working. I’m getting several errors and divergencies, but as soon as I change to the OLS formula, it all works fine. So something should be off in the model section.

data {
    int<lower=0> N;
    vector[N] x;
    vector[N] y;
    vector[N] w;
}
parameters {
    real<lower=1> alpha;
    real<lower=0> beta;
    real<lower=0> sigma;
}
model {
    y ~ normal(alpha * exp(-beta * x), sigma);
}
2

I updated my model and added some priors and everything looks good now!

data {
    int<lower=0> N;
    vector[N] x;
    vector[N] y;
    vector[N] w;
}
parameters {
    real<lower=1> alpha;
    real<lower=0> beta;
    real<lower=0> sigma;
}
model {
    alpha ~ normal(1.03, 0.03);
    beta ~ normal(0.035, 0.02);
    sigma ~ normal(0, 0.001);
    y ~ normal(alpha * exp(-beta * x), sigma);
}

Now I need to figure out how to add the weights into the model.

1 Like
3

With my fake data, your model runs without divergences.
What is the definition of weights?

y = \alpha \exp (- \beta\sum x_i w_i)?

  set.seed(1)
N <- 10 
x<- rnorm(N)
y <- rnorm(N, mean =  1*exp(-3*x),sd=1)
plot(x,y)




data <- list(
  
  x=x,
  y=y,
  N=N,
  w=x
)
library(rstan)
model <- rstan::stan_model(model_code = "
                  data {
    int<lower=0> N;
    vector[N] x;
    vector[N] y;
    vector[N] w;
}
parameters {
    real<lower=1> alpha;
    real<lower=0> beta;
    real<lower=0> sigma;
}
model {
    alpha ~ normal(0, 11);
    beta ~ normal(0, 11);
    sigma ~ normal(0, 11);
    y ~ normal(alpha * exp(-beta * x), sigma);
}
                  ")

 f<- sampling(  object = model, data = data)
  check_hmc_diagnostics(f)
  alpha <- get_posterior_mean(f)["alpha", "mean-all chains"]
  beta <- get_posterior_mean(f)["beta", "mean-all chains"]
  
  plot(x,y, xlim=c(min(x),max(x)), ylim=c(min(y),max(y)), col=2)       
  par(new=T)    
  t<-runif(1001,min(x),max(x) )
  plot(t, alpha*exp(-beta*t) ,xlim=c(min(x),max(x)), ylim=c(min(y),max(y)), col=3)

In my fake data, the parameters of truth are

\alpha = 1, \beta=3, \sigma = 1,

and the estimated posterior means are

\alpha = 1.96, \beta=2.19, \sigma = 1.13 .

Rplot01
Rplot01466×589 35.9 KB

5 Likes
4

There is a discussion about weights in this rather old google group thread:
https://groups.google.com/forum/#!topic/stan-dev/5pJdH72hoM8

2 Likes
5

I see the point that it might not be Bayesian.

In my model though, each sample has a subjective grade (RPE) on the trustworthiness of the sample(1-10, 10 being the strongest). Coupled with that, I want old samples to fade out by how many months ago the sample was taken. My weights formula is then calculated like this: w=(rpe/10)^2 \cdot e^{(-age/12)}

Is there a more Bayesian way to do this? I read elsewhere that the weights could be thought of as different sampling machine with different sigmas to the samples. How to that put into a model leaves me perplexed.

6

I do not understand enough about these things to suggest a more Bayesian way to create the effect of weights. I can see how a weight is similar to adjusting the variance allowed for the error relative to the fit but I also do not see how to translate from “I want a weight of three on this point” to some specific modification of the variance.

7

brms approaches weights in a very simple way - you multiply the likelihood of each observation by its weight. This literally means that weight of 2 is equivalent to having the same observation twice in the dataset. If you can come up with a reasonable scheme that says “This one data point is actually worth x datapoints”, this might be enough.

Alternatively, you may make sigma a function (possibly linear) of w (or even of age or rpe directly), e.g.:

parameters {
   real sigma_intercept;
   real sigma_b;
   ...
}

transformed parameters {
  vector<lower=0>[N] sigma = exp(sigma_intercept + w * sigma_b); //linear model with log link
}

Does that make sense?

1 Like