diff --git a/src/CalibrateEmulateSample.jl b/src/CalibrateEmulateSample.jl
index 247b6b9e8..9a01f3945 100644
--- a/src/CalibrateEmulateSample.jl
+++ b/src/CalibrateEmulateSample.jl
@@ -4,7 +4,6 @@ using Distributions, Statistics, LinearAlgebra, DocStringExtensions
 
 # No internal deps, light external deps
 include("Observations.jl")
-include("Priors.jl")
 include("ParameterDistribution.jl")
 
 # No internal deps, heavy external deps
diff --git a/src/EKP.jl b/src/EKP.jl
index 72042bd5b..27470faa0 100644
--- a/src/EKP.jl
+++ b/src/EKP.jl
@@ -1,6 +1,8 @@
 module EKP
 
-using ..Priors
+#using ..Priors
+using ..ParameterDistributionStorage
+
 using Random
 using Statistics
 using Distributions
@@ -46,8 +48,6 @@ $(DocStringExtensions.FIELDS)
 struct EKObj{FT<:AbstractFloat, IT<:Int, P<:Process}
     "vector of arrays of size N_ensemble x N_parameters containing the parameters (in each EK iteration a new array of parameters is added)"
      u::Vector{Array{FT, 2}}
-     "vector of parameter names"
-     unames::Vector{String}
      "vector of observations (length: N_data); mean of all observation samples"
      obs_mean::Vector{FT}
      "covariance of the observational noise, which is assumed to be normally distributed"
@@ -66,7 +66,6 @@ end
 
 # outer constructors
 function EKObj(parameters::Array{FT, 2},
-               parameter_names::Vector{String},
                obs_mean,
                obs_noise_cov::Array{FT, 2},
                process::P;
@@ -85,31 +84,21 @@ function EKObj(parameters::Array{FT, 2},
     # timestep store
     Δt = Array([Δt])
 
-    EKObj{FT, IT, P}(u, parameter_names, obs_mean, obs_noise_cov, N_ens, g,
+    EKObj{FT, IT, P}(u, obs_mean, obs_noise_cov, N_ens, g,
                      err, Δt, process)
 end
 
 
 """
-construct_initial_ensemble(N_ens::IT, priors::Array{Prior, 1}; rng_seed=42) where {IT<:Int}
+construct_initial_ensemble(prior::ParameterDistribution, N_ens::IT; rng_seed=42) where {IT<:Int}
 
 Construct the initial parameters, by sampling N_ens samples from specified
-prior distributions.
+prior distribution. Returned with parameters as rows
 """
-function construct_initial_ensemble(N_ens::IT, priors::Array{Prior, 1}; rng_seed=42) where {IT<:Int}
-    N_params = length(priors)
-    params = zeros(N_ens, N_params)
+function construct_initial_ensemble(prior::ParameterDistribution, N_ens::IT; rng_seed=42) where {IT<:Int}
     # Ensuring reproducibility of the sampled parameter values
     Random.seed!(rng_seed)
-    for i in 1:N_params
-        prior_i = priors[i].dist
-        if !(typeof(priors[i].dist) == Deterministic{Float64})
-            params[:, i] = rand(prior_i, N_ens)
-        else
-            params[:, i] = prior_i.value * ones(N_ens)
-        end
-    end
-
+    params = permutedims(sample_distribution(prior, N_ens), (2,1)) #this transpose is [N_ens x dim(param space)]
     return params
 end
 
diff --git a/src/MCMC.jl b/src/MCMC.jl
index 98f3b87b0..ab4e04fef 100644
--- a/src/MCMC.jl
+++ b/src/MCMC.jl
@@ -1,7 +1,8 @@
 module MCMC
 
 using ..GPEmulator
-using ..Priors
+#using ..Priors
+using ..ParameterDistributionStorage
 
 using Statistics
 using Distributions
@@ -34,7 +35,7 @@ struct MCMCObj{FT<:AbstractFloat, IT<:Int}
     "covariance of the observational noise"
     obs_noise_cov::Array{FT, 2}
     "array of length N_parameters with the parameters' prior distributions"
-    prior::Array{Prior, 1}
+    prior::ParameterDistribution
     "MCMC step size"
     step::Array{FT}
     "Number of MCMC steps that are considered burnin"
@@ -69,7 +70,7 @@ where max_iter is the number of MCMC steps to perform (e.g., 100_000)
 function MCMCObj(
     obs_sample::Vector{FT},
     obs_noise_cov::Array{FT, 2},
-    priors::Array{Prior, 1},
+    prior::ParameterDistribution,
     step::FT,
     param_init::Vector{FT},
     max_iter::IT,
@@ -100,7 +101,7 @@ function MCMCObj(
     end
     MCMCObj{FT,IT}(obs_sample,
                    obs_noise_cov,
-                   priors,
+                   prior,
                    [step],
                    burnin,
                    param,
@@ -124,7 +125,14 @@ end
 
 
 function get_posterior(mcmc::MCMCObj)
-    return mcmc.posterior[mcmc.burnin+1:end, :]
+    #Return a parameter distributions object
+    posterior_samples = Samples(mcmc.posterior[mcmc.burnin+1:end,:]; params_are_columns=false) 
+    parameter_constraints = get_all_constraints(mcmc.prior) #live in same space as prior
+    parameter_names = get_name(mcmc.prior) #the same parameters as in prior
+    posterior_distribution = ParameterDistribution(posterior_samples, parameter_constraints, parameter_names)
+    return posterior_distribution
+    #return mcmc.posterior[mcmc.burnin+1:end, :]
+    
 end
 
 function mcmc_sample!(mcmc::MCMCObj{FT}, g::Vector{FT}, gvar::Vector{FT}) where {FT}
@@ -174,30 +182,17 @@ end
 
 
 function log_prior(mcmc::MCMCObj{FT}) where {FT}
-    log_rho = FT[0]
-    # Assume independent priors for each parameter
-    priors = [mcmc.prior[i].dist for i in 1:length(mcmc.prior)]
-    for (param, prior_dist) in zip(mcmc.param, priors)
-        # get density at current parameter value
-        log_rho[1] += logpdf(prior_dist, param)
-    end
-
-    return log_rho[1]
+    return get_logpdf(mcmc.prior,mcmc.param)
 end
 
 
 function proposal(mcmc::MCMCObj)
 
-    variances = zeros(length(mcmc.param))
-    priors = [mcmc.prior[i].dist for i in 1:length(mcmc.prior)]
-    param_names = [mcmc.prior[i].param_name for i in 1:length(mcmc.prior)]
-    for (idx, prior) in enumerate(priors)
-        variances[idx] = var(prior)
-    end
-
+    proposal_covariance = get_cov(mcmc.prior)
+ 
     if mcmc.algtype == "rwm"
         prop_dist = MvNormal(zeros(length(mcmc.param)), 
-                             (mcmc.step[1]^2) * Diagonal(variances))
+                             (mcmc.step[1]^2) * proposal_covariance)
     end
     sample = mcmc.posterior[1 + mcmc.iter[1], :] .+ rand(prop_dist)
 
diff --git a/src/ParameterDistribution.jl b/src/ParameterDistribution.jl
index a9466d2cd..e4a567b25 100644
--- a/src/ParameterDistribution.jl
+++ b/src/ParameterDistribution.jl
@@ -2,6 +2,7 @@ module ParameterDistributionStorage
 
 ## Imports
 using Distributions
+using Statistics
 using StatsBase
 using Random
 
@@ -13,17 +14,18 @@ export ParameterDistribution
 export Constraint
 
 #functions
-export get_name, get_distribution
+export get_name, get_distribution, get_total_dimension, get_dimensions, get_all_constraints, get_n_samples
 export sample_distribution
 export no_constraint, bounded_below, bounded_above, bounded
 export transform_constrained_to_unconstrained, transform_unconstrained_to_constrained
-
+export get_logpdf, get_cov, get_var, get_mean
+    
 ## Objects
 # for the Distribution
 abstract type ParameterDistributionType end
 
 """
-    Parameterized <: ParameterDistributionType
+    struct Parameterized <: ParameterDistributionType
     
 A distribution constructed from a parametrized formula (e.g Julia Distributions.jl)
 """
@@ -32,13 +34,13 @@ struct Parameterized <: ParameterDistributionType
 end
 
 """
-    Samples{FT<:Real} <: ParameterDistributionType
+    struct Samples{FT<:Real} <: ParameterDistributionType
 
 A distribution comprised of only samples, stored as columns of parameters
 """
 struct Samples{FT<:Real} <: ParameterDistributionType
     distribution_samples::Array{FT,2} #parameters are columns
-    Samples(distribution_samples::Array{FT,2}; params_are_columns=true) where {FT <: Real} = params_are_columns ? new{FT}(distribution_samples) : new{FT}(permutedims(distribution_samples,[2,1]))
+    Samples(distribution_samples::Array{FT,2}; params_are_columns=true) where {FT <: Real} = params_are_columns ? new{FT}(distribution_samples) : new{FT}(permutedims(distribution_samples,(2,1)))
     #Distinguish 1 sample of an ND parameter or N samples of 1D parameter, and store as 2D array  
     Samples(distribution_samples::Array{FT,1}; params_are_columns=true) where {FT <: Real} = params_are_columns ? new{FT}(reshape(distribution_samples,1,:)) : new{FT}(reshape(distribution_samples,:,1))
     
@@ -48,7 +50,7 @@ end
 # For the transforms
 abstract type ConstraintType end
 """
-    Constraint <: ConstraintType
+    struct Constraint <: ConstraintType
 
 Contains two functions to map between constrained and unconstrained spaces.
 """
@@ -60,7 +62,7 @@ end
 
 
 """
-    no_constraint()
+    function no_constraint()
 
 Constructs a Constraint with no constraints, enforced by maps x -> x and x -> x.
 """
@@ -71,7 +73,7 @@ function no_constraint()
 end
     
 """
-    bounded_below(lower_bound::FT) where {FT <: Real}
+    function bounded_below(lower_bound::FT) where {FT <: Real}
 
 Constructs a Constraint with provided lower bound, enforced by maps x -> log(x - lower_bound) and x -> exp(x) + lower_bound.
 """
@@ -82,7 +84,7 @@ function bounded_below(lower_bound::FT) where {FT <: Real}
 end
 
 """
-    bounded_above(upper_bound::FT) where {FT <: Real} 
+    function bounded_above(upper_bound::FT) where {FT <: Real} 
 
 Constructs a Constraint with provided upper bound, enforced by maps x -> log(upper_bound - x) and x -> upper_bound - exp(x).
 """
@@ -94,7 +96,7 @@ end
     
 
 """
-    Bounded{FT <: Real} <: ConstraintType
+    function bounded(lower_bound::FT, upper_bound::FT) where {FT <: Real} 
 
 Constructs a Constraint with provided upper and lower bounds, enforced by maps
 x -> log((x - lower_bound) / (upper_bound - x))
@@ -112,7 +114,7 @@ function bounded(lower_bound::FT, upper_bound::FT) where {FT <: Real}
 end
 
 """
-    len(c::Array{CType})
+    function len(c::Array{CType})
 
 The number of constraints, each constraint has length 1.
 """
@@ -125,7 +127,7 @@ function len(carray::Array{CType}) where {CType <: ConstraintType}
 end
 
 """
-    dimension(d<:ParametrizedDistributionType)
+    function get_dimensions(d<:ParametrizedDistributionType)
 
 The number of dimensions of the parameter space
 """
@@ -138,13 +140,16 @@ function dimension(d::Samples)
 end
 
 """
-    n_samples(d::Samples)
+    function n_samples(d::Samples)
 
 The number of samples in the array
 """
 function n_samples(d::Samples)
     return size(d.distribution_samples)[2]
 end
+function n_samples(d::Parameterized)
+    return "Distribution stored in Parameterized form, draw samples using `sample_distribution` function" 
+end
 
 """
     struct ParameterDistribution
@@ -163,10 +168,11 @@ struct ParameterDistribution{PDType <: ParameterDistributionType, CType <: Const
                                                                        ST <: AbstractString}
         
         parameter_distributions = isa(parameter_distributions, PDType) ? [parameter_distributions] : parameter_distributions
+        n_parameter_per_dist = [dimension(pd) for pd in parameter_distributions]
+        
         constraints = isa(constraints, Union{<:ConstraintType,Array{<:ConstraintType}}) ? [constraints] : constraints #to calc n_constraints_per_dist
         names = isa(names, ST) ? [names] : names
             
-        n_parameter_per_dist = [dimension(pd) for pd in parameter_distributions]
         n_constraints_per_dist = [len(c) for c in constraints]
         n_dists = length(parameter_distributions)
         n_names = length(names)        
@@ -188,7 +194,7 @@ end
 ## Functions
 
 """
-    get_name(pd::ParameterDistribution)
+    function get_name(pd::ParameterDistribution)
 
 Returns a list of ParameterDistribution names
 """
@@ -197,7 +203,54 @@ function get_name(pd::ParameterDistribution)
 end
 
 """
-    get_distribution(pd::ParameterDistribution)
+    function get_dimensions(pd::ParameterDistribution)
+
+The number of dimensions of the parameter space
+"""
+function get_dimensions(pd::ParameterDistribution)
+    return [dimension(d) for d in pd.distributions]
+end
+function get_total_dimension(pd::ParameterDistribution)
+    return sum(dimension(d) for d in pd.distributions)
+end
+
+"""
+    function get_n_samples(pd::ParameterDistribution)
+
+The number of samples in a Samples distribution
+"""
+function get_n_samples(pd::ParameterDistribution)
+    return Dict{String,Any}(pd.names[i] => n_samples(d) for (i,d) in enumerate(pd.distributions))
+end
+
+"""
+    function get_all_constraints(pd::ParameterDistribution)
+
+returns the (flattened) array of constraints of the parameter distribution
+"""
+function get_all_constraints(pd::ParameterDistribution)
+    return pd.constraints
+end
+
+"""
+    function batch(pd:ParameterDistribution)
+
+Returns a list of contiguous [collect(1:i), collect(i+1:j),... ] used to split parameter arrays by distribution dimensions
+"""
+function batch(pd::ParameterDistribution)
+    #chunk xarray to give to the different distributions.
+    d_dim = get_dimensions(pd) #e.g [4,1,2]
+    d_dim_tmp = Array{Int64}(undef,size(d_dim)[1]+1)
+    d_dim_tmp[1] = 0
+    for i = 2:size(d_dim)[1]+1
+        d_dim_tmp[i] = sum(d_dim[1:i-1]) # e.g [0,4,5,7]
+    end
+    
+    return [collect(d_dim_tmp[i]+1:d_dim_tmp[i+1]) for i = 1:size(d_dim)[1]] # e.g [1:4, 5:5, 6:7]
+end
+
+"""
+    function get_distribution(pd::ParameterDistribution)
 
 Returns a `Dict` of `ParameterDistribution` distributions by name, (unless sample type)
 """
@@ -215,7 +268,7 @@ end
 """
     function sample_distribution(pd::ParameterDistribution)
 
-Draws samples from the parameter distributions
+Draws samples from the parameter distributions returns an array, with parameters as columns
 """
 function sample_distribution(pd::ParameterDistribution)
     return sample_distribution(pd,1)
@@ -228,99 +281,129 @@ end
 function sample_distribution(d::Samples, n_draws::IT)  where {IT <: Integer}
     n_stored_samples = n_samples(d)
     samples_idx = StatsBase.sample(collect(1:n_stored_samples), n_draws)
-    return d.distribution_samples[:,samples_idx]
-
+    if dimension(d) == 1
+        return reshape(d.distribution_samples[:,samples_idx],:,n_draws) #columns are parameters
+    else
+        return d.distribution_samples[:,samples_idx]
+    end
 end
 
 function sample_distribution(d::Parameterized, n_draws::IT) where {IT <: Integer}
-    return rand(d.distribution, n_draws)
+    if dimension(d) == 1
+        return reshape(rand(d.distribution,n_draws), :, n_draws) #columns are parameters
+    else
+        return rand(d.distribution, n_draws)
+    end
 end
 
-
-#apply transforms
-
 """
-    transform_constrained_to_unconstrained(pd::ParameterDistribution, x::Array{Real})
+    function logpdf(pd::ParameterDistribution, xarray::Array{<:Real,1})
 
-Apply the transformation to map (possibly constrained) parameters `xarray` into the unconstrained space
+Obtains the independent logpdfs of the parameter distributions at xarray (non-Samples Distributions only), and returns their sum.
 """
-function transform_constrained_to_unconstrained(pd::ParameterDistribution, xarray::Array{FT}) where {FT <: Real}
-    #split xarray into chunks 
-#    return cat([transform_constrained_to_unconstrained(c,xarray[i]) for (i,c) in enumerate(pd.constraints)]...,dims=1)
-    return cat([c.constrained_to_unconstrained(xarray[i]) for (i,c) in enumerate(pd.constraints)]...,dims=1)
+function get_logpdf(d::Parameterized, xarray::Array{FT,1}) where {FT <: Real}
+    return logpdf.(d.distribution, xarray)
+end
+
+function get_logpdf(pd::ParameterDistribution, xarray::Array{FT,1}) where {FT <: Real}
+    #first check we don't have sampled distribution
+    for d in pd.distributions
+        if typeof(d) <: Samples
+            throw(ErrorException("Cannot compute get_logpdf of Samples distribution. Consider using a Parameterized type for your prior."))
+        end
+    end
+    #assert xarray correct dim/length
+    if size(xarray)[1] != get_total_dimension(pd)
+        throw(DimensionMismatch("xarray must have dimension equal to the parameter space"))
+    end
+    
+    # get the index of xarray chunks to give to the different distributions.
+    batches = batch(pd)
+
+    # perform the logpdf of each of the distributions, and returns their sum    
+    return sum(cat([get_logpdf(d, xarray[batches[i]]) for (i,d) in enumerate(pd.distributions)]...,dims=1))
 end
 
 """
-    transform_unconstrained_to_constrained(pd::ParameterDistribution, xarray::Array{Real})
+    function get_cov(pd::ParameterDistribution)
 
-Apply the transformation to map parameters `xarray` from the unconstrained space into (possibly constrained) space
+returns a blocked covariance of the distributions
 """
-function transform_unconstrained_to_constrained(pd::ParameterDistribution, xarray::Array{FT}) where {FT <: Real}
-#    return [transform_unconstrained_to_constrained(c,xarray[i]) for (i,c) in enumerate(pd.constraints)]
-    return cat([c.unconstrained_to_constrained(xarray[i]) for (i,c) in enumerate(pd.constraints)]...,dims=1)
+function get_cov(d::Parameterized)
+    return cov(d.distribution)
 end
 
+function get_cov(d::Samples)
+    return cov(d.distribution_samples, dims=2) #parameters are columns
+end
 
+function get_var(d::Parameterized)
+    return var(d.distribution)
+end
 
+function get_var(d::Samples)
+    return var(d.distribution_samples)
+end
 
-# """
-# No constraint mapping x -> x
-# """
-# function transform_constrained_to_unconstrained(c::NoConstraint , x::FT) where {FT <: Real}
-#     return x
-# end
 
-# """
-# Bounded below -> unbounded, use mapping x -> log(x - lower_bound)
-# """
-# function transform_constrained_to_unconstrained(c::BoundedBelow, x::FT) where {FT <: Real}    
-#     return log(x - c.lower_bound)
-# end
+function get_cov(pd::ParameterDistribution)
+    #first check we don't have sampled distribution
+    
+    d_dims = get_dimensions(pd)
+    
+    # create each block (co)variance
+    block_cov = Array{Any}(undef,size(d_dims)[1]) 
+    for (i,dimension) in enumerate(d_dims)
+        if dimension == 1
+            block_cov[i] = get_var(pd.distributions[i]) 
+        else
+            block_cov[i] = get_cov(pd.distributions[i])
+        end
+    end
+   
+    return cat(block_cov...,dims=(1,2)) #build the block diagonal (dense) matrix
+    
+end
+"""
+    function get_mean(pd::ParameterDistribution)
 
-# """
-# Bounded above -> unbounded, use mapping x -> log(upper_bound - x)
-# """
-# function transform_constrained_to_unconstrained(c::BoundedAbove, x::FT) where {FT <: Real}    
-#     return log(c.upper_bound - x)
-# end
+returns a mean of the distirbutions
+"""
 
-# """
-# Bounded -> unbounded, use mapping x -> log((x - lower_bound) / (upper_bound - x)
-# """
-# function transform_constrained_to_unconstrained(c::Bounded, x::FT) where {FT <: Real}    
-#     return log( (x - c.lower_bound) / (c.upper_bound - x))
-# end
+function get_mean(d::Parameterized)
+    return mean(d.distribution)
+end
 
+function get_mean(d::Samples)
+    return mean(d.distribution_samples, dims=2) #parameters are columns
+end
 
+function get_mean(pd::ParameterDistribution)
+    return reshape(cat([get_mean(d) for d in pd.distributions]...,dims=1),:,1)
+end
 
 
-# """
-# No constraint mapping x -> x
-# """
-# function transform_unconstrained_to_constrained(c::NoConstraint , x::FT) where {FT <: Real}
-#     return x
-# end
+#apply transforms
+
+"""
+    function transform_constrained_to_unconstrained(pd::ParameterDistribution, x::Array{<:Real})
 
-# """
-# Unbounded -> bounded below, use mapping x -> exp(x) + lower_bound
-# """
-# function transform_unconstrained_to_constrained(c::BoundedBelow, x::FT) where {FT <: Real}    
-#     return exp(x) + c.lower_bound
-# end
+Apply the transformation to map (possibly constrained) parameters `xarray` into the unconstrained space
+"""
+function transform_constrained_to_unconstrained(pd::ParameterDistribution, xarray::Array{FT}) where {FT <: Real}
+    return cat([c.constrained_to_unconstrained(xarray[i]) for (i,c) in enumerate(pd.constraints)]...,dims=1)
+end
+
+"""
+    function transform_unconstrained_to_constrained(pd::ParameterDistribution, xarray::Array{Real})
+
+Apply the transformation to map parameters `xarray` from the unconstrained space into (possibly constrained) space
+"""
+function transform_unconstrained_to_constrained(pd::ParameterDistribution, xarray::Array{FT}) where {FT <: Real}
+    return cat([c.unconstrained_to_constrained(xarray[i]) for (i,c) in enumerate(pd.constraints)]...,dims=1)
+end
 
-# """
-# Unbounded -> bounded above, use mapping x -> upper_bound - exp(x)
-# """
-# function transform_unconstrained_to_constrained(c::BoundedAbove, x::FT) where {FT <: Real}    
-#     return c.upper_bound - exp(x)
-# end
 
-# """
-# Unbounded -> bounded, use mapping x -> (upper_bound * exp(x) + lower_bound) / (exp(x) + 1)
-# """
-# function transform_unconstrained_to_constrained(c::Bounded, x::FT) where {FT <: Real}    
-#     return (c.upper_bound * exp(x) + c.lower_bound) / (exp(x) + 1)
-# end
 
 
 
diff --git a/src/Priors.jl b/src/Priors.jl
deleted file mode 100644
index c634b00fb..000000000
--- a/src/Priors.jl
+++ /dev/null
@@ -1,66 +0,0 @@
-module Priors
-
-using Distributions
-using DocStringExtensions
-using Random
-
-export Prior
-export Deterministic
-
-"""
-    Deterministic(value::FT)
-
-A deterministic distribution is a distribution over only one value, which it takes with probability 1.
-It can be used when a distribution is required, but the outcome is deterministic. 
-```
-d = Deterministic(3.0)  # Deterministic "distribution" with value 3.0
-rand(d)                 # Always returns the same value
-```
-"""
-struct Deterministic{FT<:Real} <: DiscreteUnivariateDistribution
-    value::FT
-end
-
-eltype(::Deterministic{FT}) where {FT} = FT
-rand(rng::AbstractRNG, d::Deterministic) = d.value
-pdf(d::Deterministic{FT}, x::FT) where {FT} = x == d.value ? FT(1) : FT(0)
-cdf(d::Deterministic{FT}, x::FT) where {FT} = x < d.value ? FT(0) : FT(1)
-quantile(d::Deterministic{FT}, p) where {FT} = 0 <= p <= 1 ? d.value : FT(NaN)
-minimum(d::Deterministic) = d.value
-maximum(d::Deterministic) = d.value
-insupport(d::Deterministic, x) = x == d.value
-mean(d::Deterministic) = d.value
-var(d::Deterministic{FT}) where {FT} = FT(0)
-mode(d::Deterministic) = d.value
-skewness(d::Deterministic{FT}) where {FT} = FT(0)
-kurtosis(d::Deterministic{FT}) where {FT} = FT(0)
-entropy(d::Deterministic{FT}) where {FT} = FT(0)
-mgf(d::Deterministic, t) = exp(t * d.value)
-cf(d::Deterministic, t) = exp(im * t * d.value)
-
-
-"""
-    struct Prior
-
-Structure to hold information on the prior distribution of a UQ parameter
-
-# Fields
-$(DocStringExtensions.FIELDS)
-"""
-struct Prior
-    "prior distribution"
-    dist::Union{Distribution, Deterministic}
-    "name of the parameter with prior distribution dist"
-    param_name::String
-end
-
-
-"""
-    Prior(value::FT, param_name::String)
-"""
-function Prior(value::FT, param_name::String) where {FT<:Real}
-    Prior(Deterministic(value), param_name)
-end
-
-
-end # module 
diff --git a/test/EKP/runtests.jl b/test/EKP/runtests.jl
index 0db6db9ae..ac4eca61c 100644
--- a/test/EKP/runtests.jl
+++ b/test/EKP/runtests.jl
@@ -4,8 +4,7 @@ using Random
 using Test
 
 using CalibrateEmulateSample.EKP
-using CalibrateEmulateSample.Priors
-
+using CalibrateEmulateSample.ParameterDistributionStorage
 @testset "EKP" begin
 
     # Seed for pseudo-random number generator
@@ -39,13 +38,16 @@ using CalibrateEmulateSample.Priors
     @test norm(y_star - G(u_star)) < n_obs * noise_level^2
 
     #### Define prior information on parameters
-    param_names = ["u1", "u2"]
-    priors     = [Priors.Prior(Normal(1., sqrt(2)), param_names[1]),
-                  Priors.Prior(Normal(1., sqrt(2)), param_names[2])]
+    prior_distns = [Parameterized(Normal(1., sqrt(2))),
+              Parameterized(Normal(1., sqrt(2)))]
+    constraints = [[no_constraint()], [no_constraint()]]
+    prior_names = ["u1","u2"]
+    prior = ParameterDistribution(prior_distns, constraints, prior_names)
+    
+    prior_mean = reshape(get_mean(prior),:)
 
-    prior_mean = [1., 1.]
     # Assuming independence of u1 and u2
-    prior_cov = convert(Array, Diagonal([sqrt(2.), sqrt(2.)]))
+    prior_cov = get_cov(prior)#convert(Array, Diagonal([sqrt(2.), sqrt(2.)]))
 
     ###
     ###  Calibrate (1): Ensemble Kalman Inversion
@@ -53,11 +55,11 @@ using CalibrateEmulateSample.Priors
 
     N_ens = 50 # number of ensemble members
     N_iter = 20 # number of EKI iterations
-    initial_ensemble = EKP.construct_initial_ensemble(N_ens, priors;
+    initial_ensemble = EKP.construct_initial_ensemble(prior,N_ens;
                                                     rng_seed=rng_seed)
     @test size(initial_ensemble) == (N_ens, n_par)
 
-    ekiobj = EKP.EKObj(initial_ensemble, param_names,
+    ekiobj = EKP.EKObj(initial_ensemble,
                        y_obs, Γy, Inversion())
 
     # EKI iterations
@@ -88,7 +90,7 @@ using CalibrateEmulateSample.Priors
     ###
     ###  Calibrate (2): Ensemble Kalman Sampleer
     ###
-    eksobj = EKP.EKObj(initial_ensemble, param_names,
+    eksobj = EKP.EKObj(initial_ensemble,
                       y_obs, Γy,
                       Sampler(prior_mean, prior_cov))
 
diff --git a/test/MCMC/runtests.jl b/test/MCMC/runtests.jl
index a630548b9..7676a1803 100644
--- a/test/MCMC/runtests.jl
+++ b/test/MCMC/runtests.jl
@@ -5,7 +5,8 @@ using GaussianProcesses
 using Test
 
 using CalibrateEmulateSample.MCMC
-using CalibrateEmulateSample.Priors
+#using CalibrateEmulateSample.Priors
+using CalibrateEmulateSample.ParameterDistributionStorage
 using CalibrateEmulateSample.GPEmulator
 
 @testset "MCMC" begin
@@ -39,7 +40,12 @@ using CalibrateEmulateSample.GPEmulator
     ### Define prior
     umin = -1.0
     umax = 6.0
-    prior = [Priors.Prior(Uniform(umin, umax), "u")]  # prior on u
+    #prior = [Priors.Prior(Uniform(umin, umax), "u")]  # prior on u
+    prior_dist = Parameterized(Normal(0,1))
+    prior_constraint = bounded(-1.0,6.0)
+    prior_name = "u"
+    prior = ParameterDistribution(prior_dist, prior_constraint, prior_name)
+    
     obs_sample = [1.0]
 
     # MCMC parameters    
@@ -62,9 +68,9 @@ using CalibrateEmulateSample.GPEmulator
     mcmc = MCMCObj(obs_sample, σ2_y, prior, step, param_init, max_iter, 
                    mcmc_alg, burnin, svdflag=true)
     sample_posterior!(mcmc, gpobj, max_iter)
-    posterior = get_posterior(mcmc)      
-    post_mean = mean(posterior, dims=1)[1]
-    
+    posterior_distribution = get_posterior(mcmc)      
+    #post_mean = mean(posterior, dims=1)[1]
+    posterior_mean = get_mean(posterior_distribution)
     # We had svdflag=true, so the MCMCObj stores a transformed sample, 
     # 1.0/sqrt(0.05) * obs_sample ≈ 4.472
     @test mcmc.obs_sample ≈ [4.472] atol=1e-2
@@ -72,9 +78,10 @@ using CalibrateEmulateSample.GPEmulator
     @test mcmc.burnin == burnin
     @test mcmc.algtype == mcmc_alg
     @test mcmc.iter[1] == max_iter + 1
-    @test size(posterior) == (max_iter - burnin + 1, length(param_init))
+    @test get_n_samples(posterior_distribution)[prior_name] == max_iter - burnin + 1
+    @test get_total_dimension(posterior_distribution) == length(param_init)
     @test_throws Exception MCMCObj(obs_sample, σ2_y, prior, step, param_init, 
                                    max_iter, "gibbs", burnin)
-    @test post_mean ≈ π/2 atol=4e-1
+    @test isapprox(posterior_mean[1] - π/2, 0.0; atol=4e-1)
 
 end
diff --git a/test/ParameterDistribution/runtests.jl b/test/ParameterDistribution/runtests.jl
index f5fa154b9..a4780d613 100644
--- a/test/ParameterDistribution/runtests.jl
+++ b/test/ParameterDistribution/runtests.jl
@@ -93,10 +93,9 @@ using CalibrateEmulateSample.ParameterDistributionStorage
         @test u.distributions == [d1,d2]
         @test u.constraints == cat([c1,c2]...,dims=1)
         @test u.names == [name1,name2]
-
     end
 
-    @testset "get/sample functions" begin
+    @testset "getter functions" begin
         # setup for the tests:
         d1 = Parameterized(MvNormal(4,0.1))
         c1 = [no_constraint(),
@@ -112,10 +111,22 @@ using CalibrateEmulateSample.ParameterDistributionStorage
         u2 = ParameterDistribution(d2,c2,name2)
         
         u = ParameterDistribution([d1,d2], [c1,c2], [name1,name2])
+
+        # Test for get_dimension(s)
+        @test get_total_dimension(u1) == 4
+        @test get_total_dimension(u2) == 1
+        @test get_total_dimension(u) == 5
+        @test get_dimensions(u1) == [4]
+        @test get_dimensions(u2) == [1]
+        @test get_dimensions(u) == [4,1]
         
         # Tests for get_name
         @test get_name(u1) == [name1]
         @test get_name(u) == [name1, name2]
+
+        # Tests for get_n_samples
+        @test typeof(get_n_samples(u)[name1]) <: String
+        @test get_n_samples(u)[name2] == 4
         
         # Tests for get_distribution
         @test get_distribution(d1) == MvNormal(4,0.1)
@@ -127,6 +138,39 @@ using CalibrateEmulateSample.ParameterDistributionStorage
         @test d[name1] == MvNormal(4,0.1)
         @test typeof(d[name2]) <: String
 
+        # Test for get_all_constraints
+        @test get_all_constraints(u) == cat([c1,c2]...,dims=1)
+    end
+
+    @testset "statistics functions" begin
+
+        # setup for the tests:
+        d1 = Parameterized(MvNormal(4,0.1))
+        c1 = [no_constraint(),
+              bounded_below(-1.0),
+              bounded_above(0.4),
+              bounded(-0.1,0.2)]
+        name1 = "constrained_mvnormal"
+        u1 = ParameterDistribution(d1,c1,name1)
+
+        d2 = Samples([1 2 3 4])
+        c2 = [bounded(10,15)]
+        name2 = "constrained_sampled"
+        u2 = ParameterDistribution(d2,c2,name2)
+        
+        d3 = Parameterized(Beta(2,2))
+        c3 = [no_constraint()]
+        name3 = "unconstrained_beta"
+        u3 = ParameterDistribution(d3,c3,name3)
+
+        u = ParameterDistribution([d1,d2],[c1,c2],[name1,name2])
+
+        d4 = Samples([1 2 3 4 5 6 7 8; 8 7 6 5 4 3 2 1])
+        c4 = [no_constraint(),
+              no_constraint()]
+        name4 = "constrained_MVsampled"
+        v = ParameterDistribution([d1,d2,d3,d4],[c1,c2,c3,c4],[name1,name2,name3,name4])
+        
         # Tests for sample distribution
         seed=2020
         Random.seed!(seed)
@@ -157,6 +201,23 @@ using CalibrateEmulateSample.ParameterDistributionStorage
         Random.seed!(seed)
         s = sample_distribution(u,3)
         @test s == cat([s1,s2]...,dims=1)
+
+        #Test for get_logpdf
+        @test_throws ErrorException get_logpdf(u,zeros(get_total_dimension(u))) 
+        x_in_bd = [0.5]
+        Random.seed!(seed)
+        lpdf3 = logpdf.(Beta(2,2),x_in_bd)[1] #throws deprecated warning without "."
+        Random.seed!(seed)
+        @test isapprox(get_logpdf(u3,x_in_bd) - lpdf3 , 0.0; atol=1e-6)
+        @test_throws DimensionMismatch get_logpdf(u3, [0.5,0.5])
+
+        #Test for get_cov, get_var        
+        block_cov = cat([get_cov(d1),get_var(d2),get_var(d3),get_cov(d4)]..., dims=(1,2)) 
+        @test isapprox(get_cov(v) - block_cov, zeros(get_total_dimension(v),get_total_dimension(v)); atol=1e-6)
+        #Test for get_mean
+        means = reshape(cat([get_mean(d1), get_mean(d2), get_mean(d3), get_mean(d4)]...,dims=1),:,1)
+        @test isapprox(get_mean(v) - means, zeros(get_total_dimension(v)); atol=1e-6)
+        
     end
 
     @testset "transform functions" begin
diff --git a/test/Priors/runtests.jl b/test/Priors/runtests.jl
deleted file mode 100644
index 66585153d..000000000
--- a/test/Priors/runtests.jl
+++ /dev/null
@@ -1,21 +0,0 @@
-using Test
-using Distributions
-
-using CalibrateEmulateSample.Priors
-
-@testset "Priors" begin
-
-    # u1 is a parameter with a Gamma prior
-    u1_prior = Prior(Gamma(2.0, 0.8), "u1")
-    @test u1_prior.param_name == "u1"
-    @test u1_prior.dist == Gamma(2.0, 0.8)
-
-    # u2 is a parameter with a Deterministic prior, i.e., it always takes on
-    # the same value, val
-    val = 3.0
-    u2_prior = Prior(Deterministic(val), "u2")
-    @test u2_prior.param_name == "u2"
-    @test u2_prior.dist == Deterministic(val)
-    @test u2_prior.dist.value == val
-
-end
diff --git a/test/runtests.jl b/test/runtests.jl
index ed3c4f384..97c2388d9 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -24,8 +24,7 @@ end
         end
     end
 
-    for submodule in ["Priors",
-		      "EKP",
+    for submodule in ["EKP",
 		      "GPEmulator",
 		      "MCMC",
 		      "Observations",