I’m using probit model to do regression.

My model are as following,

y_{i} = \begin{equation}
\left\{
\begin{array}{**lr**}
1, z_{i} > 0 & \\
0, otherwise &
\end{array}
\right.
\end{equation}

z_{i} = x_{i}*beta + \epsilon_{i} where \epsilon_{i} \sim N(0,1)

In order to make selection on \beta 's element, I place bayesian lasso prior on \beta, which is as following,

\pi(\beta \mid \lambda) = \prod_{j=1}^{p} \frac{\lambda}{2}e^{-\lambda \mid \beta_{j} \mid}

\pi(\lambda^{2}) \sim gamma(a,b).

My code is as following,

```
data{
int<lower=0> N; // number of obs
int<lower=0> p; // number of predictors
int<lower=0,upper=1> y[N];
matrix[N, p] X; // predictor variables
real<lower=0> a;
real<lower=0> b;
}
parameters{
vector[p] beta;
real<lower=0> lambda;
}
transformed parameters{
vector[N] mu;
mu = X*beta;
}
model{
for(n in 1:N){
y[n] ~ bernoulli(Phi(mu[n]));
}
for(j in 1:p){
beta[j] ~ double_exponential( 0 , inv(sqrt(lambda)));
}
lambda ~ gamma(a,b);
}
```

My problem is that I can recover the true beta.

Here is how I generate data and my result,

```
N <- 400
p <- 4
rho = .35
Sig_beta <- matrix(rho, p, p)
Sig_beta <- Sig_beta^(toeplitz(1:p)-1)
diag(Sig_beta) <- 1
X <- array(0,dim=c(N,p))
X <- mvrnorm(N, rep(0,p), Sig_beta)
beta <- numeric(p)
index <- c(1,0,1,1)
beta[which(index==1)] <- rnorm(length(which(index==1)),3,10)
z <- numeric(N)
for(n in 1:N){
z <- rnorm(N,beta*X[n,],1)
}
y <- ifelse(z>0,1,0)
dat <- list("N" = N , "p" = p, "X" = X , "y" = y , "a" = 1 , "b" = 1)
```

Rplot01.pdf (17.1 KB)

and

beta

[1] -1.978069 0.000000 -5.755343 -1.668030

It is clear that all components of beta are concentrates near 0.

It’s very strange strange and I don’t know where is wrong.