2
$\begingroup$

When fitting models using abc, data is simulated using parameters drawn from the prior. The distance between the simulated data and the observed data is calculated, and typically if less than a certain threshold, the parameters are accepted as being part of the posterior.

Ive seen instead of using an arbitrary distance measure, that the distance between the observed and simulated data has been assumed to have a distribution, the likelihood of this calculated, and this likelihood used to derive the probability of acceptance for the proposed parameters (https://arxiv.org/abs/0811.3355).

Can non-parametric likelihood measures be used instead such as the Kolmogorov–Smirnov test?

In the below example I fit the same model using stan and abc rejection to data which is a mixture of a normal distribution and a binary proportion. As the distance for abc I used the likelihood derived from the Kolmogorov–Smirnov p value and proportion test. Both methods seem to give similar results.

rm(list = ls(all = TRUE))
library(rstan)
Sys.setenv(LOCAL_CPPFLAGS = '-march=corei7 -mtune=corei7')
options(mc.cores = parallel::detectCores())

##rejecton abc
#draw from prior
priordraw = function() {
  theta =vector()
  theta[1] = runif(1,0,1)
  theta[2] = runif(1,-5,5)
  theta[3] = runif(1,0,10)
  theta
}

#calculate pval as for dist between observed and simulated data
disfn = function(observed,simulated) {
  prop_p= log(prop.test(c(N-observed[[1]],N-simulated[[1]]),c(N,N))$p.value)
  if (length(simulated[[2]])>1){
    ks_p = log(ks.test(observed[[2]],simulated[[2]])$p.value)
  } else {
    ks_p = 0
  }
  prop_p+ks_p
}

#simulate data using using params drawn from prior and calculate distance from observed data
simdata = function(theta){
  pos = rbinom(1,N,theta[1])
  mes = rnorm(pos,theta[2],theta[3])
  simulated = list(pos=pos,mes=mes)
  dist = disfn(observed,simulated)
}

#observed data
set.seed(1234)
N = 1e3
observed = list(pos=400,mes = rnorm(600,1,2))

#do abc rejection
reps = 1e5
pval = vector()
priors = matrix(NA,reps,3)
for (x in 1:reps){
  priors[x,] = priordraw()
  pval[x] = simdata(priors[x,])
}
#sample prop to lik
inds = 1:length(pval)
selected = sample(inds,size=500,replace=F,prob=exp(pval))
par(mar=c(2,4,2,1),mfrow=c(3,1)) 
hist(priors[selected,1])
hist(priors[selected,2])
hist(priors[selected,3])

quantile(priors[selected,1],c(0.05,0.5,0.95))
quantile(priors[selected,2],c(0.05,0.5,0.95))
quantile(priors[selected,3],c(0.05,0.5,0.95))

###fit with stan
stancode <- 'data {
  int N;
  int pos;
  real mes[N-pos];
}
parameters {
  real<lower=0,upper=1> prop;
  real mu;
  real<lower=0> sigma;
}
model {
  prop ~ uniform(0,1);
  mu ~ uniform(-5,5);
  sigma ~ uniform(0,10);
  target += binomial_lpmf(pos | N, prop);

  for (i in 1:pos) {
      target += normal_lpdf(mes[i]| mu, sigma);
  }
}'

#mixture of proportion and normal
dat = list(pos=observed[[1]],mes = observed[[2]], N=1e3)
fit <- stan(model_code = stancode, data = dat, iter=2000,chains=4)
print(fit)
```
Improve this question
$\endgroup$
4
  • $\begingroup$ We wrote a paper using Wasserstein distances to compare the empirical distributions of the actual sample and the observed sample. $\endgroup$
    – Xi'an
    Commented May 23, 2020 at 6:30
  • $\begingroup$ Thanks for the reference. Initially I tried using the Kolmogorov–Smirnov or Wasserstein distance alone. Since the data is a mixture of a binomial proportion and a normal distribution (e.g. rainfall in the last day) this did not work so well, I think because the empirical cdf is dominated by the proportion, and, if using the Kolmogorov–Smirnov distance as the metric a e.g 5% difference in the proportion is effectively given the same weight as a 5% difference in the ecdf of the normally distributed bit. $\endgroup$
    – hugh
    Commented May 23, 2020 at 22:11
  • $\begingroup$ Hence I wanted to use a different metric for each part of the mixture. A logical way to assign a weight to differences in each part of the mixture seemed to be multiplying the Kolmogorov–Smirnov p value with the proportion test p-val. It seems to work but not sure it is 'correct'. $\endgroup$
    – hugh
    Commented May 23, 2020 at 22:12
  • $\begingroup$ I will update the example to include comparison with using the Wasserstein distance $\endgroup$
    – hugh
    Commented May 24, 2020 at 7:06

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.