diff --git a/src/GaussianEnsembles.jl b/src/GaussianEnsembles.jl
index 0e509c8..d9fe874 100644
--- a/src/GaussianEnsembles.jl
+++ b/src/GaussianEnsembles.jl
@@ -40,8 +40,9 @@ julia> rand(Wigner(2), 3)
   -0.313208+0.330435im   -0.131337-0.0904235im  -0.481758+0.0im
 ```
 """
-struct GaussianHermite{β} <: ContinuousMatrixDistribution end
-GaussianHermite(β) = GaussianHermite{β}()
+struct GaussianHermite{B} <: ContinuousMatrixDistribution 
+  beta::B
+end
 
 """
 Synonym for GaussianHermite{β}
@@ -57,30 +58,27 @@ The Dyson index `β` is restricted to `β = 1,2` or `4`, for real, complex, and
 """
 rand(d::Type{Wigner{β}}, dims...) where {β} = rand(d(), dims...)
 
-function rand(d::Wigner{1}, n::Int)
-    A = randn(n, n)
-    normalization = 1 / √(2n)
-    return Symmetric((A + A') / 2 * normalization)
-end
-
-function rand(d::Wigner{2}, n::Int)
-    A = randn(n, n) + im*randn(n, n)
-    normalization = √(4*n)
-    return Hermitian((A + A') / normalization)
-end
-
-function rand(d::Wigner{4}, n::Int)
-    #Employs 2x2 matrix representation of quaternions
-    X = randn(n, n) + im*randn(n, n)
-    Y = randn(n, n) + im*randn(n, n)
-    A = [X Y; -conj(Y) conj(X)]
-    normalization = √(8*n)
-    return Hermitian((A + A') / normalization)
+function rand(d::Wigner, n::Int)
+    if d.beta == 1
+      A = randn(n, n)
+      normalization = 1 / √(2n)
+      return Symmetric((A + A') / 2 * normalization)
+    elseif d.beta == 2
+      A = randn(n, n) + im*randn(n, n)
+      normalization = √(4*n)
+      return Hermitian((A + A') / normalization)
+    elseif d.beta == 4
+      #Employs 2x2 matrix representation of quaternions
+      X = randn(n, n) + im*randn(n, n)
+      Y = randn(n, n) + im*randn(n, n)
+      A = [X Y; -conj(Y) conj(X)]
+      normalization = √(8*n)
+      return Hermitian((A + A') / normalization)
+    else
+      throw(ArgumentError("Cannot sample random matrix of size $n x $n for β=$(d.beta)"))
+    end
 end
 
-rand(d::Wigner{β}, n::Int) where {β} =
-    throw(ArgumentError("Cannot sample random matrix of size $n x $n for β=$β"))
-
 function rand(d::Wigner{β}, dims::Int...) where {β}
     if length(dims)==2 && dims[1] == dims[2]
 	return rand(d, dims[1])
@@ -102,14 +100,14 @@ The `β == ∞` case is defined in Edelman, Persson and Sutton, 2012.
 """
 function tridrand(d::Wigner{β}, n::Int) where {β}
     χ(df::Real) = rand(Distributions.Chi(df))
-    if β≤0
-        throw(ArgumentError("β = $β cannot be nonpositive"))
-    elseif isinf(β)
+    if d.beta≤0
+        throw(ArgumentError("β = $(d.beta) cannot be nonpositive"))
+    elseif isinf(d.beta)
         return tridrand(Wigner{Inf}, n)
     else
         normalization = 1 / √(2n)
         Hd = rand(Distributions.Normal(0,2), n)./√2
-        He = [χ(β*i)/√2 for i=n-1:-1:1]
+        He = [χ(d.beta*i)/√2 for i=n-1:-1:1]
         return normalization * SymTridiagonal(Hd, He)
     end
 end
@@ -118,7 +116,7 @@ function tridrand(d::Wigner{β}, dims...) where {β}
     if length(dims)==2 && dims[1] == dims[2]
 	return rand(d, dims[1])
     else
-        throw(ArgumentError("Cannot sample random matrix of size $dims for β=$β"))
+        throw(ArgumentError("Cannot sample random matrix of size $dims for β=$(d.beta)"))
     end
 end
 
@@ -147,10 +145,10 @@ function eigvaljpdf(d::Wigner{β}, λ::AbstractVector{Eigenvalue}) where {β,Eig
     #Calculate normalization constant
     c = (2π)^(-n/2)
     for j=1:n
-        c *= gamma(1 + β/2)/gamma(1 + β*j/2)
+        c *= gamma(1 + d.beta/2)/gamma(1 + d.beta*j/2)
     end
     Energy = sum(λ.^2/2) #Calculate argument of exponential
-    VandermondeDeterminant(λ, β) * exp(-Energy)
+    VandermondeDeterminant(λ, d.beta) * exp(-Energy)
 end
 
 #####################
@@ -186,11 +184,11 @@ julia> rand(GaussianLaguerre(4, 8), 2)
 ## References:
 - Edelman and Rao, 2005
 """
-mutable struct GaussianLaguerre <: ContinuousMatrixDistribution
-  beta::Real
-  a::Real
-end
-const Wishart = GaussianLaguerre
+struct GaussianLaguerre{B,A} <: ContinuousMatrixDistribution 
+  beta::B
+  a::A
+end 
+const Wishart{B,A} = GaussianLaguerre{B,A}
 
 #TODO Check - the eigenvalue distribution looks funky
 #TODO The appropriate matrix size should be calculated from a and one matrix dimension
@@ -300,12 +298,12 @@ Represents a Gaussian-Jacobi ensemble with Dyson index `β`, while
 ## References:
 - Edelman and Rao, 2005
 """
-mutable struct GaussianJacobi <: ContinuousMatrixDistribution
-  beta::Real
-  a::Real
-  b::Real
-end
-const MANOVA = GaussianJacobi
+struct GaussianJacobi{B,A} <: ContinuousMatrixDistribution 
+  beta::B
+  a::A
+  b::A
+end 
+const MANOVA{B,A} = GaussianJacobi{B,A}
 
 """
     rand(d::GaussianJacobi, n::Int)
diff --git a/src/Ginibre.jl b/src/Ginibre.jl
index 21f9df2..0608519 100644
--- a/src/Ginibre.jl
+++ b/src/Ginibre.jl
@@ -2,19 +2,18 @@ export rand, Ginibre
 import Base.rand
 
 """
-    Ginibre(β::Int, N::Int) <: ContinuousMatrixDistribution
+    Ginibre(β::Int) <: ContinuousMatrixDistribution
 
 Represents a Ginibre ensemble with Dyson index `β` living in `GL(N, F)`, the set
 of all invertible `N × N` matrices over the field `F`. 
 
 ## Fields
 - `beta`: Dyson index
-- `N`: Matrix dimension over the field `F`.
 
 ## Examples
 
 ```@example
-julia> rand(Ginibre(2, 3))
+julia> rand(Ginibre(2), 3)
 3×3 Matrix{ComplexF64}:
  0.781329+2.00346im   0.0595122+0.488652im  -0.323494-0.35966im
   1.11089+0.935174im  -0.384457+1.71419im    0.114358-0.360676im
@@ -24,21 +23,20 @@ julia> rand(Ginibre(2, 3))
 ## References:
 - Edelman and Rao, 2005
 """
-struct Ginibre <: ContinuousMatrixDistribution
-   beta::Float64
-   N::Integer
-end
+struct Ginibre{B} <: ContinuousMatrixDistribution 
+    beta::B
+end 
 
 """
-    rand(W::Ginibre)
+    rand(W::Ginibre, n::Int)
 
 Samples a matrix from the Ginibre ensemble.
 
 For `β = 1,2,4`, generates matrices randomly sampled from the real, complex, and quaternion
 Ginibre ensemble, respectively.
 """
-function rand(W::Ginibre)
-    beta, n = W.beta, W.N
+function rand(W::Ginibre, n::Int)
+    beta = W.beta
     if beta==1
         randn(n,n)
     elseif beta==2
diff --git a/src/Haar.jl b/src/Haar.jl
index 03a8459..f495183 100644
--- a/src/Haar.jl
+++ b/src/Haar.jl
@@ -88,9 +88,9 @@ julia> rand(Haar(2), 3)
 ## References:
 - Edelman and Rao, 2005
 """
-mutable struct Haar <: ContinuousMatrixDistribution
-    beta::Real
-end
+struct Haar{B} <: ContinuousMatrixDistribution 
+    beta::B
+end 
 
 # In random matrix theory one often encounters expressions of the form
 #
diff --git a/src/HaarMeasure.jl b/src/HaarMeasure.jl
index 417302c..3908907 100644
--- a/src/HaarMeasure.jl
+++ b/src/HaarMeasure.jl
@@ -26,7 +26,7 @@ implemented in most versions of LAPACK.
 """
 function rand(W::Haar, n::Int, doCorrection::Int=1)
     beta = W.beta
-    M=rand(Ginibre(beta,n))
+    M=rand(Ginibre(beta), n)
     q,r=qr(M)
     if doCorrection==0
         q