diff --git a/src/MarkovProcesses.jl b/src/MarkovProcesses.jl
index 243ccee291878858c51d76457bda3044753d9fa1..4e0964d0c69ae9e077e8780765708ec997ac2356 100644
--- a/src/MarkovProcesses.jl
+++ b/src/MarkovProcesses.jl
@@ -10,11 +10,11 @@ import DelimitedFiles: readdlm, writedlm
 import Distributed: @everywhere, @distributed, @sync, @async, nworkers, nprocs, workers
 import Distributed: nworkers, nprocs, workers, remotecall_fetch
 import DistributedArrays: DArray, dzeros, convert, localpart
-import Distributions: Product, Uniform, Normal, MvNormal, Categorical
+import Distributions: Uniform, Normal, MvNormal, Categorical
 import Distributions: Distribution, Univariate, Continuous, 
                       UnivariateDistribution, DiscreteUnivariateDistribution,
                       MultivariateDistribution, product_distribution
-import Distributions: insupport, isbounded, ncategories, pdf
+import Distributions: insupport, isbounded, ncategories, pdf, product_distribution
 import FunctionWrappers: FunctionWrapper
 import Logging: @info
 using LinearAlgebra
@@ -39,8 +39,8 @@ function __init__()
 end
 
 ## Exports
-export Distribution, Product, Uniform, Normal
-export @everywhere
+export Distribution, Uniform, Normal
+export product_distribution, @everywhere
 
 # Common types and constructors
 export SVector, @SVector
diff --git a/test/unit/draw_pm.jl b/test/unit/draw_pm.jl
index 0a9e6f635a962b40bcad38bf5539db690ecc2298..c8ca7e2a7e265e993d71717f18caf2006dc371e5 100644
--- a/test/unit/draw_pm.jl
+++ b/test/unit/draw_pm.jl
@@ -3,7 +3,7 @@ using MarkovProcesses
 load_model("ER")
 
 k1 = ER[:k1]
-dist_mv_unif = Product(Uniform.([2.5,6.0], [3.5,7.0]))
+dist_mv_unif = product_distribution(Uniform.([2.5,6.0], [3.5,7.0]))
 pm = ParametricModel(ER, [:k3,:k2], dist_mv_unif)
 draw_model!(pm)
 test1 = 2.5 <= ER[:k3] <= 3.5 && 6.0 <= ER[:k2] <= 7.0 && pm.df == 2
diff --git a/test/unit/model_prior.jl b/test/unit/model_prior.jl
index dd246066cdcac0807afbf6241c19c53b6423e01b..aec19461325e5d6d06cf17f5fa224df9088f228e 100644
--- a/test/unit/model_prior.jl
+++ b/test/unit/model_prior.jl
@@ -9,7 +9,7 @@ pm1 = ParametricModel(ER, (:k2, Uniform(2.0, 4.0)))
 draw_model!(pm1)
 test_all = test_all && 2.0 <= ER[:k2] <= 4.0 && pm1.df == 1
 
-pm2 = ParametricModel(ER, [:k3,:k2], Product(Uniform.([2.5,6.0], [3.5,7.0])))
+pm2 = ParametricModel(ER, [:k3,:k2], product_distribution(Uniform.([2.5,6.0], [3.5,7.0])))
 draw_model!(pm2)
 test_all = test_all && 2.5 <= ER[:k3] <= 3.5 && 6.0 <= ER[:k2] <= 7.0 && pm2.df == 2