Bernoulli decision model consistently underestimates parameter

Modeling
1

I am modelling choice data from a task in which participants made two consecutive right/left decisions about stimuli (s1 and s2) and, if the first decision was correct, the second stimulus would go to the right. They would then use this information to help them in the second decision (choice), by shifting their decision criterion (crit). x1 and x2 are internal signals (in the brain) from s1 and s2 that we cannot measure. The parameter w captures to what extent they shift their criterion (optimal observer would have w=1), and is my parameter of interest.

Here is my current model, which runs well and samples nicely, but always underestimates w compared to simulated data, simulated from code below. I guess the phi going into the likelihood is too low for some reason.
When I run a model using the analytically worked-out likelihood of a rightward decision given the stimulus and previous choice (without sampling x1 and x2), it works well but is very inefficient and seems like a strange approach for STAN.

Any ideas for why this is happening? I really appreciate any pointers or resources, thanks so much!

data {
  int<lower=0> N; 
  vector[N] s1; 
  vector[N] s2; 
  int<lower=0,upper=1> choice[N]; 
}

parameters {
  real<lower=0> w; 
  vector[N] x1; //internal sample from s1 (latent variable)
  vector[N] x2; //internal sample from s2 (latent variable)
}

transformed parameters {
  vector[N] crit; 
  vector<lower=0,upper=1>[N] phi;

  crit = -fabs(x1)/w; 
  phi = Phi_approx(x2-crit);
}

model {
  w ~ lognormal(1,1);

  x1 ~ normal(s1, 1);
  x2 ~ normal(s2, 1);
  choice ~ bernoulli(phi); 
}

Simulation Code:
simulateChoiceData ← function(w, s1, s2, n)
{
data ← data.frame(s1 = numeric(n),
s2 = numeric(n),
x1 = numeric(n),
x2 = numeric(n))

data$direction1 = sample(c(-1,1), n, replace=TRUE)
data$s1 = data$direction1s1
data$x1 = rnorm(n,data$s1,1)
data$choseRight1 = data$x1>0
data$correct1 = (data$x1data$direction1)>0

data$direction2 = ifelse(data$correct1==1,1,-1)
data$s2 = data$direction2*s2
data$x2 = rnorm(n, data$s2, 1)

data$crit = -abs(data$x1)/w
data$choice = data$x2>data$crit
data$correct2 = data$correct1==data$choice

return(data)
}
data = simulateChoiceData(w=2,s1=0.5,s2=1,n=720)