diff --git a/src/bibtex/all.bib b/src/bibtex/all.bib
index 83c3d2479..6d18607fd 100644
--- a/src/bibtex/all.bib
+++ b/src/bibtex/all.bib
@@ -358,6 +358,24 @@ @article{Marsaglia:1972
   year={1972}
 }
 
+@article{Muller:1959,
+  author = {Muller, Mervin E.},
+  title = {A Note on a Method for Generating Points Uniformly on n-Dimensional Spheres},
+  year = {1959},
+  issue_date = {April 1959},
+  publisher = {Association for Computing Machinery},
+  address = {New York, NY, USA},
+  volume = {2},
+  number = {4},
+  issn = {0001-0782},
+  url = {https://doi.org/10.1145/377939.377946},
+  doi = {10.1145/377939.377946},
+  journal = {Commun. ACM},
+  month = {apr},
+  pages = {19–20},
+  numpages = {2}
+}
+
 @book{HopcroftMotwani:2006,
   title={Introduction to Automata Theory, Languages, and Computation},
   author={Hopcroft, John E. and Rajeev Motwani},
diff --git a/src/reference-manual/transforms.Rmd b/src/reference-manual/transforms.Rmd
index e7226bcdc..4e19bbd55 100644
--- a/src/reference-manual/transforms.Rmd
+++ b/src/reference-manual/transforms.Rmd
@@ -663,7 +663,7 @@ $$
 To generate a unit vector, Stan generates points at random in
 $\mathbb{R}^n$ with independent unit normal distributions, which are then
 standardized by dividing by their Euclidean length.
-@Marsaglia:1972 showed this generates points uniformly at random
+@Muller:1959 showed this generates points uniformly at random
 on $S^{n-1}$.  That is, if we draw $y_n \sim \mathsf{Normal}(0, 1)$
 for $n \in 1{:}n$, then $x = \frac{y}{\Vert y \Vert}$ has a uniform
 distribution over $S^{n-1}$.  This allows us to use an $n$-dimensional
diff --git a/src/stan-users-guide/hyperspherical-models.Rmd b/src/stan-users-guide/hyperspherical-models.Rmd
index d089a32a9..06ed9661a 100644
--- a/src/stan-users-guide/hyperspherical-models.Rmd
+++ b/src/stan-users-guide/hyperspherical-models.Rmd
@@ -93,7 +93,7 @@ compact in that the distance between any two points is bounded.
 
 Stan (inverse) transforms arbitrary points in $\mathbb{R}^{K+1}$ to points
 in $S^K$ using the auxiliary variable approach of
-@Marsaglia:1972.  A point $y \in \mathbb{R}^K$ is transformed to a
+@Muller:1959.  A point $y \in \mathbb{R}^K$ is transformed to a
 point $x \in S^{K-1}$ by
 $$
 x = \frac{y}{\sqrt{y^{\top} y}}.
@@ -101,7 +101,7 @@ $$
 
 The problem with this mapping is that it's many to one; any point
 lying on a vector out of the origin is projected to the same point on
-the surface of the sphere.  @Marsaglia:1972 introduced an
+the surface of the sphere.  @Muller:1959 introduced an
 auxiliary variable interpretation of this mapping that provides the
 desired properties of uniformity; the reference manual contains the
 precise definitions used in the chapter on constrained parameter