# Linear Regression with Outlier accounting in Bugs

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,

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$):