Suppose you know $Y \sim N(\mu_1, \sigma_1^2)$ or $Y \sim N(\mu_2, \sigma_2^2)$. You observe $Y=y$, some realization of the random variable $Y$. What is the probability that $Y \sim N(\mu_1, \sigma_1^2)$?

My intuition is to compare $p$-values from each distribution. Let $p_i$ be the $p$-value for $y$ under $N(\mu_i, \sigma_i^2)$. Here I am thinking of the two-sided $p$-value, $p_i = 2\Phi(-|y-\mu_i|/\sigma_i)$ where $\Phi(x)$ is the standard normal distribution function. I would answer my own question as $p_1/(p_1+p_2)$. But I cannot find any reference that would support this (or even treats this problem).

Here's an unintuitive fact - you're not actually supposed to give glmnet a single value of lambda. From the documentation here: 


  Do not supply a single value for lambda (for predictions after CV use predict() instead).
  Supply instead a decreasing sequence of lambda values. glmnet relies on its
  warms starts for speed, and its often faster to ﬁt a whole path than compute a
  single ﬁt.


cv.glmnet will help you choose lambda, as you alluded to in your examples. The authors of the glmnet package suggest cv$lambda.1se instead of cv$lambda.min, but in practice I've had success with the latter.

After running cv.glmnet, you don't have to rerun glmnet! Every lambda in the grid (cv$lambda) has already been run. This technique is called "Warm Start" and you can read more about it here. Paraphrasing from the introduction, the Warm Start technique reduces running time of iterative methods by using the solution of a different optimization problem (e.g., glmnet with a larger lambda) as the starting value for a later optimization problem (e.g., glmnet with a smaller lambda).

To extract the desired run from cv.glmnet.fit, try this:

small.lambda.index <- which(cv$lambda == cv$lambda.min)
small.lambda.betas <- cv$glmnet.fit$beta[, small.lambda.index]


Revision (1/28/2017)

No need to hack to the glmnet object like I did above; take @alex23lemm's advice below and pass the s = "lambda.min", s = "lambda.1se" or some other number (e.g., s = .007) to both coef and predict. Note that your coefficients and predictions depend on this value which is set by cross validation. Use a seed for reproducibility! And don't forget that if you don't supply an "s" in coef and predict, you'll be using the default of s = "lambda.1se". I have warmed up to that default after seeing it work better in a small data situation. s = "lambda.1se" also tends to provide more regularization, so if you're working with alpha > 0, it will also tend towards a more parsimonious model. You can also choose a numerical value of s with the help of plot.glmnet to get to somewhere in between (just don't forget to exponentiate the values from the x axis!).
