diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0cf723e5..b749db7d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -15,6 +15,7 @@ This project adheres to [Semantic Versioning](http://semver.org/).
     to compute the rotation matrix between two 2D **unit** vectors.
   * Add methods `.axis_angle()` to `UnitComplex` and `UnitQuaternion` in order to retrieve both the
     unit rotation axis and the rotation angle simultaneously. 
+  * Add functions to construct a random matrix with a user-defined distribution: `::from_distribution(...)`.
 
 ## [0.14.0]
 ### Modified
diff --git a/src/core/construction.rs b/src/core/construction.rs
index fa3854cc..52391147 100644
--- a/src/core/construction.rs
+++ b/src/core/construction.rs
@@ -1,19 +1,19 @@
 #[cfg(feature = "arbitrary")]
-use quickcheck::{Arbitrary, Gen};
-#[cfg(feature = "arbitrary")]
 use core::storage::Owned;
+#[cfg(feature = "arbitrary")]
+use quickcheck::{Arbitrary, Gen};
 
-use std::iter;
 use num::{Bounded, One, Zero};
-use rand::{self, Rand, Rng};
+use rand::{self, distributions::Sample, Rand, Rng};
+use std::iter;
 use typenum::{self, Cmp, Greater};
 
 use alga::general::{ClosedAdd, ClosedMul};
 
-use core::{DefaultAllocator, Matrix, MatrixMN, MatrixN, Scalar, Unit, Vector, VectorN};
-use core::dimension::{Dim, DimName, Dynamic, U1, U2, U3, U4, U5, U6};
 use core::allocator::Allocator;
+use core::dimension::{Dim, DimName, Dynamic, U1, U2, U3, U4, U5, U6};
 use core::storage::Storage;
+use core::{DefaultAllocator, Matrix, MatrixMN, MatrixN, Scalar, Unit, Vector, VectorN};
 
 /*
  *
@@ -234,6 +234,17 @@ where
     {
         Self::from_fn_generic(nrows, ncols, |_, _| rand::random())
     }
+
+    /// Creates a matrix filled with random values from the given distribution.
+    #[inline]
+    pub fn from_distribution_generic<Distr: Sample<N> + ?Sized, G: Rng>(
+        nrows: R,
+        ncols: C,
+        distribution: &mut Distr,
+        rng: &mut G,
+    ) -> Self {
+        Self::from_fn_generic(nrows, ncols, |_, _| distribution.sample(rng))
+    }
 }
 
 impl<N, D: Dim> MatrixN<N, D>
@@ -356,6 +367,16 @@ macro_rules! impl_constructors(
                 where N: Zero {
                 Self::from_partial_diagonal_generic($($gargs, )* elts)
             }
+
+            /// Creates a matrix filled with random values from the given distribution.
+            #[inline]
+            pub fn from_distribution<Distr: Sample<N> + ?Sized, G: Rng>(
+                $($args: usize,)*
+                distribution: &mut Distr,
+                rng: &mut G,
+            ) -> Self {
+                Self::from_distribution_generic($($gargs, )* distribution, rng)
+            }
         }
 
         impl<N: Scalar + Rand, $($DimIdent: $DimBound, )*> MatrixMN<N $(, $Dims)*>