Weighting observations and measurement uncertainty in bayes

Question

I am working on using MCMC (via STAN) to estimate model parameters for a bunch of observations with measurement uncertainty. I'm having problems with weighting each observation, and have reduced the problem to something more simple which still has the pathology I'm trying to fix -- a Pareto distribution.

Let's say I generate data in R in the following way:

rpareto = function(n,location,shape){location/runif(n)^(1/shape)}

N = 1000
std = 0.2
underlying = rpareto(N,location=10,shape=2)
logx = log(underlying) + rnorm(N,0,std)    

So that I have Pareto variates (with $y>10$ and a slope of -2), with an uncertainty of 0.2 in logscale. I want to recover the location and shape.

The following STAN code will do this fine:

data {
  int<lower=0> N;
  vector[N] logx;
  real<lower=0> std;
}
parameters {
  real<lower=1,upper=4> logy_min;
  real<lower=0,upper=4> alpha;
  vector<lower=logy_min,upper=10>[N] log_xtrue;
}
model {
  exp(log_xtrue) ~ pareto(exp(logy_min), alpha);
  logx ~ wnormal(log_xtrue, std);
}

Here, I have the variates in log-space because that makes more sense in my final application (it could go either way here).

However, in my final application it makes sense to weight each observation so that higher variates are upweighted (ie. we flatten out the contribution at different scales). This can mathematically be easily done by giving each observation $x_{true}$ a weight of $w_i=x^s_{i,true}/\bar{x^s}_{true}$. In this case, the whole power-law is tilted by the scale $s$, and we need to estimate $\alpha+s$ in the model.

Conceptually, what's happening is that we're replacing every observation by $x^s_{true}$ of the same observation, as if we had data with a PDF that was proportional to $\alpha+s$. Typically I want to use $s\sim 1$.

In the case that there are no uncertainties on the observations, this works fine, that is, the following model (used on data with no uncertainties), works perfectly to give alpha=2:

functions {
  real wpareto_log(vector y,real y_min,real alpha,vector weight){
    return sum((log(alpha)+alpha*log(y_min)-(alpha+1)*log(y)) .*weight);
  }
}
data {
  int<lower=0> N;
  vector[N] logx;
  real s; //scale 
}
transformed data {
  vector[N] weight;
  for (i in 1:N) weight[i] <- pow(exp(logx[i]),s);
  weight <- weight/mean(weight);
} 
parameters {
  real<lower=1,upper=4> logy_min;
  real<lower=0,upper=4> alpha;
}
model {
  exp(logx) ~ wpareto(exp(logy_min),alpha+s,weight);
}

Here I had to add the weighted pareto as a new function, and calculated the weight as transformed data.

The problem comes in when using the scaling, $s$, with uncertainties as well. The obvious change that needs to be made is to weight the uncertainty distribution contribution in the same way, so that we have

functions {
  real wpareto_log(vector y,real y_min,real alpha,vector weight){
    return sum((log(alpha)+alpha*log(y_min)-(alpha+1)*log(y)) .*weight);
  }
  real wnormal_log(vector x, vector mu, real std, vector weight){
     vector[num_elements(x)] y;
     y <- x-mu;
     return sum(-y .*y .*weight ./(2*std *std));
  }
}
data {
  int<lower=0> N;
  vector[N] logx;
  real<lower=0> std;
  real s; //scale
}
transformed data {
  vector[N] weight;
  for (i in 1:N) weight[i] <- pow(exp(logx[i]),s);
  weight <- weight/mean(weight);
}
parameters {
  real<lower=1,upper=4> logy_min;
  real<lower=0,upper=4> alpha;
  vector<lower=logy_min,upper=10>[N] log_xtrue;
}
model {
  exp(log_xtrue) ~ wpareto(exp(logy_min), alpha+s,weight);
  logx ~ wnormal(log_xtrue, std,weight);
}

(so we just added the wnormal_log function, which added the weighting in).

However, this doesn't work. The mean $alpha$ is estimated as 0.09 (ie. it's just pushing against its lower boundary).

The one thing I can see that could be wrong with this model is that when replacing $x_i$ with $x_i$ copies of itself, we should perhaps be using the standard error on the mean of the $x_i$ (as computed by how many copies there are), ie. we should have $\sigma_{true} = \frac{\sigma}{\sqrt{N}}$, where $N$ here is just the weight of the variate. However, implementing this doesn't help.

Can anyone tell me what I'm doing wrong?

EDIT:

I found a bug in the stan code that I had which was causing this not to work. Working code is:

data {
  int<lower=0> N;            //Number of obs.
  vector[N] l10_x_meas;      //log10 of measurements
  vector<lower=0>[N] sd_dex; //uncertainty of measurements in dex
  real scale;
}
transformed data {
  vector[N] weight; //weighting from scale
  for (i in 1:N) weight[i] <- pow(10,l10_x_meas[i]*scale);
  weight <- weight/mean(weight);
}
parameters {
  real<lower=scale,upper=4> rate;       //slope of pareto
  real<lower=10,upper=14> l10_x_min; //lower truncation of log10(x)
  vector<lower=l10_x_min,upper=18>[N] l10_x_true; //estimates of true x
}
transformed parameters {
  vector<lower=0>[N] y;      //will be transformation of x: y log10(x/xmin) 
  y <- l10_x_true-l10_x_min; //transformation
}
model {
  for (i in 1:N){
    increment_log_prob(exponential_log(y[i],(rate-scale)*log10())*weight[i]);
    increment_log_prob(normal_log(l10_x_true[i],l10_x_meas[i],sd_dex[i])*weight[i]);
  }
}

(Note that I replaced alpha+s with alpha-s). Note also that in this toy example I would not actually use the weighting, but I need it for my actual use-case. Using the weighting here reduces the accuracy of the result.

Answer 1

We have had some discussions about weights on the Stan-development Google group. Although I don't think we reached a resolution, I would say while what you are proposing is somewhat accepted in the frequentist literature, it seems difficult to reconcile it with a Bayesian approach. The weights just do not fit into Bayes' theorem very well because if the weight is greater than 1, you are conditioning on additional data that you have not actually observed (and conversely if the weight is less than 1).

I believe what you are implying is that the standard deviation of the measurement error is a decreasing function of the true measurement. If so, it should be possible to model that directly without resorting to weights.

We know that if the true value has a Pareto distribution, then its logarithm has a shifted exponential distribution. In your case, $\log{x} - \log{x_{min}} = \log{x_{true}} - \log{x_{min}} + \log{e}$ where $e$ is log-normal and $\log{x_{true}} - \log{x_{min}}$ is exponential with rate $\alpha$. So the .stan program could be something like

data {
  int<lower=1> N;
  vector[N] log_x;
}
parameters {
  real<lower=0> alpha;
  vector[N] log_e;
  real<upper=min(log_x - log_e)> log_minimum;
}
transformed parameters {
  vector[N] truth;
  truth <- exp(log_x - log_e);
}
model {
  // can ignore Jacobian warning in this case because it's a linear transformation
  log_x - log_e - log_minimum ~ exponential(alpha);
  {
    vector[N] std_dev;
    for (i in 1:N) std_dev[i] <- inv(truth[i]);
    // might want to estimate a scaling parameter on std_dev here
    log_e ~ normal(0, std_dev);
  }
  // prior on alpha
}

That said, I'm not sure if Stan would be able to sample from the posterior distribution very well due to the rapidly changing curvature of each element of log_e.