I am trying to define the Normal-Laplace density as stated here at the bottom of Page 5 https://arxiv.org/pdf/1712.07216.pdf but also supporting a location in the Normal component. I will attach my code here, but find the examples in the documentation hard to understand in this context so there may be some obvious issues with it:

```
functions {
real mills_ratio(real z) {
// return (1 - normal_cdf(z, 0, 1)) / exp(normal_lpdf(z | 0, 1));
return (1 - Phi_approx(z)) / exp(std_normal_lpdf(z));
}
real normal_laplace_lpdf(vector y, real mu, real sigma, real scale) {
real total = 0.0;
for (i in 1:num_elements(y)) {
total += log((1 / sqrt(2) * scale) * std_normal_lpdf((y[i] - mu) / sigma) * (mills_ratio((sqrt(2) * sigma) / scale - (y[i] - mu) / sigma) + mills_ratio((sqrt(2) * sigma) / scale + (y[i] - mu) / sigma)));
}
return total;
}
}
```

I am looking to call it like this:

`y2 ~ normal_laplace(mu, sqrt(sigma2), scale);`

Currently this results in the following working with pystan:

Rejecting initial value:

Log probability evaluates to log(0), i.e. negative infinity.

Stan canâ€™t start sampling from this initial value.

Rejecting initial value:

Log probability evaluates to log(0), i.e. negative infinity.

Stan canâ€™t start sampling from this initial value.

Rejecting initial value:

Log probability evaluates to log(0), i.e. negative infinity.

Stan canâ€™t start sampling from this initial value.

Rejecting initial value:

Log probability evaluates to log(0), i.e. negative infinity.

Stan canâ€™t start sampling from this initial value.

Rejecting initial value:

Log probability evaluates to log(0), i.e. negative infinity.

Stan canâ€™t start sampling from this initial value.

Rejecting initial value:

Log probability evaluates to log(0), i.e. negative infinity.

Stan canâ€™t start sampling from this initial value.

Rejecting initial value:

Log probability evaluates to log(0), i.e. negative infinity.

Stan canâ€™t start sampling from this initial value.Initialization between (-2, 2) failed after 100 attempts.

Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

y2 is an input of samples taken from a Normal(5,2) distribution with some Laplace(0,1.5) noise applied to them, scale is given as a known quantity to the stan model in the data section.

Any help would be much appreciated.