Bayesian Mixture Model Gibbs Sampler for two linear relationships I am attempting to use a Gibbs Sampler to model a mixture of two groups, where the group membership is defined by a linear relationship conditional on x. Both groups have the same slope and intercept, but one group has an additional, positive constant added to the common intercept. My question is that my set up does not seem to be correctly separating out the two groups. (see plot with result below)
$$
z_i = \begin{cases}
  1, & \text{if observation $i$'s response is in group 1 (with positive shift)}, \\
  0, & \text{if observation $i$'s response is in group 2}.
\end{cases}
$$
$$
y_i = \begin{cases}
  \text{N} (\gamma + \beta x_i + \tau z_i\, , \, \sigma^2 )  & \text{ with prob. } \lambda, \\
  \text{N} (\gamma + \beta x_i \,,\, \sigma^2)  & \text{ with prob. } 1- \lambda.
\end{cases}
$$
The (simulated) data look like this:

And here is my code:
#synthetic data a synthetic dataset of 100 points
truebeta <-0.003
truegamma <-0
truetau <-15
truesigma2 <-1
truez <- rbinom(100, 1, 0.50)
truex <- seq(from=1, to= 1000, by = 10)
truey <- rnorm(100, truegamma+truebeta*truex+truetau*truez, sqrt(truesigma2+truesigma2*truez))
plot(truex, truey)

xs <-truex
ys <-truey
n <-length(xs)

niters <-1000

logcondfxn <- function(y, z, gamma, beta, x, tau, sigma2){
   ndist1 <-dnorm(y[z==1], gamma + beta*x[z==1] + tau, sqrt(sigma2))
   ndist2 <-dnorm(y[z==0],gamma + beta*x[z==0], sqrt(sigma2))
   sum(log(ndist1)[is.finite(log(ndist1))]) + sum(log(ndist2)[is.finite(log(ndist2))])
}

library(MCMCpack)

#initialize values
lambda <- rep(0.5, niters)
beta <- rep(0.001, niters)
tau <- rep(0, niters)
gamma <- rep(0, niters)
sigma2 <- rep(5, niters)

### Gibbs Sampler with M-H ####
for (i in 2: niters) {
   #sample zs from bernoulli with prob 
   znum <-lambda[i-1]*(dnorm(ys, mean=gamma[i-1]+beta[i-1]*xs+tau[i-1], sd=sqrt(sigma2[i-1])))
   zdenom1 <- lambda[i-1]*(dnorm(ys, mean=gamma[i-1]+beta[i-1]*xs+tau[i-1], sd=sqrt(sigma2[i-1])))
   zdenom2 <- (1-lambda[i-1])*dnorm(ys, mean=gamma[i-1]+beta[i-1]*xs, sd=sqrt(sigma2[i-1]))
   prob <-znum/(zdenom1+zdenom2)
   prob[is.na(prob)] <-0
   z <- rbinom(n,1, prob)

   #sample gamma
   ystar <- ys - beta[i-1]*xs- tau[i-1]*z
   #what is left from this is the intercept gamma regress ystar on x
   gamma[i] <- rnorm(1, summary(lm(ystar~1))$coef[1], sqrt(sigma2[i-1]/n ))

   #sample beta
   ystar <- ys - gamma[i] - tau[i-1]*z
   beta[i] <- rnorm(1, summary(lm(ystar~xs))$coef[2], sqrt(sigma2[i-1]/n ))

   #sample tau
   ystar <- ys[z==1] - beta[i]*xs[z==1] - gamma[i]
   tau[i] <- rnorm(1, min(mean(ystar), 0, na.rm = TRUE), sqrt(sigma2[i-1]/n))

   #sample lambda
   lambda[i] <-rbeta(1, sum(z)+1, n - sum(z)+1)

   #sample sigma2 using metropolis-hastings
   sigma2.prop <- rgamma(1, sigma2[i-1], 1)
   Mnum <- logcondfxn(ys, z, gamma[i], beta[i], xs, tau[i], sigma2.prop)
   Mnum[is.na(Mnum)] <-0
   Mden <- logcondfxn(ys, z, gamma[i], beta[i], xs, tau[i], sigma2[i-1])
   Mden[is.na(Mden)] <-0

   #Hastings step to correct of assymmetry of proposed sigma2 dist
   logHr <-log(dgamma(sigma2[i-1], sigma2.prop, 1)) - log(dgamma(sigma2.prop, sigma2[i-1], 1))

   #Calc logr
   logr <- Mnum - Mden + logHr

   logu <- log(runif(1, 0, 1))
   sigma2[i] <- ifelse(logr>=logu, sigma2.prop, sigma2[i-1])
}

I tried plotting the iterations to gauge convergence, and the chain is not converging.
plot(sigma2, type = "l")
plot(lambda, type = "l", xlim = c(0,100))
plot(beta, type = "l")
plot(gamma, type = "l", xlim= c(0,50))
plot(tau, type = "l", xlim = c(0,50))

And to diagnose autocorrelation:
acf(beta)
acf(lambda)
acf(gamma)
acf(tau)

I haven't thinned or thrown out burn-in yet, but from these preliminary plots, I am not getting what I hope to (regression lines through each group of points). Instead, group membership goes virtually all to one group (0) in the z vector, and both lines are between the two groups.
plot(xs,ys)
abline(mean(gamma), mean(beta))
abline(mean(gamma) + mean(tau), mean(beta))

Here is a plot of what I am getting right now:

 A: Here is the corrected code to run this problem:
lambda<- rep(0.5, niters)
beta<- rep(0.001, niters)
tau<- rep(0, niters)
gamma<- rep(0, niters)
sigma2<- rep(5, niters)


for (i in 2: niters){
  ## 1. Sample z. sample zs from bernoulli with prob 
  znum<-lambda[i-1]*(dnorm(ys, mean=gamma[i-1]+beta[i-1]*xs+tau[i-1], sd=sqrt(sigma2[i-1])))
  zdenom1<- lambda[i-1]*(dnorm(ys, mean=gamma[i-1]+beta[i-1]*xs+tau[i-1], sd=sqrt(sigma2[i-1])))
  zdenom2<- (1-lambda[i-1])*dnorm(ys, mean=gamma[i-1]+beta[i-1]*xs, sd=sqrt(sigma2[i-1]))

  prob<-znum/(zdenom1+zdenom2)

  #prob[is.na(prob)]<-0

  z<- rbinom(n,1, prob)

  ## 2. Sample lamda from Beta dist
  lambda[i]<-rbeta(1, sum(z)+1, n - sum(z)+1)

  ## 3. Sample mus (beta and gamma are common, tau is only for group1)
  #### 3a. GAMMA
      ystar<- ys - beta[i-1]*xs- tau[i-1]*z
   #what is left from this is the intercept gamma regress ystar on x
   gamma[i]<- rnorm(1, summary(lm(ystar~1))$coef[1], sqrt(sigma2[i-1]/n ))
      #### 3b. BETA
      ystar<- ys - gamma[i] - tau[i-1]*z
      beta[i]<- rnorm(1, summary(lm(ystar~xs-1))$coef[1],vcov(lm(ystar~xs-1)))
  #### 3c. TAU
  n1<- max(sum(z), 0.00001)
  ystar<- ys[z==1] - beta[i]*xs[z==1] - gamma[i]
  tau[i]<- rnorm(1, max(mean(ystar), 0, na.rm = TRUE), sqrt(sigma2[i-1]/n1))

  ## 4. Sample sigma2 from invgamma
  sigma2[i]<- 1/rgamma(1, n/2 +1 , sum((ys - gamma[i]- beta[i]*xs - tau[i]*z)^2, na.rm = TRUE)/ 2)
}

