What seems to work, and what I ended up doing, was to define a new sampling distribution using the "ones trick" described in the [BUGS manual](http://users.aims.ac.za/~mackay/BUGS/Manuals/Tricks.html) and with and example for jags given [here](https://sourceforge.net/p/mcmc-jags/discussion/610037/thread/5d0047ab). 

For the "ones" trick to work I need to define the [likelihood function](http://en.wikipedia.org/wiki/Likelihood_function) for my new distribution, which is the same as the [probability density function] but with the data as the variable argument instead of the parameters. 

Using the following equality, where $X$ and $Y$ are random variables and $x$ is a constant: 

$ prob(min(X, Y)>x) = prob(X>x\ and\ Y>x) = prob(X>x) \cdot prob(Y>x) $

I can construct the [cumulative density function](http://en.wikipedia.org/wiki/Cumulative_distribution_function) for my special distribution as follows:

$ prob(min(norm1, norm2) < x) \Rightarrow 1 - prob(min(norm1, norm2) > x) \Rightarrow  1 - prob(norm1 > x\ and\ norm2 > x) \Rightarrow 1 - prob(norm1>x) \cdot prob(norm2>x) \Rightarrow 1 - (1 - prob(norm1<x)) \cdot (1 - prob(norm2<x))$

As $prob(norm1 < z)$ is the cumulative density function of a normal distribution the final "special" cumulative density function is:

$1 - (1 - pnorm(x,\mu_1, \sigma_1)) \cdot (1-pnorm(x,\mu_2, \sigma_2))$

To get the probability density function I deviate this using the equality:

 $ h(x)=f_1(x)*f_2(x) \Rightarrow h'(x) = f_2'(x)*f_1(x) + f_1'(x)*f_2(x)$ 

and end up with the following R expression:

    -(-dnorm(x, m1, s1) * (1 - pnorm(x, m2, s2)) + -dnorm(x, m2, s2) * (1 - pnorm(x, m1, s1)))

In jags the final model specification became (notice that s1 has become 1/pow(s1, 2) due to jags using precision instead of SD):

    model{
        for (i in 1:n){
        p[i] <- -(-dnorm(x[i], m1, 1/pow(s1, 2)) * (1 - pnorm(x[i], m2, 1/pow(s2, 2))) + -dnorm(x[i], m2, 1/pow(s2, 2)) * (1 - pnorm(x[i], m1, 1/pow(s1, 2))))
        ones[i] ~ dbern(p[i])
        }
        m1_sd <- 1000
        m1 ~ dnorm(0, 1/pow(m1_sd, 2))
        m2_sd <- 1000
        m2 ~ dnorm(400, 1/pow(m2_sd, 2))
        s1_m <- 400 
        s1_s <- 1000
        s1 ~ dgamma(pow(s1_m,2)/pow(s1_s,2), s1_m/pow(s1_s,2))
        s2_m <- 100 
        s2_s <- 1000
        s2 ~ dgamma(pow(s2_m,2)/pow(s2_s,2), s2_m/pow(s2_s,2))
    }

This model seems to (and should) retrieve the original parameters as the following 10000 sample posteriors show:
![posterior plot][1]


  [1]: https://i.sstatic.net/OJNlR.png