1
$\begingroup$

I'm trying to redo an exercise in BUGS from this webpage: a linear regression over a data set with some outliers, using a model that accounts for them. This model uses a mixture of signal and noise (the outliers). The likelihood is:

$$p(x_i,y_i,e_i | \theta, g_i,\sigma_B) = \frac{g_i}{\sqrt{2\pi e_i^2}} \exp \Big\{ -\frac{1}{2e_i^2} (y - \hat{y}(x|\theta))^2 \Big\} \\ + \frac{1-g_i}{\sqrt{2\pi \sigma_B^2}} \exp \Big\{ -\frac{1}{2 \sigma_B^2} (y - \hat{y}(x|\theta))^2 \Big\}$$

where ($x_i,y_i$) is the available data, and $e_i$ the respective error, $\hat{y}(x|\theta)$ is the estimate $\theta_0 + \theta_1 x$. The nuisance parameters $g_i$ range from $0$ to $1$ and state the degree to which the i-th data point fits the model (so, $g_i \approx 0$ indicates an outlier) and $\sigma_B$ is the variance of the Gaussian that models the noise (herein, just some big number, say $50$).

I've tried to code this in BUGS:

model {
  for (i in 1:n) {
    tau[i] <- 1/pow(e[i],2)
    part1[i] ~ dnorm(mu[i], tau[i])

    part2[i] ~ dnorm(mu[i], tauOutliers)

    y[i] <- g[i] * part1[i] + (1-g[i]) * part2[i]

    mu[i] <- theta0 + theta1 * x[i]
    g[i] ~ dunif(0,1)
  }

  theta0 ~ dflat()
  theta1 ~ dflat()
  tauOutliers <- 1/sigmaB
}

However I get a multiple definitions of node y[1] error. I think BUGS does not accept deterministic assignments to observable data.

It also does not accept the more direct y[i] ~ g[i] * dnorm(mu[i], tau[i]) + (1-g[i]) * dnorm(mu[i], tauOutliers).

Does anyone know how to solve this problem in BUGS?

Thanks,

$\endgroup$

1 Answer 1

1
$\begingroup$

Ok, after some extra searching, let me try to answer my own question :-)

Since BUGS cannot sample from an arbitrary distribution, we can use the zeros trick to plug the likelihood directly. For this specific problem:

model {
  for (i in 1:n) {

    phi[i] <- -log( (g[i]/sqrt(2*pi*pow(e[i],2))) 
                          * exp(-0.5*pow(y[i]-mu[i],2)/pow(e[i],2)) 
                  + ((1-g[i])/sqrt(2*pi*pow(sigmaB,2))) 
                          * exp(-0.5*pow(y[i]-mu[i],2)/pow(sigmaB,2)) ) + C
    dummy[i] <- 0
    dummy[i] ~ dpois( phi[i] )

    mu[i] <- theta0 + theta1 * x[i]
    g[i] ~ dunif(0,1)

  }

  theta0 ~ dflat()
  theta1 ~ dflat()

  C <- 10000    # for the zeros trick
  pi <- 3.14159
}

Using this model BUGS was able to find a good solution and the 'suitable' outliers (with $g_i \approx 1/3$):

Outlier and linear fit

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.