Computationally verifying the equivalence of ridge regression estimates and Bayesian regression estimates I'm trying to show that the numerical estimates of ridge regression's parameter estimates are the same as the MAP parameter estimates of a Bayesian regression model with normal prior distributions. So far I have not been able to achieve numerical agreement. 
I'm just toying with the mtcars data.
Ridge estimates:
> library(MASS)
> fit_ridge <- lm.ridge(mpg ~ cyl + hp + disp + wt,
+ data = mtcars,
+ lambda = 1)
> fit_ridge
                      cyl           hp         disp           wt 
39.133455180 -1.070747007 -0.019314307  0.003855685 -3.255577449 

Bayesian model fit:
> library(rstanarm)
> fit_bayes <- stan_glm(mpg ~ cyl + hp + disp + wt,
+                       family = "gaussian",
+                       data = mtcars,
+                       prior = normal(0, 1),
+                       algorithm = "optimizing",
+                       prior_aux = NULL,
+                       prior_intercept = NULL)
+ > summary(fit_bayes)

Model Info:

 function:     stan_glm
 family:       gaussian [identity]
 formula:      mpg ~ cyl + hp + disp + wt
 algorithm:    optimizing
 priors:       see help('prior_summary')
 observations: 32
 predictors:   5

Estimates:
            Median MAD_SD 2.5% 97.5%
(Intercept) 40.3    2.5   35.5 45.2 
cyl         -1.2    0.6   -2.3 -0.1 
hp           0.0    0.0    0.0  0.0 
disp         0.0    0.0    0.0  0.0 
wt          -3.7    0.9   -5.5 -2.0 
sigma        2.3    0.3    1.8  3.0 

Clearly there is some disagreement between the parameter estimates. Is my model specification incorrect? Or my understanding that these distinct estimators ought to be numerically the same?
 A: I've been able to show numerical agreement using stan (instead of rstanarm) and my own matrix math (instead of lm.ridge or glmnet). The stan model is 
$$y \sim N(X\beta, \sigma^2_e)$$
$$ \beta \sim N(0, 1/\lambda) $$
which is implemented with the stan code:
data {
  int<lower=0> n;
  int<lower=0> p;
  matrix[n, p] mat;
  vector[n] y;
  real<lower=0> lambda;
  real<lower=0> sigma_e;
}

parameters {
  vector[p] beta;
}

model {
  // ridge prior on beta
  beta ~ normal(0,  sqrt(1/(lambda)));
  // likelihood model
  y ~ normal(mat * beta, sigma_e);
}

The MAP is given by :
out <- optimizing(stan_fit,
                  data = list(mat = Xmat,
                              n = nrow(Xmat),
                              p = ncol(Xmat),
                              lambda = lambda,
                              sigma_e = 1,
                              y = y))

where Xmat is the fixed effects design matrix and y is a vector of responses.
The ridge estimator is given by 
ridge_beta <- solve(t(Xmat) %*% Xmat + lambda * diag(ncol(Xmat))) %*% t(Xmat) %*% y

I'm assuming that $\sigma_e$ is a known constant. 
In the examples I've tried, the resulting MAPs and ridge point estimates agree. 
Edit: This shrinks the intercept coefficient too.
