diff --git a/src/GaussianEnsembles.jl b/src/GaussianEnsembles.jl
index 0e509c8..404f9d1 100644
--- a/src/GaussianEnsembles.jl
+++ b/src/GaussianEnsembles.jl
@@ -24,11 +24,12 @@ export GaussianHermite, GaussianLaguerre, GaussianJacobi,
 #####################
 
 """
-    GaussianHermite(β::Int) <: ContinuousMatrixDistribution
+    GaussianHermite{β} <: ContinuousMatrixDistribution
+    GaussianHermite(β::Real) -> GaussianHermite{β}()
 
 Represents a Gaussian Hermite ensemble with Dyson index `β`.
 
-`Wigner(β)` is a synonym.
+`Wigner{β}` is a synonym.
 
 ## Examples
 
@@ -41,7 +42,7 @@ julia> rand(Wigner(2), 3)
 ```
 """
 struct GaussianHermite{β} <: ContinuousMatrixDistribution end
-GaussianHermite(β) = GaussianHermite{β}()
+GaussianHermite(β::Real) = GaussianHermite{β}()
 
 """
 Synonym for GaussianHermite{β}
@@ -49,7 +50,7 @@ Synonym for GaussianHermite{β}
 const Wigner{β} = GaussianHermite{β}
 
 """
-    rand(d::Wigner, n::Int)
+    rand(d::Wigner{β}, n::Int)
 
 Generates an `n × n` matrix randomly sampled from the Gaussian-Hermite ensemble (also known as the Wigner ensemble).
 
@@ -64,16 +65,24 @@ function rand(d::Wigner{1}, n::Int)
 end
 
 function rand(d::Wigner{2}, n::Int)
-    A = randn(n, n) + im*randn(n, n)
+    A = randn(ComplexF64, 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)]
+    X = randn(ComplexF64, n, n)
+    Y = randn(ComplexF64, n, n)
+    A = Matrix{ComplexF64}(undef, 2n, 2n)
+    @inbounds for j in 1:n, i in 1:n
+        x = X[i, j]
+        y = Y[i, j]
+        A[i, j] = x
+        A[i+n, j] = -conj(y)
+        A[i, j+n] = y
+        A[i+n, j+n] = conj(x)
+    end
     normalization = √(8*n)
     return Hermitian((A + A') / normalization)
 end
@@ -90,7 +99,7 @@ function rand(d::Wigner{β}, dims::Int...) where {β}
 end
 
 """
-    tridand(d::Wigner, n::Int)
+    tridand(d::Wigner{β}, n::Int)
 
 Generates an `n × n` symmetric tridiagonal matrix from the Gaussian-Hermite ensemble (also known as the Wigner ensemble).
 
@@ -158,15 +167,15 @@ end
 #####################
 
 """
-    GaussianLaguerre(β::Real, a::Real)` <: ContinuousMatrixDistribution
+    GaussianLaguerre{β}(a::Real)` <: ContinuousMatrixDistribution
+    GaussianLaguerre(β::Real, a::Real) -> GaussianLaguerre{β}(a)
 
 Represents a Gaussian-Laguerre ensemble with Dyson index `β` and `a` parameter
 used to control the density of eigenvalues near `λ = 0`.
 
-`Wishart(β, a)` is a synonym.
+`Wishart{β}(a)` is a synonym.
 
 ## Fields
-- `beta`: Dyson index
 - `a`: Parameter used for weighting the joint probability density function of the ensemble
 
 ## Examples
@@ -186,16 +195,16 @@ julia> rand(GaussianLaguerre(4, 8), 2)
 ## References:
 - Edelman and Rao, 2005
 """
-mutable struct GaussianLaguerre <: ContinuousMatrixDistribution
-  beta::Real
+struct GaussianLaguerre{β} <: ContinuousMatrixDistribution
   a::Real
 end
-const Wishart = GaussianLaguerre
+GaussianLaguerre(β::Real, a::Real) = GaussianLaguerre{β}(a::Real)
+const Wishart{β} = GaussianLaguerre{β}
 
 #TODO Check - the eigenvalue distribution looks funky
 #TODO The appropriate matrix size should be calculated from a and one matrix dimension
 """
-    rand(d::GaussianLaguerre, n::Int)
+    rand(d::GaussianLaguerre{β}, n::Int)
 
 Generate a random matrix sampled from the Gaussian Laguerre ensemble (also known as the Wishart ensemble)
 with parameters defined in `d`.
@@ -203,15 +212,24 @@ with parameters defined in `d`.
 The Dyson index `β` is restricted to `β = 1,2` (`n × n` matrix) or `4` (`2n × 2n` block matrix representation),
 for real, complex, and quaternionic fields, respectively.
 """
-function rand(d::GaussianLaguerre, n::Int)
-  a, beta = d.a, d.beta
-  a >= beta * n / 2 || throw(ArgumentError("the minimum value of `a` must be `βn/2`."))
-  m = Int(2*a/beta)
-  if beta == 1 # real
+function rand(d::GaussianLaguerre{1}, n::Int)
+    a = d.a
+    a >=  n / 2 || throw(ArgumentError("the minimum value of `a` must be `βn/2`."))
+    m = Int(2a)
     A = randn(m, n)
-  elseif beta == 2 # complex
+    return (A' * A) / n
+end
+function rand(d::GaussianLaguerre{2}, n::Int)
+    a = d.a
+    a >= n || throw(ArgumentError("the minimum value of `a` must be `βn/2`."))
+    m = Int(2a)
     A = randn(ComplexF64, m, n)
-  elseif beta == 4 # quaternion
+    return (A' * A) / n
+end
+function rand(d::GaussianLaguerre{4}, n::Int)
+    a = d.a
+    a >= 2n || throw(ArgumentError("the minimum value of `a` must be `βn/2`."))
+    m = Int(2a)
     # employs 2x2 matrix representation of quaternions
     X = randn(ComplexF64, m, n)
     Y = randn(ComplexF64, m, n)
@@ -224,26 +242,24 @@ function rand(d::GaussianLaguerre, n::Int)
         A[i, j+n] = y
         A[i+m, j+n]   = conj(x)
     end
-  else
-    error("beta = $(beta) is not implemented")
-  end
-  return (A' * A) / n
+    return (A' * A) / n
 end
+rand(d::GaussianLaguerre{β}, n::Int) where {β} = throw(ArgumentError("beta = $(β) is not implemented"))
 
 """
-    bidrand(d::GaussianLaguerre, n::Int)
+    bidrand(d::GaussianLaguerre{β}, n::Int)
 
 Generate an `n × n` bidiagonal matrix sampled from the Gaussian Laguerre ensemble (also known as the Wishart ensemble).
 """
-function bidrand(d::GaussianLaguerre, m::Integer)
-  if d.a <= d.beta*(m-1)/2.0
-      error("Given your choice of m and beta, a must be at least $(d.beta*(m-1)/2.0) (You said a = $(d.a))")
+function bidrand(d::GaussianLaguerre{β}, m::Integer) where {β}
+  if d.a <= β*(m-1)/2.0
+      error("Given your choice of m and beta, a must be at least $(β*(m-1)/2.0) (You said a = $(d.a))")
   end
-  Bidiagonal([chi(2*d.a-i*d.beta) for i=0:m-1], [chi(d.beta*i) for i=m-1:-1:1], true)
+  Bidiagonal([chi(2*d.a-i*β) for i=0:m-1], [chi(β*i) for i=m-1:-1:1], true)
 end
 
 """
-    tridrand(d::GaussianLaguerre, n::Int)
+    tridrand(d::GaussianLaguerre{β}, n::Int)
 
 Generate an `n × n` tridiagonal matrix sampled from the Gaussian Laguerre ensemble (also known as the Wishart ensemble).
 """
@@ -257,25 +273,25 @@ end
 eigvalrand(d::GaussianLaguerre, m::Integer) = eigvals(tridrand(d, m))
 
 #TODO Check m and ns
-function eigvaljpdf(d::GaussianLaguerre, lambda::Vector{Eigenvalue}) where {Eigenvalue<:Number}
+function eigvaljpdf(d::GaussianLaguerre{β}, lambda::Vector{Eigenvalue}) where {β,Eigenvalue<:Number}
   m = length(lambda)
   #Laguerre parameters
-  p = 0.5*d.beta*(m-1) + 1.0
+  p = 0.5*β*(m-1) + 1.0
   #Calculate normalization constant
   c = 2.0^-(m*d.a)
-  z = (d.a - d.beta*(m)*0.5)
+  z = (d.a - β*(m)*0.5)
   for j=1:m
-     z += 0.5*d.beta
+     z += 0.5*β
      if z < 0 && (int(z) - z) < eps()
          #Pole of gamma function, there is no density here no matter what
          return 0.0
      end
-     c *= gamma(1 + beta/2)/(gamma(1 + beta*j/2)*gamma(z))
+     c *= gamma(1 + β/2)/(gamma(1 + β*j/2)*gamma(z))
   end
 
-  Prod = prod(lambda.^(a-p)) #Calculate Laguerre product term
+  Prod = prod(lambda.^(d.a-p)) #Calculate Laguerre product term
   Energy = sum(lambda)/2 #Calculate argument of exponential
-  return c * VandermondeDeterminant(lambda, beta) * Prod * exp(-Energy)
+  return c * VandermondeDeterminant(lambda, β) * Prod * exp(-Energy)
 end
 
 
@@ -285,37 +301,37 @@ end
 ###################
 
 """
-    GaussianJacobi(β::Real, a::Real, a::Real)` <: ContinuousMatrixDistribution
+    GaussianJacobi{β}(a::Real, b::Real) <: ContinuousMatrixDistribution
+    GaussianJacobi(β::Real, a::Real, b::Real) -> GaussianJacobi{β}(a, b)
 
 Represents a Gaussian-Jacobi ensemble with Dyson index `β`, while
 `a`and `b` are parameters used to weight the joint probability density function of the ensemble.
 
-`MANOVA(β, a, b)` is a synonym.
+`MANOVA{β}(a, b)` is a synonym.
 
 ## Fields
-- `beta`: Dyson index
 - `a`: Parameter used for shaping the joint probability density function near `λ = 0`
 - `b`: Parameter used for shaping the joint probability density function near `λ = 1`
 
 ## References:
 - Edelman and Rao, 2005
 """
-mutable struct GaussianJacobi <: ContinuousMatrixDistribution
-  beta::Real
+struct GaussianJacobi{β} <: ContinuousMatrixDistribution
   a::Real
   b::Real
 end
-const MANOVA = GaussianJacobi
+GaussianJacobi(β::Real, a::Real, b::Real) = GaussianJacobi{β}(a, b)
+const MANOVA{β} = GaussianJacobi{β}
 
 """
-    rand(d::GaussianJacobi, n::Int)
+    rand(d::GaussianJacobi{β}, n::Int)
 
 Generate an `n × n` random matrix sampled from the Gaussian-Jacobi ensemble (also known as the MANOVA ensemble)
 with parameters defined in `d`.
 """
-function rand(d::GaussianJacobi, m::Integer)
-  w1 = Wishart(m, int(2.0*d.a/d.beta), d.beta)
-  w2 = Wishart(m, int(2.0*d.b/d.beta), d.beta)
+function rand(d::GaussianJacobi{β}, n::Int) where {β}
+  w1 = rand(Wishart(β, d.a), n)
+  w2 = rand(Wishart(β, d.b), n)
   return (w1 + w2) \ w1
 end
 
@@ -346,12 +362,12 @@ end
 #   and generalized singular value problems", Foundations of Computational Mathematics,
 #   vol. 8 iss. 2 (2008), pp 259-285.
 #TODO check normalization
-function sprand(d::GaussianJacobi, n::Integer, a::Real, b::Real)
+function sprand(d::GaussianJacobi{β}, n::Integer, a::Real, b::Real) where {β}
   CoordI = zeros(8n-4)
   CoordJ = zeros(8n-4)
   Values = zeros(8n-4)
 
-  c, s, cp, sp = SampleCSValues(n, a, b, d.beta)
+  c, s, cp, sp = SampleCSValues(n, a, b, β)
 
   #Diagonals of each block
   for i=1:n
@@ -392,9 +408,9 @@ function sprand(d::GaussianJacobi, n::Integer, a::Real, b::Real)
 end
 
 #Return n eigenvalues distributed according to the Jacobi ensemble
-function eigvalrand(d::GaussianJacobi, n::Integer)
+function eigvalrand(d::GaussianJacobi{β}, n::Integer) where {β}
   #Generate just the upper left quadrant of the matrix
-  c, s, cp, sp = SampleCSValues(n, d.a, d.b, d.beta)
+  c, s, cp, sp = SampleCSValues(n, d.a, d.b, β)
   dv = [i==1 ? c[n] : c[n+1-i] * sp[n+1-i] for i=1:n]
   ev = [-s[n+1-i]*cp[n-i] for i=1:n-1]
 
@@ -404,28 +420,28 @@ function eigvalrand(d::GaussianJacobi, n::Integer)
 end
 
 #TODO Check m and ns
-function eigvaljpdf(d::GaussianJacobi, lambda::Vector{Eigenvalue}) where {Eigenvalue<:Number}
+function eigvaljpdf(d::GaussianJacobi{β}, lambda::Vector{Eigenvalue}) where {β,Eigenvalue<:Number}
   m = length(lambda)
   #Jacobi parameters
   a1, a2 = d.a, d.b
-  p = 1.0 + d.beta*(m-1)/2.0
+  p = 1.0 + β*(m-1)/2.0
   #Calculate normalization constant
   c = 1.0
   for j=1:m
-    z1 = (a1 - beta*(m-j)/2.0)
+    z1 = (a1 - β*(m-j)/2.0)
     if z1 < 0 && (int(z1) - z1) < eps()
       return 0.0 #Pole of gamma function, there is no density here no matter what
     end
-    z2 = (a2 - beta*(m-j)/2.0)
+    z2 = (a2 - β*(m-j)/2.0)
     if z2 < 0 && (int(z2) - z2) < eps()
       return 0.0 #Pole of gamma function, there is no density here no matter what
     end
-    c *= gamma(1 + beta/2)*gamma(a1+a2-beta*(m-j)/2)
-    c /= gamma(1 + beta*j/2)*gamma(z1)*gamma(z2)
+    c *= gamma(1 + β/2)*gamma(a1+a2-β*(m-j)/2)
+    c /= gamma(1 + β*j/2)*gamma(z1)*gamma(z2)
   end
 
   Prod = prod(lambda.^(a1-p))*prod((1-lambda).^(a2-p)) #Calculate Laguerre product term
   Energy = sum(lambda/2) #Calculate argument of exponential
 
-  return c * VandermondeDeterminant(lambda, beta) * Prod * exp(-Energy)
+  return c * VandermondeDeterminant(lambda, β) * Prod * exp(-Energy)
 end
diff --git a/src/Ginibre.jl b/src/Ginibre.jl
index 21f9df2..154ee6b 100644
--- a/src/Ginibre.jl
+++ b/src/Ginibre.jl
@@ -2,19 +2,16 @@ export rand, Ginibre
 import Base.rand
 
 """
-    Ginibre(β::Int, N::Int) <: ContinuousMatrixDistribution
+    Ginibre{β} <: ContinuousMatrixDistribution
+    Ginibre(β::Real) -> Ginibre{β}()
 
 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,35 +21,28 @@ julia> rand(Ginibre(2, 3))
 ## References:
 - Edelman and Rao, 2005
 """
-struct Ginibre <: ContinuousMatrixDistribution
-   beta::Float64
-   N::Integer
-end
+struct Ginibre{β} <: ContinuousMatrixDistribution end
+Ginibre(β::B) where {B} = Ginibre{β}()
 
 """
-    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
-    if beta==1
-        randn(n,n)
-    elseif beta==2
-        randn(n,n)+im*randn(n,n)
-    elseif beta==4
-        Q0=randn(n,n)
-        Q1=randn(n,n)
-        Q2=randn(n,n)
-        Q3=randn(n,n)
-        [Q0+im*Q1 Q2+im*Q3;-Q2+im*Q3 Q0-im*Q1]
-    else 
-        error(string("beta = ", beta, " not implemented"))
-    end
+rand(W::Ginibre{1}, n::Int) = randn(n, n)
+rand(W::Ginibre{2}, n::Int) = randn(ComplexF64, n, n)
+function rand(W::Ginibre{4}, n::Int)
+    Q0=randn(n,n)
+    Q1=randn(n,n)
+    Q2=randn(n,n)
+    Q3=randn(n,n)
+    return [Q0+im*Q1 Q2+im*Q3;-Q2+im*Q3 Q0-im*Q1]
 end
+rand(W::Ginibre{β}, n::Int) where {β} = throw(ArgumentError("Cannot sample random matrix of size $n x $n for β=$β"))
+
 
 function jpdf(Z::AbstractMatrix{z}) where {z<:Complex}
     pi^(size(Z,1)^2)*exp(-trace(Z'*Z))
diff --git a/src/Haar.jl b/src/Haar.jl
index 03a8459..2694a94 100644
--- a/src/Haar.jl
+++ b/src/Haar.jl
@@ -65,16 +65,14 @@ data(P::Ptr{gsl_permutation}) = [convert(Int64, x)+1 for x in
     pointer_to_array(permutation_data(P), (convert(Int64, permutation_size(P)) ,))]
 
 """
-    Haar(β::Int) <: ContinuousMatrixDistribution
+    Haar{β} <: ContinuousMatrixDistribution
+    Haar(β::Int) -> Haar{β}()
 
 Represents a Haar measure with Dyson index `β`, in which values of `β = 1,2` or `4`
 correspond to matrices are distributed with uniform Haar measure over the
 classical orthogonal, unitary and symplectic groups `O(n)`, `U(n)` and
 `Sp(n)~USp(2n)` respectively.
 
-## Fields
-- `beta`: Dyson index
-
 ## Examples
 
 ```@example
@@ -88,9 +86,8 @@ julia> rand(Haar(2), 3)
 ## References:
 - Edelman and Rao, 2005
 """
-mutable struct Haar <: ContinuousMatrixDistribution
-    beta::Real
-end
+struct Haar{β} <: ContinuousMatrixDistribution end
+Haar(β::Real) = Haar{β}()
 
 # In random matrix theory one often encounters expressions of the form
 #
diff --git a/src/HaarMeasure.jl b/src/HaarMeasure.jl
index 417302c..5e11395 100644
--- a/src/HaarMeasure.jl
+++ b/src/HaarMeasure.jl
@@ -24,42 +24,55 @@ implemented in most versions of LAPACK.
 - Edelman and Rao, 2005
 - Mezzadri, 2006, math-ph/0609050
 """
-function rand(W::Haar, n::Int, doCorrection::Int=1)
-    beta = W.beta
-    M=rand(Ginibre(beta,n))
-    q,r=qr(M)
+function rand(W::Haar{1}, n::Int, doCorrection::Int=1)
+    M = rand(Ginibre(1), n)
+    q, r = qr(M)
     if doCorrection==0
-        q
+        return q
     elseif doCorrection==1
-        if beta==1
-            L = sign.(rand(n).-0.5)
-        elseif beta==2
-            L = exp.(im*rand(n)*2pi)
-        elseif beta==4
-            L = exp.(im*rand(2n)*2pi)
-        else
-            error(string("beta = ",beta, " not implemented."))
-        end
-        q*Diagonal(L)
+        L = sign.(rand(n).-0.5)
+        return q * Diagonal(L)
     elseif doCorrection==2
-        if beta==1
-            L=sign.(diag(r))
-        elseif (beta==2 || beta==4)
-            L=diag(r)
-            L=L./abs.(L)
-        else
-            error(string("beta = ",beta, " not implemented."))
-        end
-        q*Diagonal(L)
+        L = sign.(diag(r))
+        return q * Diagonal(L)
+    end
+end
+function rand(W::Haar{2}, n::Int, doCorrection::Int=1)
+    M = rand(Ginibre(2), n)
+    q, r = qr(M)
+    if doCorrection==0
+        return q
+    elseif doCorrection==1
+        L = exp.(im*rand(n)*2pi)
+        return q * Diagonal(L)
+    elseif doCorrection==2
+        L = diag(r)
+        L = L ./ abs.(L)
+        return q * Diagonal(L)
+    end
+end
+function rand(W::Haar{4}, n::Int, doCorrection::Int=1)
+    M = rand(Ginibre(4), n)
+    q, r = qr(M)
+    if doCorrection==0
+        return q
+    elseif doCorrection==1
+        L = exp.(im*rand(2n)*2pi)
+        return q * Diagonal(L)
+    elseif doCorrection==2
+        L = diag(r)
+        L = L ./ abs.(L)
+        return q * Diagonal(L)
     end
 end
+rand(W::Haar{β}, n::Int, doCorrection::Int=1) where {β} = throw(ArgumentError("Cannot sample random matrix of size $n x $n for β=$β"))
 
 #A utility method to check if the piecewise correction is needed
 #This checks the R part of the QR factorization; if correctly done,
 #the diagonals are all chi variables so are non-negative
 function NeedPiecewiseCorrection()
     n=20
-    R=qr(randn(Ginibre(2,n)))[2]
+    R=qr(randn(Ginibre(2), n))[2]
     return any([x<0 for x in diag(R)])
 end
 
@@ -140,12 +153,7 @@ function Stewart(::Type{T}, n) where {T<:Union{Float64,ComplexF64}}
 end
 
 export randfast
-function randfast(W::Haar, n::Int)
-    if W.beta==1
-        Stewart(Float64, n)
-    elseif W.beta==2
-        Stewart(ComplexF64, n)
-    else
-        error("beta = $beta not implemented")
-    end
-end
+randfast(W::Haar{1}, n::Int) = Stewart(Float64, n)
+randfast(W::Haar{2}, n::Int) = Stewart(ComplexF64, n)
+randfast(W::Haar{β}, n::Int) where {β} = throw(ArgumentError("Cannot sample random matrix of size $n x $n for β=$β"))
+
diff --git a/test/GaussianEnsembles.jl b/test/GaussianEnsembles.jl
index 342cd48..7d80823 100644
--- a/test/GaussianEnsembles.jl
+++ b/test/GaussianEnsembles.jl
@@ -4,7 +4,7 @@ using Test
 @testset "GaussianEnsembles" begin
 
 @test Wigner{3} == GaussianHermite{3}
-@test Wishart == GaussianLaguerre
+@test Wishart{3} == GaussianLaguerre{3}
 
 n = 25
 
@@ -47,8 +47,16 @@ for (β, T, N) in [(1, Real, n), (2, Complex, n), (4, Complex, 2n)]
         #vd = RandomMatrices.VandermondeDeterminant(vals, β)
         #@test isa(vd, Real)
 
-        #ed = eigvaljpdf(d, vals)
-        #@test isa(ed, Real)
+        ed = eigvaljpdf(d, rand(3))
+        @test isa(ed, Real)
+    end
+    @testset "MANOVA (β = $(β))" begin
+        a = 2.0 * (rand(1:5) + β * n)
+        b = a * 3.5
+        d = MANOVA(β, a, b)
+        A = rand(d, n)
+        @test eltype(A) <: T
+        @test size(A) == (N, N)
     end
 end
 
diff --git a/test/Ginibre.jl b/test/Ginibre.jl
new file mode 100644
index 0000000..ec8dc51
--- /dev/null
+++ b/test/Ginibre.jl
@@ -0,0 +1,13 @@
+using RandomMatrices
+using Test
+
+n = 25
+
+for (β, T, N) in [(1, Real, n), (2, Complex, n), (4, Complex, 2n)]
+    @testset "Ginibre (β = $(β))" begin
+        d = Ginibre(β)
+        A = rand(d, n)
+        @test eltype(A) <: T
+        @test size(A) == (N, N)
+    end
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 0212a0f..62176e7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,7 +7,9 @@ using Test
 	include("GaussianEnsembles.jl")
 	include("FormalPowerSeries.jl")
 	include("Haar.jl")
+	include("Ginibre.jl")
 	include("StochasticProcess.jl")
+	include("test_throws.jl")
 	include("doctests.jl")
 	
 	@testset "densities" begin
diff --git a/test/test_throws.jl b/test/test_throws.jl
new file mode 100644
index 0000000..2894105
--- /dev/null
+++ b/test/test_throws.jl
@@ -0,0 +1,12 @@
+using RandomMatrices
+using Test
+
+@testset "throws" begin
+
+@test_throws ArgumentError rand(GaussianLaguerre(10, rand(Float64)), rand(2:5))
+@test_throws ErrorException bidrand(GaussianLaguerre(10, rand(Float64)), rand(2:5))
+@test_throws ArgumentError rand(Ginibre(10), rand(2:5))
+@test_throws ArgumentError rand(Haar(10), rand(2:5))
+@test_throws ArgumentError randfast(Haar(10), rand(2:5))
+
+end # testset