diff --git a/NEWS.md b/NEWS.md
index decc020..91683ee 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -32,6 +32,8 @@
   rather than P(X > q) for consistency with all other `cdf()` methods.
 * The `quantile()` method for `dist_multivariate_normal()` now correctly gives
   the boundaries when `p=0` or `p=1` when `type="equicoordinate"`.
+* The `quantile()` method for `dist_multivariate_normal()` now only square roots
+  the diagonal elements when `type="marginal"`.
 
 
 ## Breaking changes
diff --git a/R/dist_multivariate_normal.R b/R/dist_multivariate_normal.R
index f9ea8d1..1593a53 100644
--- a/R/dist_multivariate_normal.R
+++ b/R/dist_multivariate_normal.R
@@ -60,7 +60,7 @@ quantile.dist_mvnorm <- function(x, p, type = c("equicoordinate", "marginal"),
   type <- match.arg(type)
   q <- if (type == "marginal") {
     stats::qnorm(p, mean = rep(x[["mu"]], each = length(p)),
-                 sd = rep(diag(sqrt(x[["sigma"]])), each = length(p)), ...)
+                 sd = rep(sqrt(diag(x[["sigma"]])), each = length(p)), ...)
   } else {
     require_package("mvtnorm")
     vapply(p, function(p, ...) {