diff --git a/.github/workflows/Breakage.yml b/.github/workflows/Breakage.yml
index f9bf7970..0cafac23 100644
--- a/.github/workflows/Breakage.yml
+++ b/.github/workflows/Breakage.yml
@@ -15,8 +15,8 @@ jobs:
       fail-fast: false
       matrix:
         pkg: [
-          "cscherrer/Soss.jl",
-          "mschauer/Mitosis.jl"
+          # "cscherrer/Soss.jl",
+          # "mschauer/Mitosis.jl"
         ]
         pkgversion: [latest, stable]
 
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 193fa517..4876ceef 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -22,8 +22,8 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.6'
-          - '1.7'
+          - '1.8'
+          - '1.9'
         os:
           - ubuntu-latest
           - macOS-latest
diff --git a/Project.toml b/Project.toml
index 7dfa46d1..1437c744 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "MeasureTheory"
 uuid = "eadaa1a4-d27c-401d-8699-e962e1bbc33b"
 authors = ["Chad Scherrer <chad.scherrer@gmail.com> and contributors"]
-version = "0.18.4"
+version = "0.19.0"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
@@ -10,16 +10,15 @@ Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 DensityInterface = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
-DistributionMeasures = "35643b39-bfd4-4670-843f-16596ca89bf3"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DynamicIterators = "6c76993d-992e-5bf1-9e63-34920a5a5a38"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 Infinities = "e1ba4f0e-776d-440f-acd9-e1d2e9742647"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
 KeywordCalls = "4d827475-d3e4-43d6-abe3-9688362ede9f"
-LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
 MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
@@ -47,30 +46,29 @@ Compat = "3.42, 4"
 ConcreteStructs = "0.2"
 ConstructionBase = "1.3"
 DensityInterface = "0.4"
-DistributionMeasures = "0.2"
 Distributions = "0.25"
 DynamicIterators = "0.4"
-FillArrays = "0.12, 0.13"
+FillArrays = "1"
+ForwardDiff = "0.10"
 IfElse = "0.1"
 Infinities = "0.1"
 InverseFunctions = "0.1"
 KeywordCalls = "0.2"
-LazyArrays = "0.22"
 LogExpFunctions = "0.3.3"
 MLStyle = "0.4"
 MacroTools = "0.5"
 MappedArrays = "0.4"
-MeasureBase = "0.13"
+MeasureBase = "0.14"
 NamedTupleTools = "0.13, 0.14"
 PositiveFactorizations = "0.2"
 PrettyPrinting = "0.3, 0.4"
 Reexport = "1"
 SpecialFunctions = "1, 2"
-Static = "0.5, 0.6"
+Static = "0.8"
 StaticArraysCore = "1"
-StatsBase = "0.32, 0.33"
+StatsBase = "0.34"
 StatsFuns = "0.9, 1"
-TransformVariables = "0.5, 0.6, 0.7"
+TransformVariables = "0.8"
 Tricks = "0.1"
 julia = "1.6"
 
diff --git a/src/MeasureTheory.jl b/src/MeasureTheory.jl
index e43642d7..96b3a7f7 100644
--- a/src/MeasureTheory.jl
+++ b/src/MeasureTheory.jl
@@ -8,12 +8,9 @@ using MLStyle
 import TransformVariables
 const TV = TransformVariables
 
-using DistributionMeasures
 using TransformVariables: asℝ₊, as𝕀, asℝ, transform
 
 import Base
-import Distributions
-const Dists = Distributions
 
 export TV
 export transform
@@ -27,14 +24,12 @@ export Lebesgue
 export ℝ, ℝ₊, 𝕀
 export ⊙
 export SpikeMixture
-export CountingMeasure
 export TrivialMeasure
 export Likelihood
 export testvalue
 export basekernel
 
 using Infinities
-using DynamicIterators
 using KeywordCalls
 using ConstructionBase
 using Accessors
@@ -58,9 +53,12 @@ import MeasureBase:
     paramnames,
     ∫,
     𝒹,
-    ∫exp
+    ∫exp,
+    smf,
+    invsmf,
+    massof
 import MeasureBase: ≪
-using MeasureBase: BoundedInts, BoundedReals, CountingMeasure, IntegerDomain, IntegerNumbers
+using MeasureBase: BoundedInts, BoundedReals, CountingBase, IntegerDomain, IntegerNumbers
 using MeasureBase: weightedmeasure, restrict
 using MeasureBase: AbstractTransitionKernel
 
@@ -97,13 +95,16 @@ using MeasureBase: kernel
 using MeasureBase: Returns
 import MeasureBase: proxy, @useproxy
 import MeasureBase: basemeasure_depth
-using MeasureBase: LebesgueMeasure
+using MeasureBase: LebesgueBase
 
 import DensityInterface: logdensityof
 import DensityInterface: densityof
 import DensityInterface: DensityKind
 using DensityInterface
 
+using ForwardDiff
+using ForwardDiff: Dual
+
 gentype(μ::AbstractMeasure) = typeof(testvalue(μ))
 
 # gentype(μ::AbstractMeasure) = gentype(basemeasure(μ))
@@ -117,20 +118,22 @@ xlogy(x, y) = x * log(y)
 xlog1py(x::Number, y::Number) = LogExpFunctions.xlog1py(x, y)
 xlog1py(x, y) = x * log(1 + y)
 
+using MeasureBase: Φ, Φinv
 as(args...; kwargs...) = TV.as(args...; kwargs...)
 
+# Type piracy until https://github.com/JuliaMath/MeasureBase.jl/issues/127 is fixed
+MeasureBase.rand(::FixedRNG, ::Type{Bool}) = true
+
 include("utils.jl")
 include("const.jl")
 include("combinators/for.jl")
-# include("traits.jl")
-include("parameterized.jl")
 
 include("macros.jl")
 include("combinators/affine.jl")
 include("combinators/weighted.jl")
 include("combinators/product.jl")
 include("combinators/transforms.jl")
-include("combinators/exponential-families.jl")
+# include("combinators/exponential-families.jl")
 include("resettable-rng.jl")
 include("realized.jl")
 include("combinators/chain.jl")
@@ -165,5 +168,7 @@ include("combinators/ifelse.jl")
 include("transforms/corrcholesky.jl")
 include("transforms/ordered.jl")
 
+include("parameterized.jl")
+
 include("distproxy.jl")
 end # module
diff --git a/src/combinators/affine.jl b/src/combinators/affine.jl
index f4ccd80e..291ce628 100644
--- a/src/combinators/affine.jl
+++ b/src/combinators/affine.jl
@@ -128,7 +128,7 @@ logjac(f::AffineTransform{(:μ, :σ)}) = logjac(f.σ)
 logjac(f::AffineTransform{(:μ, :λ)}) = -logjac(f.λ)
 logjac(f::AffineTransform{(:σ,)}) = logjac(f.σ)
 logjac(f::AffineTransform{(:λ,)}) = -logjac(f.λ)
-logjac(f::AffineTransform{(:μ,)}) = 0.0
+logjac(f::AffineTransform{(:μ,)}) = static(0.0)
 
 ###############################################################################
 
@@ -161,7 +161,7 @@ basemeasure(d::OrthoLebesgue) = d
 
 logdensity_def(::OrthoLebesgue, x) = static(0.0)
 
-struct AffinePushfwd{N,M,T} <: AbstractMeasure
+struct AffinePushfwd{N,M,T} <: MeasureBase.AbstractPushforward
     f::AffineTransform{N,T}
     parent::M
 end
@@ -172,6 +172,10 @@ function Pretty.tile(d::AffinePushfwd)
     Pretty.list_layout([pars, Pretty.tile(d.parent)]; prefix = :AffinePushfwd)
 end
 
+@inline MeasureBase.transport_origin(d::AffinePushfwd) = d.parent
+@inline MeasureBase.to_origin(d::AffinePushfwd, y) = inverse(getfield(d, :f))(y)
+@inline MeasureBase.from_origin(d::AffinePushfwd, x) = getfield(d, :f)(x)
+
 @inline function testvalue(d::AffinePushfwd)
     f = getfield(d, :f)
     z = testvalue(parent(d))
@@ -282,7 +286,7 @@ end
 @inline function basemeasure(d::AffinePushfwd{N,L}) where {N,L<:Lebesgue}
     weightedmeasure(-logjac(d), d.parent)
 end
-@inline function basemeasure(d::AffinePushfwd{N,L}) where {N,L<:LebesgueMeasure}
+@inline function basemeasure(d::AffinePushfwd{N,L}) where {N,L<:LebesgueBase}
     weightedmeasure(-logjac(d), d.parent)
 end
 
@@ -313,19 +317,6 @@ supportdim(nt::NamedTuple{(:σ,),T}) where {T}    = colsize(nt.σ)
 supportdim(nt::NamedTuple{(:λ,),T}) where {T}    = rowsize(nt.λ)
 supportdim(nt::NamedTuple{(:μ,),T}) where {T}    = size(nt.μ)
 
-asparams(::AffinePushfwd, ::StaticSymbol{:μ}) = asℝ
-asparams(::AffinePushfwd, ::StaticSymbol{:σ}) = asℝ₊
-asparams(::Type{A}, ::StaticSymbol{:μ}) where {A<:AffinePushfwd} = asℝ
-asparams(::Type{A}, ::StaticSymbol{:σ}) where {A<:AffinePushfwd} = asℝ₊
-
-function asparams(d::AffinePushfwd{N,M,T}, ::StaticSymbol{:μ}) where {N,M,T<:AbstractArray}
-    as(Array, asℝ, size(d.μ))
-end
-
-function asparams(d::AffinePushfwd{N,M,T}, ::StaticSymbol{:σ}) where {N,M,T<:AbstractArray}
-    as(Array, asℝ, size(d.σ))
-end
-
 function Base.rand(rng::Random.AbstractRNG, ::Type{T}, d::AffinePushfwd) where {T}
     z = rand(rng, T, parent(d))
     f = getfield(d, :f)
@@ -336,8 +327,8 @@ end
     insupport(d.parent, inverse(d.f)(x))
 end
 
-@inline function Distributions.cdf(d::AffinePushfwd, x)
-    cdf(parent(d), inverse(d.f)(x))
+@inline function MeasureBase.smf(d::AffinePushfwd, x)
+    smf(parent(d), inverse(d.f)(x))
 end
 
 @inline function mean(d::AffinePushfwd)
@@ -364,3 +355,67 @@ end
 @inline function std(d::AffinePushfwd{(:μ, :λ)})
     std(parent(d)) / d.λ
 end
+
+###############################################################################
+# smf
+
+@inline function smf(d::AffinePushfwd{(:μ,)}, x)
+    f = getfield(d, :f)
+    smf(parent(d), inverse(f)(x))
+end
+
+@inline function smf(d::AffinePushfwd{(:μ, :σ)}, x)
+    f = getfield(d, :f)
+    p = smf(parent(d), inverse(f)(x))
+    d.σ > 0 ? p : one(p) - p
+end
+
+@inline function smf(d::AffinePushfwd{(:σ,)}, x)
+    f = getfield(d, :f)
+    p = smf(parent(d), inverse(f)(x))
+    d.σ > 0 ? p : one(p) - p
+end
+
+@inline function smf(d::AffinePushfwd{(:λ,)}, x)
+    f = getfield(d, :f)
+    p = smf(parent(d), inverse(f)(x))
+    d.λ > 0 ? p : one(p) - p
+end
+
+@inline function smf(d::AffinePushfwd{(:μ, :λ)}, x)
+    f = getfield(d, :f)
+    p = smf(parent(d), inverse(f)(x))
+    d.λ > 0 ? p : one(p) - p
+end
+
+###############################################################################
+# invsmf
+
+@inline function invsmf(d::AffinePushfwd{(:μ,)}, p)
+    f = getfield(d, :f)
+    f(invsmf(parent(d), p))
+end
+
+@inline function invsmf(d::AffinePushfwd{(:μ, :σ)}, p)
+    p = d.σ > 0 ? p : one(p) - p
+    f = getfield(d, :f)
+    f(invsmf(parent(d), p))
+end
+
+@inline function invsmf(d::AffinePushfwd{(:σ,)}, p)
+    p = d.σ > 0 ? p : one(p) - p
+    f = getfield(d, :f)
+    f(invsmf(parent(d), p))
+end
+
+@inline function invsmf(d::AffinePushfwd{(:λ,)}, p)
+    p = d.λ > 0 ? p : one(p) - p
+    f = getfield(d, :f)
+    f(invsmf(parent(d), p))
+end
+
+@inline function invsmf(d::AffinePushfwd{(:μ, :λ)}, p)
+    p = d.λ > 0 ? p : one(p) - p
+    f = getfield(d, :f)
+    f(invsmf(parent(d), p))
+end
diff --git a/src/combinators/exponential-families.jl b/src/combinators/exponential-families.jl
deleted file mode 100644
index e3e169a7..00000000
--- a/src/combinators/exponential-families.jl
+++ /dev/null
@@ -1,119 +0,0 @@
-export ExponentialFamily
-using LazyArrays
-
-@concrete terse struct ExponentialFamily <: AbstractTransitionKernel
-    support_contains
-    base
-    mdim
-    pdim
-    t
-    x
-    a
-end
-
-MeasureBase.insupport(fam::ExponentialFamily, x) = fam.support_contains(x)
-
-function ExponentialFamily(support_contains, base, mdim, pdim, t, a)
-    return ExponentialFamily(support_contains, base, mdim, pdim, t, I, a)
-end
-
-function MeasureBase.powermeasure(fam::ExponentialFamily, dims::NTuple)
-    support_contains(x) = all(xj -> fam.support_contains(xj), x)
-    t = Tuple((y -> f.(y) for f in fam.t))
-    a(η) = LazyArrays.BroadcastArray(fam.a, η)
-    p = prod(dims)
-    ExponentialFamily(
-        support_contains,
-        fam.base^dims,
-        fam.mdim * p,
-        fam.pdim * p,
-        t,
-        fam.x,
-        a,
-    )
-end
-
-powermeasure(fam::ExponentialFamily, ::Tuple{}) = fam
-
-@concrete terse struct ExpFamMeasure <: AbstractMeasure
-    fam
-    η # instantiated to a value
-    a # instantiated to a value
-end
-
-MeasureBase.insupport(μ::ExpFamMeasure, x) = μ.fam.support_contains(x)
-
-@inline function (fam::ExponentialFamily)(β)
-    η = fam.x * β
-    a = fam.a(η)
-    ExpFamMeasure(fam, η, a)
-end
-
-MeasureBase.basemeasure(d::ExpFamMeasure) = d.fam.base
-
-tracedot(a::AbstractVector, b::AbstractVector) = dot(a, b)
-
-tracedot(a::AbstractVector, x, b::AbstractVector) = dot(a, x, b)
-
-tracedot(a, b) = sum((dot(view(a, :, j), view(b, :, j)) for j in 1:size(a, 2)))
-
-tracedot(a, x, b) =
-    sum(1:size(a, 2)) do j
-        dot(view(a, :, j), x, view(b, :, j))
-    end
-
-# @inline function tracedot(a::BlockDiag, b::BlockDiag)
-#     numblocks = length(a.blocks)
-#     sum(tracedot(a.blocks[j], b.blocks[j]) for j in 1:length(a.blocks))
-# end
-
-# @inline function tracedot(a::BlockDiag, x::BlockDiag, b::BlockDiag)
-#     numblocks = length(x.blocks)
-#     sum(tracedot(a.blocks[j], x.blocks[j], b.blocks[j]) for j in 1:length(x.blocks))
-# end
-
-function logdensity_def(d::ExpFamMeasure, y)
-    t = ApplyArray(vcat, (f.(y) for f in d.fam.t)...)
-    η = d.η
-    dot(t, η)
-end
-
-function withX(fam::ExponentialFamily, x)
-    @inline t(y) = fam.t.(y)
-    newx = ApplyArray(kron, x, fam.x)
-    η(β) = fam.η.(β)
-    a(β) = sum(fam.a, β)
-    ExponentialFamily(fam.base^size(x, 1), t, x, η, a)
-end
-
-@concrete terse struct ExpFamLikelihood <: AbstractLikelihood
-    fam
-    y
-    tᵀx
-    c
-end
-
-# function regression(fam, uᵀ, vᵀ)
-
-# end
-
-function MeasureBase.likelihoodof(fam::ExponentialFamily, y)
-    c = logdensityof(fam.base, y)
-    t = ApplyArray(vcat, (f.(y) for f in fam.t)...)
-    tᵀx = t' * fam.x
-    ExpFamLikelihood(fam, y, tᵀx, c)
-end
-
-@inline function logdensityof(ℓ::ExpFamLikelihood, β)
-    xβ = ApplyArray(*, ℓ.fam.x, β)
-    a = sum(ℓ.fam.a(xβ))
-    # a = sum(ℓ.fam.a, ApplyArray(*, ℓ.fam.uᵀ', ℓ.fam.vᵀ, β))
-    ℓ.c + dot(ℓ.tᵀx, β) - a
-end
-
-basemeasure(fam::ExponentialFamily) = fam.base
-
-# function stack_functions(funs, inds)
-#     function(x::AbstractArray{T,N}) where {T,N}
-#         ApplyArray(cat, )
-# end
diff --git a/src/combinators/for.jl b/src/combinators/for.jl
index aa95206a..299a51e1 100644
--- a/src/combinators/for.jl
+++ b/src/combinators/for.jl
@@ -69,6 +69,17 @@ as(d::For) = as(Array, as(first(marginals(d))), size(first(d.inds))...)
     end
 end
 
+@inline function logdensity_def(
+    d::For{T,F,I},
+    x::AbstractVector{X};
+) where {X,T,F,I<:Tuple{<:AbstractUnitRange}}
+    ind = only(d.inds)
+    Δj = firstindex(x) - first(ind)
+    sum(ind) do j
+        @inbounds logdensity_def(d.f(j), x[j+Δj])
+    end
+end
+
 function logdensity_def(d::For, x::AbstractVector)
     get_i(j) = map(Base.Fix2(getindex, j), d.inds)
     # get_i(j) =  (getindex(ind, j) for ind in d.inds)
@@ -104,6 +115,15 @@ function basemeasure(d::For{T,F,I}) where {T,F,I}
     _basemeasure(d, B, sing)
 end
 
+@inline function _basemeasure(
+    d::For{T,F,I},
+    ::Type{<:WeightedMeasure{StaticFloat64{0.0},B}},
+    ::True,
+) where {T,F,I,B}
+    dim = size(first(d.inds))
+    instance(B)^dim
+end
+
 @inline function _basemeasure(
     d::For{T,F,I},
     ::Type{<:WeightedMeasure{StaticFloat64{N},B}},
diff --git a/src/combinators/ifelse.jl b/src/combinators/ifelse.jl
index a0cb89ec..3157601b 100644
--- a/src/combinators/ifelse.jl
+++ b/src/combinators/ifelse.jl
@@ -4,7 +4,9 @@ struct IfElseMeasure{B,T,F} <: AbstractMeasure
     f::F
 end
 
-insupport(d::IfElseMeasure, x) = insupport(d.t, x) || insupport(d.f, x)
+using MeasureBase: istrue
+
+insupport(d::IfElseMeasure, x) = istrue(insupport(d.t, x)) || istrue(insupport(d.f, x))
 
 function logdensity_def(d::IfElseMeasure, x)
     p = mean(d.b)
diff --git a/src/combinators/product.jl b/src/combinators/product.jl
index f0c03e0b..2e3e6909 100644
--- a/src/combinators/product.jl
+++ b/src/combinators/product.jl
@@ -1,40 +1,6 @@
-function as(d::PowerMeasure)
-    as(Array, as(d.parent), length.(d.axes)...)
-end
-
-function as(d::ProductMeasure{<:AbstractArray{<:Dirac}})
-    return asConst(testvalue.(marginals(d)))
-end
-
-function as(d::ProductMeasure{A}) where {A<:AbstractArray}
-    mar = marginals(d)
-    ts = map(as, mar)
-    if allequal(ts)
-        return as(Array, first(ts), size(ts))
-    else
-        error("Not yet implemented")
-    end
-end
 
 using MappedArrays
 
-function as(d::ProductMeasure{A}) where {A<:MappedArrays.ReadonlyMappedArray}
-    d1 = marginals(d).f(first(marginals(d).data))
-    as(Array, as(d1), size(marginals(d))...)
-end
-
-function as(d::ProductMeasure{T}) where {T<:Tuple}
-    as(map(as, d.marginals))
-end
-
-###############################################################################
-# I <: Base.Generator
-
-function as(d::ProductMeasure{<:Base.Generator})
-    d1 = marginals(d).f(first(marginals(d).iter))
-    as(Array, as(d1), size(marginals(d))...)
-end
-
 # function as(d::ProductMeasure{Returns{T},F,A}) where {T,F,A<:AbstractArray}
 #     as(Array, as(d.f.f.value), size(d.xs))
 # end
diff --git a/src/combinators/tweedie.jl b/src/combinators/tweedie.jl
deleted file mode 100644
index b18a3fb8..00000000
--- a/src/combinators/tweedie.jl
+++ /dev/null
@@ -1,159 +0,0 @@
-export Tweedie
-
-abstract type AbstractEDM <: AbstractTransitionKernel end
-
-"""
-From https://en.wikipedia.org/wiki/Tweedie_distribution:
-
-> The Tweedie distributions include a number of familiar distributions as well as
-some unusual ones, each being specified by the domain of the index parameter. We
-have the
-
-    extreme stable distribution, p < 0,
-    normal distribution, p = 0,
-    Poisson distribution, p = 1,
-    compound Poisson-gamma distribution, 1 < p < 2,
-    gamma distribution, p = 2,
-    positive stable distributions, 2 < p < 3,
-    Inverse Gaussian distribution, p = 3,
-    positive stable distributions, p > 3, and
-    extreme stable distributions, p = ∞.
-
-For 0 < p < 1 no Tweedie model exists. Note that all stable distributions mean
-actually generated by stable distributions. 
-"""
-struct Tweedie{S,B,D,P} <: AbstractEDM
-    support_contains::S
-    base::B
-    dim::D
-    p::P
-end
-
-struct TweedieMeasure{B,Θ,P,S,C} <: AbstractMeasure
-    fam::Tweedie{B,Θ,P}
-    θ::Θ
-    σ::S
-    cumulant::C
-end
-
-mean(d::TweedieMeasure) = tweedie_mean(d.fam.p, d.θ)
-
-var(d::TweedieMeasure) = d.σ^2 * mean(d)^d.fam.p
-
-###############################################################################
-# Tweedie cumulants
-
-@inline function tweedie_cumulant(p::P, θ) where {P}
-    if p == zero(P)
-        return 0.5 * θ^2
-    elseif p == one(P)
-        return exp(θ)
-    elseif p == 2
-        return -log(-θ)
-    else
-        α = (p - 2) / (p - 1)
-        coeff = (α - 1) / α
-        return coeff * (θ / (α - 1))^α
-    end
-end
-
-@inline function tweedie_cumulant(::StaticFloat64{0.0}, θ)
-    return 0.5 * θ^2
-end
-
-@inline function tweedie_cumulant(::StaticFloat64{1.0}, θ)
-    return exp(θ)
-end
-
-@inline function tweedie_cumulant(::StaticFloat64{2.0}, θ)
-    return -log(-θ)
-end
-
-@generated function tweedie_cumulant(::StaticFloat64{p}, θ) where {p}
-    α = (p - 2) / (p - 1)
-    coeff = (α - 1) / α
-
-    quote
-        $(Expr(:meta, :inline))
-        coeff * (θ / (α - 1))^α
-    end
-end
-
-@inline function (fam::Tweedie)(par)
-    base = fam.base(par.σ)
-    θ = fam.θ(par)
-    η = fam.η(θ)
-    t = fam.t
-    a = fam.a(θ)
-    TweedieMeasure(base, θ, p, σ)
-end
-
-###############################################################################
-# Tweedie mean function
-
-@inline function tweedie_mean(p::P, θ) where {P}
-    if p == zero(P)
-        return θ
-    elseif p == one(P)
-        return exp(θ)
-    elseif p == 2
-        return inv(log(-θ))
-    else
-        α_minus_1 = (p - 2) / (p - 1) - 1
-        return (θ / α_minus_1)^α_minus_1
-    end
-end
-
-@inline function tweedie_mean(::StaticFloat64{0.0}, θ)
-    return θ
-end
-
-@inline function tweedie_mean(::StaticFloat64{1.0}, θ)
-    return exp(θ)
-end
-
-@inline function tweedie_mean(::StaticFloat64{2.0}, θ)
-    return inv(log(-θ))
-end
-
-@generated function tweedie_mean(::StaticFloat64{p}, θ) where {p}
-    α_minus_1 = (p - 2) / (p - 1) - 1
-
-    quote
-        $(Expr(:meta, :inline))
-        (θ / α_minus_1)^α_minus_1
-    end
-end
-
-basemeasure(d::TweedieMeasure) = d.base
-
-function logdensity_def(d::TweedieMeasure, x)
-    mydot(x, d.θ) - d.cumulant
-end
-
-function MeasureBase.powermeasure(fam::Tweedie, dims::NTuple{N,I}) where {N,I}
-    base(σ) = fam.base(σ)^dims
-    a = AffineTransform((σ = prod(dims),)) ∘ fam.a
-    Tweedie(fam.base^dims, fam.θ, fam.η, t, a)
-end
-
-struct TweedieLikelihood{C,Θ,H,T,A} <: AbstractLikelihood
-    c::C
-    θ::Θ
-    η::H
-    t::T
-    a::A
-end
-
-function likelihoodof(fam::Tweedie, x)
-    c = logdensityof(fam.base, x)
-    t = fam.t(x)
-    TweedieLikelihood(c, fam.θ, fam.η, t, fam.a)
-end
-
-@inline function logdensity_def(ℓ::TweedieLikelihood, par)
-    θ = ℓ.θ(par)
-    mydot(θ, ℓ.t) - ℓ.a(θ) + ℓ.c
-end
-
-basemeasure(fam::Tweedie) = fam.base
diff --git a/src/distributions.jl b/src/distributions.jl
index ce79f520..43e85f58 100644
--- a/src/distributions.jl
+++ b/src/distributions.jl
@@ -39,4 +39,3 @@ function as(d::Dists.Product, _data::NamedTuple = NamedTuple())
     as(Vector, as(v[1]), n)
 end
 
-as(m::DistributionMeasures.DistributionMeasure) = as(m.d)
diff --git a/src/parameterized.jl b/src/parameterized.jl
index 08a4bf10..9ae51e22 100644
--- a/src/parameterized.jl
+++ b/src/parameterized.jl
@@ -75,3 +75,127 @@ function asparams(μ::M, nt::NamedTuple = NamedTuple()) where {M<:ParameterizedM
 end
 
 as(::Half) = asℝ₊
+
+asparams(::AffinePushfwd, ::StaticSymbol{:μ}) = asℝ
+asparams(::AffinePushfwd, ::StaticSymbol{:σ}) = asℝ₊
+asparams(::Type{A}, ::StaticSymbol{:μ}) where {A<:AffinePushfwd} = asℝ
+asparams(::Type{A}, ::StaticSymbol{:σ}) where {A<:AffinePushfwd} = asℝ₊
+
+function asparams(d::AffinePushfwd{N,M,T}, ::StaticSymbol{:μ}) where {N,M,T<:AbstractArray}
+    as(Array, asℝ, size(d.μ))
+end
+
+function asparams(d::AffinePushfwd{N,M,T}, ::StaticSymbol{:σ}) where {N,M,T<:AbstractArray}
+    as(Array, asℝ, size(d.σ))
+end
+
+asparams(::Type{<:Bernoulli}, ::StaticSymbol{:p}) = as𝕀
+asparams(::Type{<:Bernoulli}, ::StaticSymbol{:logitp}) = asℝ
+
+asparams(::Type{<:Beta}, ::StaticSymbol{:α}) = asℝ₊
+asparams(::Type{<:Beta}, ::StaticSymbol{:β}) = asℝ₊
+
+asparams(::Type{<:BetaBinomial}, ::StaticSymbol{:α}) = asℝ₊
+asparams(::Type{<:BetaBinomial}, ::StaticSymbol{:β}) = asℝ₊
+
+asparams(::Type{<:Binomial}, ::StaticSymbol{:p}) = as𝕀
+asparams(::Type{<:Binomial}, ::StaticSymbol{:logitp}) = asℝ
+asparams(::Type{<:Binomial}, ::StaticSymbol{:probitp}) = asℝ
+
+asparams(::Type{<:Exponential}, ::StaticSymbol{:β}) = asℝ₊
+asparams(::Type{<:Exponential}, ::StaticSymbol{:logβ}) = asℝ
+asparams(::Type{<:Exponential}, ::StaticSymbol{:λ}) = asℝ₊
+asparams(::Type{<:Exponential}, ::StaticSymbol{:logλ}) = asℝ
+
+asparams(::Type{<:LKJCholesky}, ::StaticSymbol{:η}) = asℝ₊
+asparams(::Type{<:LKJCholesky}, ::StaticSymbol{:logη}) = asℝ
+
+asparams(::Type{<:NegativeBinomial}, ::StaticSymbol{:p}) = as𝕀
+asparams(::Type{<:NegativeBinomial}, ::StaticSymbol{:logitp}) = asℝ
+asparams(::Type{<:NegativeBinomial}, ::StaticSymbol{:r}) = asℝ₊
+asparams(::Type{<:NegativeBinomial}, ::StaticSymbol{:λ}) = asℝ₊
+asparams(::Type{<:NegativeBinomial}, ::StaticSymbol{:logλ}) = asℝ
+
+asparams(::Type{<:Normal}, ::StaticSymbol{:μ}) = asℝ
+asparams(::Type{<:Normal}, ::StaticSymbol{:σ}) = asℝ₊
+asparams(::Type{<:Normal}, ::StaticSymbol{:σ²}) = asℝ₊
+asparams(::Type{<:Normal}, ::StaticSymbol{:τ}) = asℝ₊
+asparams(::Type{<:Normal}, ::StaticSymbol{:logτ}) = asℝ
+
+asparams(::Type{<:Poisson}, ::StaticSymbol{:λ}) = asℝ₊
+asparams(::Type{<:Poisson}, ::StaticSymbol{:logλ}) = asℝ
+
+asparams(::Type{<:SnedecorF}, ::StaticSymbol{:ν1}) = asℝ₊
+asparams(::Type{<:SnedecorF}, ::StaticSymbol{:ν2}) = asℝ₊
+
+asparams(::Type{<:StudentT}, ::StaticSymbol{:ν}) = asℝ₊
+
+function as(d::PowerMeasure)
+    as(Array, as(d.parent), length.(d.axes)...)
+end
+
+function as(d::ProductMeasure{<:AbstractArray{<:Dirac}})
+    return asConst(testvalue.(marginals(d)))
+end
+
+function as(d::ProductMeasure{A}) where {A<:AbstractArray}
+    mar = marginals(d)
+    ts = map(as, mar)
+    if allequal(ts)
+        return as(Array, first(ts), size(ts))
+    else
+        error("Not yet implemented")
+    end
+end
+
+function as(d::ProductMeasure{A}) where {A<:MappedArrays.ReadonlyMappedArray}
+    d1 = marginals(d).f(first(marginals(d).data))
+    as(Array, as(d1), size(marginals(d))...)
+end
+
+function as(d::ProductMeasure{T}) where {T<:Tuple}
+    as(map(as, d.marginals))
+end
+
+function as(d::ProductMeasure{<:Base.Generator})
+    d1 = marginals(d).f(first(marginals(d).iter))
+    as(Array, as(d1), size(marginals(d))...)
+end
+
+as(::Beta) = as𝕀
+as(::Cauchy) = asℝ
+as(d::Dirichlet{(:α,)}) = TV.UnitSimplex(length(d.α))
+as(::Exponential) = asℝ₊
+as(::Gamma) = asℝ₊
+as(::Gumbel) = asℝ
+as(::InverseGaussian) = asℝ₊
+as(::Laplace) = asℝ
+as(d::MvNormal{(:μ,)}) = as(Array, length(d.μ))
+
+as(d::MvNormal{(:Σ,),Tuple{C}}) where {C<:Cholesky} = as(Array, size(d.Σ, 1))
+as(d::MvNormal{(:Λ,),Tuple{C}}) where {C<:Cholesky} = as(Array, size(d.Λ, 1))
+as(d::MvNormal{(:μ, :Σ),<:Tuple{T,C}}) where {T,C<:Cholesky} = as(Array, size(d.Σ, 1))
+as(d::MvNormal{(:μ, :Λ),<:Tuple{T,C}}) where {T,C<:Cholesky} = as(Array, size(d.Λ, 1))
+
+function as(d::MvNormal{(:σ,),Tuple{M}}) where {M<:Triangular}
+    σ = d.σ
+    if @inbounds all(i -> σ[i] ≠ 0, diagind(σ))
+        return as(Array, size(σ, 1))
+    else
+        @error "Not implemented yet"
+    end
+end
+
+function as(d::MvNormal{(:λ,),Tuple{M}}) where {M<:Triangular}
+    λ = d.λ
+    if @inbounds all(i -> λ[i] > 0, diagind(λ))
+        return as(Array, size(λ, 1))
+    else
+        @error "Not implemented yet"
+    end
+end
+
+as(::Normal) = asℝ
+as(::SnedecorF) = asℝ₊
+as(::StudentT) = asℝ
+as(::Uniform{()}) = as𝕀
diff --git a/src/parameterized/bernoulli.jl b/src/parameterized/bernoulli.jl
index f49485ff..b1d0f484 100644
--- a/src/parameterized/bernoulli.jl
+++ b/src/parameterized/bernoulli.jl
@@ -13,7 +13,7 @@ import Base
 
 Bernoulli(p) = Bernoulli((p = p,))
 
-basemeasure(::Bernoulli) = CountingMeasure()
+basemeasure(::Bernoulli) = CountingBase()
 
 testvalue(::Bernoulli) = true
 
@@ -37,6 +37,9 @@ function density_def(d::Bernoulli{(:p,)}, y)
     return 2 * p * y - p - y + 1
 end
 
+densityof(d::Bernoulli, y) = density_def(d, y)
+unsafe_densityof(d::Bernoulli, y) = density_def(d, y)
+
 @inline function logdensity_def(d::Bernoulli{(:logitp,)}, y)
     x = d.logitp
     return y * x - log1pexp(x)
@@ -57,8 +60,17 @@ function Base.rand(rng::AbstractRNG, T::Type, d::Bernoulli{(:logitp,)})
     rand(rng, T) < logistic(d.logitp)
 end
 
-asparams(::Type{<:Bernoulli}, ::StaticSymbol{:p}) = as𝕀
-asparams(::Type{<:Bernoulli}, ::StaticSymbol{:logitp}) = asℝ
+function smf(b::B, x::X) where {B<:Bernoulli,X}
+    T = Core.Compiler.return_type(densityof, Tuple{B,X})
+    x < zero(x) && return zero(T)
+    x ≥ one(x) && return one(T)
+    return densityof(b, zero(x))
+end
+
+function invsmf(b::Bernoulli, p)
+    p0 = densityof(b, 0)
+    p ≤ p0 ? 0 : 1
+end
 
 proxy(d::Bernoulli{(:p,)}) = Dists.Bernoulli(d.p)
 proxy(d::Bernoulli{(:logitp,)}) = Dists.Bernoulli(logistic(d.logitp))
diff --git a/src/parameterized/beta.jl b/src/parameterized/beta.jl
index f6520652..744f79f4 100644
--- a/src/parameterized/beta.jl
+++ b/src/parameterized/beta.jl
@@ -13,23 +13,26 @@ export Beta
     beta => β
 ]
 
-as(::Beta) = as𝕀
-
 @inline function logdensity_def(d::Beta{(:α, :β),Tuple{A,B}}, x::X) where {A,B,X}
     return xlogy(d.α - 1, x) + xlog1py(d.β - 1, -x)
 end
 
 @inline function basemeasure(d::Beta{(:α, :β)})
     ℓ = -logbeta(d.α, d.β)
-    weightedmeasure(ℓ, LebesgueMeasure())
+    weightedmeasure(ℓ, LebesgueBase())
 end
 
 Base.rand(rng::AbstractRNG, T::Type, μ::Beta) = rand(rng, Dists.Beta(μ.α, μ.β))
 
-proxy(d::Beta{(:α, :β)}) = Dists.Beta(d.α, d.β)
-
-asparams(::Type{<:Beta}, ::StaticSymbol{:α}) = asℝ₊
-asparams(::Type{<:Beta}, ::StaticSymbol{:β}) = asℝ₊
-
 insupport(::Beta, x) = in𝕀(x)
 insupport(::Beta) = in𝕀
+
+function smf(d::Beta{(:α, :β)}, x::Real)
+    StatsFuns.betacdf(d.α, d.β, x)
+end
+
+function invsmf(d::Beta{(:α, :β)}, p)
+    StatsFuns.betainvcdf(d.α, d.β, p)
+end
+
+proxy(d::Beta{(:α, :β)}) = Dists.Beta(d.α, d.β)
diff --git a/src/parameterized/betabinomial.jl b/src/parameterized/betabinomial.jl
index 00d48947..16c10d9b 100644
--- a/src/parameterized/betabinomial.jl
+++ b/src/parameterized/betabinomial.jl
@@ -6,17 +6,13 @@ using SpecialFunctions
 
 @parameterized BetaBinomial(n, α, β)
 
-basemeasure(d::BetaBinomial) = CountingMeasure()
+basemeasure(d::BetaBinomial) = CountingBase()
 
 testvalue(::BetaBinomial) = 0
 
 @kwstruct BetaBinomial(n, α, β)
 
-function Base.rand(
-    rng::AbstractRNG,
-    ::Type{T},
-    d::BetaBinomial{(:n, :α, :β)},
-) where {T}
+function Base.rand(rng::AbstractRNG, ::Type{T}, d::BetaBinomial{(:n, :α, :β)}) where {T}
     k = rand(rng, T, Beta(d.α, d.β))
     return rand(rng, T, Binomial(d.n, k))
 end
@@ -33,5 +29,4 @@ end
     return logbinom + lognum - logdenom
 end
 
-asparams(::Type{<:BetaBinomial}, ::StaticSymbol{:α}) = asℝ₊
-asparams(::Type{<:BetaBinomial}, ::StaticSymbol{:β}) = asℝ₊
+proxy(d::BetaBinomial{(:n, :α, :β)}) = Dists.BetaBinomial(d.n, d.α, d.β)
diff --git a/src/parameterized/binomial.jl b/src/parameterized/binomial.jl
index 1825cb69..51f89921 100644
--- a/src/parameterized/binomial.jl
+++ b/src/parameterized/binomial.jl
@@ -4,12 +4,11 @@ export Binomial
 import Base
 using SpecialFunctions
 
-probit(p) = sqrt2 * erfinv(2p - 1)
-Φ(z) = (1 + erf(invsqrt2 * z)) / 2
+probit(p) = Φinv(p)
 
 @parameterized Binomial(n, p)
 
-basemeasure(d::Binomial) = CountingMeasure()
+basemeasure(d::Binomial) = CountingBase()
 
 testvalue(::Binomial) = 0
 
@@ -84,13 +83,19 @@ function Base.rand(
 end
 
 proxy(d::Binomial{(:n, :p),Tuple{I,A}}) where {I<:Integer,A} = Dists.Binomial(d.n, d.p)
+
 function proxy(d::Binomial{(:n, :logitp),Tuple{I,A}}) where {I<:Integer,A}
     Dists.Binomial(d.n, logistic(d.logitp))
 end
+
 function proxy(d::Binomial{(:n, :probitp),Tuple{I,A}}) where {I<:Integer,A}
     Dists.Binomial(d.n, Φ(d.probitp))
 end
 
-asparams(::Type{<:Binomial}, ::StaticSymbol{:p}) = as𝕀
-asparams(::Type{<:Binomial}, ::StaticSymbol{:logitp}) = asℝ
-asparams(::Type{<:Binomial}, ::StaticSymbol{:probitp}) = asℝ
+function smf(μ::Binomial{(:n, :p)}, k)
+    α = max(0, floor(k) + 1)
+    β = max(0, μ.n - floor(k))
+
+    beta_smf = smf(Beta(α, β), μ.p)
+    one(beta_smf) - beta_smf
+end
diff --git a/src/parameterized/cauchy.jl b/src/parameterized/cauchy.jl
index f37c23f9..d49fa48a 100644
--- a/src/parameterized/cauchy.jl
+++ b/src/parameterized/cauchy.jl
@@ -14,7 +14,7 @@ export Cauchy, HalfCauchy
 
 logdensity_def(d::Cauchy, x) = logdensity_def(proxy(d), x)
 
-basemeasure(d::Cauchy) = WeightedMeasure(static(-logπ), LebesgueMeasure())
+basemeasure(d::Cauchy) = WeightedMeasure(static(-logπ), LebesgueBase())
 
 # @affinepars Cauchy
 
@@ -28,6 +28,8 @@ end
 
 Base.rand(rng::AbstractRNG, T::Type, μ::Cauchy{()}) = randn(rng, T) / randn(rng, T)
 
+Base.rand(::FixedRNG, ::Type{T}, ::Cauchy{()}) where {T} = zero(T)
+
 for N in AFFINEPARS
     @eval begin
         proxy(d::Cauchy{$N}) = affine(params(d), Cauchy())
@@ -40,12 +42,20 @@ end
 
 ≪(::Cauchy, ::Lebesgue{X}) where {X<:Real} = true
 
-as(::Cauchy) = asℝ
-
 @half Cauchy
 
 HalfCauchy(σ) = HalfCauchy(σ = σ)
 
-proxy(d::Cauchy{()}) = Dists.Cauchy()
+Base.rand(::FixedRNG, ::Type{T}, μ::Half{<:Cauchy}) where {T} = one(T)
 
 insupport(::Cauchy, x) = true
+
+function smf(::Cauchy{()}, x)
+    invπ * atan(x) + 0.5
+end
+
+function invsmf(::Cauchy{()}, p)
+    tan(π * (p - 0.5))
+end
+
+proxy(d::Cauchy{()}) = Dists.Cauchy()
diff --git a/src/parameterized/dirichlet.jl b/src/parameterized/dirichlet.jl
index bd2940b0..b2debe52 100644
--- a/src/parameterized/dirichlet.jl
+++ b/src/parameterized/dirichlet.jl
@@ -7,12 +7,9 @@ using FillArrays
 
 @parameterized Dirichlet(α)
 
-as(d::Dirichlet{(:α,)}) = TV.UnitSimplex(length(d.α))
-
 @inline function basemeasure(μ::Dirichlet{(:α,)})
     α = μ.α
-    t = as(μ)
-    d = TV.dimension(t)
+    # TODO: We can make this traverse α only once. Is that faster?
     logw = loggamma(sum(α)) - sum(loggamma, α)
     return WeightedMeasure(logw, Lebesgue(Simplex()))
 end
@@ -22,23 +19,25 @@ end
 Dirichlet(k::Integer, α) = Dirichlet(Fill(α, k))
 
 @inline function logdensity_def(d::Dirichlet{(:α,)}, x)
-    α = d.α
-    s = 0.0
-    for j in eachindex(x)
-        s += xlogy(α[j] - 1, x[j])
+    sum(zip(d.α, x)) do (αj, xj)
+        xlogy(αj - 1, xj)
     end
-    return s
+
+    # `mapreduce` is slow for this case
+    # mapreduce(+, d.α, x) do αj, xj
+    #     xlogy(αj - 1, xj)
+    # end
 end
 
 Base.rand(rng::AbstractRNG, T::Type, μ::Dirichlet) = rand(rng, Dists.Dirichlet(μ.α))
 
-proxy(d::Dirichlet{(:α,)}) = Dists.Dirichlet(d.α)
-
-function testvalue(d::Dirichlet{(:α,)})
+function testvalue(::Type{T}, d::Dirichlet{(:α,)}) where {T}
     n = length(d.α)
-    Fill(1 / n, n)
+    Fill(inv(convert(T, n)), n)
 end
 
 @inline function insupport(d::Dirichlet{(:α,)}, x)
     length(x) == length(d.α) && sum(x) ≈ 1.0
 end
+
+as(μ::Dirichlet) = TV.UnitSimplex(length(μ.α))
diff --git a/src/parameterized/exponential.jl b/src/parameterized/exponential.jl
index cfd40bee..b2ac6c7a 100644
--- a/src/parameterized/exponential.jl
+++ b/src/parameterized/exponential.jl
@@ -6,7 +6,7 @@ export Exponential
 @parameterized Exponential(β)
 
 insupport(::Exponential, x) = x ≥ 0
-basemeasure(::Exponential) = LebesgueMeasure()
+basemeasure(::Exponential) = LebesgueBase()
 
 @kwstruct Exponential()
 
@@ -16,8 +16,6 @@ end
 
 Base.rand(rng::AbstractRNG, T::Type, μ::Exponential{()}) = randexp(rng, T)
 
-as(::Exponential) = asℝ₊
-
 ##########################
 # Scale β
 
@@ -34,8 +32,6 @@ end
 
 proxy(d::Exponential{(:β,)}) = Dists.Exponential(d.β)
 
-asparams(::Type{<:Exponential}, ::StaticSymbol{:β}) = asℝ₊
-
 ##########################
 # Log-Scale logβ
 
@@ -52,8 +48,6 @@ end
 
 proxy(d::Exponential{(:logβ,)}) = Dists.Exponential(exp(d.logβ))
 
-asparams(::Type{<:Exponential}, ::StaticSymbol{:logβ}) = asℝ
-
 ##########################
 # Rate λ
 
@@ -70,8 +64,6 @@ end
 
 proxy(d::Exponential{(:λ,)}) = Dists.Exponential(1 / d.λ)
 
-asparams(::Type{<:Exponential}, ::StaticSymbol{:λ}) = asℝ₊
-
 ##########################
 # Log-Rate logλ
 
@@ -88,4 +80,6 @@ end
 
 proxy(d::Exponential{(:logλ,)}) = Dists.Exponential(exp(-d.logλ))
 
-asparams(::Type{<:Exponential}, ::StaticSymbol{:logλ}) = asℝ
+smf(::Exponential{()}, x) = smf(StdExponential(), x)
+
+invsmf(::Exponential{()}, p) = invsmf(StdExponential(), p)
diff --git a/src/parameterized/gamma.jl b/src/parameterized/gamma.jl
index 2e40375f..90c99af7 100644
--- a/src/parameterized/gamma.jl
+++ b/src/parameterized/gamma.jl
@@ -19,7 +19,7 @@ end
 
 function basemeasure(d::Gamma{(:k,)})
     ℓ = -loggamma(d.k)
-    weightedmeasure(ℓ, LebesgueMeasure())
+    weightedmeasure(ℓ, LebesgueBase())
 end
 
 @kwstruct Gamma(k, σ)
@@ -28,7 +28,7 @@ Gamma(k, σ) = Gamma((k = k, σ = σ))
 
 @useproxy Gamma{(:k, :σ)}
 function proxy(d::Gamma{(:k, :σ)})
-    affine(NamedTuple{(:σ,)}(d.σ), Gamma((k = d.k,)))
+    affine((σ = d.σ,), Gamma((k = d.k,)))
 end
 
 @kwstruct Gamma(k, λ)
@@ -42,8 +42,6 @@ Base.rand(rng::AbstractRNG, T::Type, μ::Gamma{()}) = rand(rng, T, Exponential()
 
 Base.rand(rng::AbstractRNG, T::Type, μ::Gamma{(:k,)}) = rand(rng, Dists.Gamma(μ.k))
 
-as(::Gamma) = asℝ₊
-
 insupport(::Gamma, x) = x > 0
 
 @kwstruct Gamma(μ, ϕ)
@@ -57,11 +55,11 @@ function basemeasure(d::Gamma{(:μ, :ϕ)})
     ϕ = d.ϕ
     ϕinv = inv(ϕ)
     ℓ = -ϕinv * log(ϕ) - first(logabsgamma(ϕinv))
-    weightedmeasure(ℓ, LebesgueMeasure())
+    weightedmeasure(ℓ, LebesgueBase())
 end
 
 function basemeasure(d::Gamma{(:μ, :ϕ),Tuple{M,StaticFloat64{ϕ}}}) where {M,ϕ}
     ϕinv = inv(ϕ)
     ℓ = static(-ϕinv * log(ϕ) - first(logabsgamma(ϕinv)))
-    weightedmeasure(ℓ, LebesgueMeasure())
+    weightedmeasure(ℓ, LebesgueBase())
 end
diff --git a/src/parameterized/gumbel.jl b/src/parameterized/gumbel.jl
index c23f4b5a..984f9b28 100644
--- a/src/parameterized/gumbel.jl
+++ b/src/parameterized/gumbel.jl
@@ -4,7 +4,7 @@ export Gumbel
 
 @parameterized Gumbel()
 
-basemeasure(::Gumbel{()}) = LebesgueMeasure()
+basemeasure(::Gumbel{()}) = LebesgueBase()
 
 @kwstruct Gumbel()
 
@@ -28,15 +28,13 @@ end
 
 import Base
 
-function Base.rand(rng::AbstractRNG, d::Gumbel{()})
-    u = rand(rng)
+function Base.rand(rng::AbstractRNG, ::Type{T}, d::Gumbel{()}) where {T}
+    u = rand(rng, T)
     -log(-log(u))
 end
 
-as(::Gumbel) = asℝ
-
 ≪(::Gumbel, ::Lebesgue{X}) where {X<:Real} = true
 
-proxy(::Gumbel{()}) = Dists.Gumbel()
-
 insupport(::Gumbel, x) = true
+
+proxy(::Gumbel{()}) = Dists.Gumbel()
diff --git a/src/parameterized/inverse-gaussian.jl b/src/parameterized/inverse-gaussian.jl
index ee34e455..5c10ac4e 100644
--- a/src/parameterized/inverse-gaussian.jl
+++ b/src/parameterized/inverse-gaussian.jl
@@ -35,8 +35,6 @@ export InverseGaussian
 
 Base.rand(rng::AbstractRNG, T::Type, d::InverseGaussian) = rand(rng, proxy(d))
 
-as(::InverseGaussian) = asℝ₊
-
 insupport(::InverseGaussian, x) = x > 0
 
 # GLM parameterization
@@ -49,5 +47,5 @@ end
 
 function basemeasure(d::InverseGaussian{(:μ, :ϕ)})
     ℓ = static(-0.5) * (static(float(log2π)) + log(d.ϕ))
-    weightedmeasure(ℓ, LebesgueMeasure())
+    weightedmeasure(ℓ, LebesgueBase())
 end
diff --git a/src/parameterized/laplace.jl b/src/parameterized/laplace.jl
index d5f9a112..9a691fe8 100644
--- a/src/parameterized/laplace.jl
+++ b/src/parameterized/laplace.jl
@@ -25,17 +25,18 @@ insupport(::Laplace, x) = true
     return -abs(x)
 end
 
-logdensity_def(d::Laplace, x) = logdensity_def(proxy(d), x)
+Laplace(μ, σ) = Laplace((μ = μ, σ = σ))
 
-basemeasure(::Laplace{()}) = WeightedMeasure(static(-logtwo), LebesgueMeasure())
+basemeasure(::Laplace{()}) = WeightedMeasure(static(-logtwo), LebesgueBase())
 
 # @affinepars Laplace
 
 function Base.rand(rng::AbstractRNG, ::Type{T}, μ::Laplace{()}) where {T}
-    rand(rng, Dists.Laplace())
+    sign = rand(rng, Bool)
+    absx = randexp(rng, T)
+    sign == true ? absx : -absx
 end
+
 Base.rand(rng::AbstractRNG, ::Type{T}, μ::Laplace) where {T} = Base.rand(rng, T, proxy(μ))
 
 ≪(::Laplace, ::Lebesgue{X}) where {X<:Real} = true
-
-as(::Laplace) = asℝ
diff --git a/src/parameterized/lkj-cholesky.jl b/src/parameterized/lkj-cholesky.jl
index 72e8d2e5..11c3eb3a 100644
--- a/src/parameterized/lkj-cholesky.jl
+++ b/src/parameterized/lkj-cholesky.jl
@@ -38,9 +38,6 @@ https://github.com/tpapp/AltDistributions.jl
 
 LKJCholesky(k::Integer) = LKJCholesky(k, 1.0)
 
-asparams(::Type{<:LKJCholesky}, ::StaticSymbol{:η}) = asℝ₊
-asparams(::Type{<:LKJCholesky}, ::StaticSymbol{:logη}) = asℝ
-
 # Modified from
 # https://github.com/tpapp/AltDistributions.jl
 
@@ -81,14 +78,14 @@ as(d::LKJCholesky) = CorrCholesky(d.k)
 
 @inline function basemeasure(d::LKJCholesky{(:k, :η)})
     t = as(d)
-    base = Pushforward(t, LebesgueMeasure()^TV.dimension(t), False())
+    base = Pushforward(t, LebesgueBase()^TV.dimension(t), False())
     WeightedMeasure(Dists.lkj_logc0(d.k, d.η), base)
 end
 
 @inline function basemeasure(d::LKJCholesky{(:k, :logη)})
     t = as(d)
     η = exp(d.logη)
-    base = Pushforward(t, LebesgueMeasure()^TV.dimension(t), False())
+    base = Pushforward(t, LebesgueBase()^TV.dimension(t), False())
     WeightedMeasure(Dists.lkj_logc0(d.k, η), base)
 end
 
diff --git a/src/parameterized/multinomial.jl b/src/parameterized/multinomial.jl
index 781de7d7..4119f71d 100644
--- a/src/parameterized/multinomial.jl
+++ b/src/parameterized/multinomial.jl
@@ -4,7 +4,7 @@ export Multinomial
 
 @parameterized Multinomial(n, p)
 
-basemeasure(d::Multinomial) = CountingMeasure()
+basemeasure(d::Multinomial) = CountingBase()
 
 @inline function insupport(d::Multinomial{(:n, :p)}, x)
     length(x) == length(d.p) || return false
diff --git a/src/parameterized/mvnormal.jl b/src/parameterized/mvnormal.jl
index 31d7ebc9..57cdae47 100644
--- a/src/parameterized/mvnormal.jl
+++ b/src/parameterized/mvnormal.jl
@@ -20,31 +20,6 @@ export MvNormal
 @kwstruct MvNormal(μ, Σ)
 @kwstruct MvNormal(μ, Λ)
 
-as(d::MvNormal{(:μ,)}) = as(Array, length(d.μ))
-
-as(d::MvNormal{(:Σ,),Tuple{C}}) where {C<:Cholesky} = as(Array, size(d.Σ, 1))
-as(d::MvNormal{(:Λ,),Tuple{C}}) where {C<:Cholesky} = as(Array, size(d.Λ, 1))
-as(d::MvNormal{(:μ, :Σ),<:Tuple{T,C}}) where {T,C<:Cholesky} = as(Array, size(d.Σ, 1))
-as(d::MvNormal{(:μ, :Λ),<:Tuple{T,C}}) where {T,C<:Cholesky} = as(Array, size(d.Λ, 1))
-
-function as(d::MvNormal{(:σ,),Tuple{M}}) where {M<:Triangular}
-    σ = d.σ
-    if @inbounds all(i -> σ[i] ≠ 0, diagind(σ))
-        return as(Array, size(σ, 1))
-    else
-        @error "Not implemented yet"
-    end
-end
-
-function as(d::MvNormal{(:λ,),Tuple{M}}) where {M<:Triangular}
-    λ = d.λ
-    if @inbounds all(i -> λ[i] > 0, diagind(λ))
-        return as(Array, size(λ, 1))
-    else
-        @error "Not implemented yet"
-    end
-end
-
 for N in setdiff(AFFINEPARS, [(:μ,)])
     @eval begin
         function as(d::MvNormal{$N})
diff --git a/src/parameterized/negativebinomial.jl b/src/parameterized/negativebinomial.jl
index dcca6741..227a6435 100644
--- a/src/parameterized/negativebinomial.jl
+++ b/src/parameterized/negativebinomial.jl
@@ -7,7 +7,7 @@ import Base
 
 insupport(::NegativeBinomial, x) = isinteger(x) && x ≥ 0
 
-basemeasure(::NegativeBinomial) = CountingMeasure()
+basemeasure(::NegativeBinomial) = CountingBase()
 
 testvalue(::NegativeBinomial) = 0
 
@@ -76,9 +76,3 @@ function proxy(d::NegativeBinomial{(:r, :logλ)})
 end
 
 ###############################################################################
-
-asparams(::Type{<:NegativeBinomial}, ::StaticSymbol{:p}) = as𝕀
-asparams(::Type{<:NegativeBinomial}, ::StaticSymbol{:logitp}) = asℝ
-asparams(::Type{<:NegativeBinomial}, ::StaticSymbol{:r}) = asℝ₊
-asparams(::Type{<:NegativeBinomial}, ::StaticSymbol{:λ}) = asℝ₊
-asparams(::Type{<:NegativeBinomial}, ::StaticSymbol{:logλ}) = asℝ
diff --git a/src/parameterized/normal.jl b/src/parameterized/normal.jl
index aa4d37a1..fd82dd66 100644
--- a/src/parameterized/normal.jl
+++ b/src/parameterized/normal.jl
@@ -25,6 +25,8 @@ export Normal, HalfNormal
 #
 @parameterized Normal()
 
+massof(::Normal) = static(1.0)
+
 for N in AFFINEPARS
     @eval begin
         proxy(d::Normal{$N,T}) where {T} = affine(params(d), Normal())
@@ -32,12 +34,12 @@ for N in AFFINEPARS
     end
 end
 
-insupport(d::Normal, x) = true
+insupport(d::Normal, x) = True()
 
-insupport(d::Normal) = Returns(true)
+insupport(d::Normal) = Returns(True())
 
-@inline logdensity_def(d::Normal{()}, x) = -x^2 / 2
-@inline basemeasure(::Normal{()}) = WeightedMeasure(static(-0.5 * log2π), LebesgueMeasure())
+@inline logdensity_def(d::Normal{()}, x) = -muladd(x, x, log2π) / 2;
+@inline basemeasure(::Normal{()}) = LebesgueBase()
 
 @kwstruct Normal(μ)
 @kwstruct Normal(σ)
@@ -51,8 +53,6 @@ Normal(μ::M, σ::S) where {M,S} = Normal((μ = μ, σ = σ))::Normal{(:μ, :σ)
 
 # Normal(nt::NamedTuple{N,Tuple{Vararg{AbstractArray}}}) where {N} = MvNormal(nt)
 
-as(::Normal) = asℝ
-
 # `@kwalias` defines some alias names, giving users flexibility in the names
 # they use. For example, σ² is standard notation for the variance parameter, but
 # it's a lot to type. Some users might prefer to just use `var` and have us do
@@ -81,10 +81,6 @@ as(::Normal) = asℝ
 # (μ = -0.4548087051528626, σ² = 11.920775478312793)
 #
 # And of course, you can apply `Normal` to any one of the above.
-#
-asparams(::Type{<:Normal}, ::StaticSymbol{:σ²}) = asℝ₊
-asparams(::Type{<:Normal}, ::StaticSymbol{:τ}) = asℝ₊
-asparams(::Type{<:Normal}, ::StaticSymbol{:logτ}) = asℝ
 
 # Rather than try to reimplement everything in Distributions, measures can have
 # a `proxy` method. This just delegates some methods to the corresponding
@@ -147,7 +143,7 @@ end
 
 @inline function basemeasure(d::Normal{(:σ²,)})
     ℓ = static(-0.5) * (static(float(log2π)) + log(d.σ²))
-    weightedmeasure(ℓ, LebesgueMeasure())
+    weightedmeasure(ℓ, LebesgueBase())
 end
 
 proxy(d::Normal{(:μ, :σ²)}) = affine((μ = d.μ,), Normal((σ² = d.σ²,)))
@@ -163,7 +159,7 @@ end
 
 @inline function basemeasure(d::Normal{(:τ,)})
     ℓ = static(-0.5) * (static(float(log2π)) - log(d.τ))
-    weightedmeasure(ℓ, LebesgueMeasure())
+    weightedmeasure(ℓ, LebesgueBase())
 end
 
 proxy(d::Normal{(:μ, :τ)}) = affine((μ = d.μ,), Normal((τ = d.τ,)))
@@ -175,7 +171,12 @@ proxy(d::Normal{(:μ, :τ)}) = affine((μ = d.μ,), Normal((τ = d.τ,)))
 @inline function logdensity_def(d::Normal{(:μ, :logσ)}, x)
     μ = d.μ
     logσ = d.logσ
-    -logσ - 0.5(exp(-2logσ) * ((x - μ)^2))
+    -0.5(exp(-2 * logσ) * ((x - μ)^2))
+end
+
+function basemeasure(d::Normal{(:μ, :logσ)})
+    ℓ = static(-0.5) * (static(float(log2π)) + static(2.0) * d.logσ)
+    weightedmeasure(ℓ, LebesgueBase())
 end
 
 function logdensity_def(p::Normal, q::Normal, x)
@@ -184,18 +185,26 @@ function logdensity_def(p::Normal, q::Normal, x)
     μq = mean(q)
     σq = std(q)
 
-    # Result is (((x - μq) / σq)^2 - ((x - μp) / σp)^2) / 2 
+    # Result is (((x - μq) / σq)^2 - ((x - μp) / σp)^2 + log(abs(σq / σp))) / 2 
 
     # We'll write the difference of squares as sqdiff, then divide that by 2 at
     # the end
 
-    sqdiff = if σp == σq
-        (2x - μq - μp) * (μq - μp) / σp^2
+    if σp == σq
+        return (2x - μq - μp) * (μp - μq) / (2 * σp^2)
     else
         zp = (x - μp) / σp
         zq = (x - μq) / σq
-        (zq + zp) * (zq - zp)
+        return ((zq + zp) * (zq - zp)) / 2 + log(abs(σq / σp))
     end
-
-    return sqdiff / 2
 end
+
+MeasureBase.transport_origin(::Normal) = StdNormal()
+
+MeasureBase.to_origin(::Normal{()}, y) = y
+MeasureBase.from_origin(::Normal{()}, x) = x
+
+MeasureBase.smf(::Normal{()}, x) = Φ(x)
+MeasureBase.invsmf(::Normal{()}, p) = Φinv(p)
+
+@smfAD Normal{()}
diff --git a/src/parameterized/poisson.jl b/src/parameterized/poisson.jl
index 7135fd70..8bae5f45 100644
--- a/src/parameterized/poisson.jl
+++ b/src/parameterized/poisson.jl
@@ -28,9 +28,6 @@ end
     return y * d.logλ - exp(d.logλ) - loggamma(1 + y)
 end
 
-asparams(::Type{<:Poisson}, ::StaticSymbol{:λ}) = asℝ₊
-asparams(::Type{<:Poisson}, ::StaticSymbol{:logλ}) = asℝ
-
 gentype(::Poisson) = Int
 
 Base.rand(rng::AbstractRNG, T::Type, d::Poisson{(:λ,)}) = rand(rng, Dists.Poisson(d.λ))
diff --git a/src/parameterized/snedecorf.jl b/src/parameterized/snedecorf.jl
index 2f0e31a6..e50a530d 100644
--- a/src/parameterized/snedecorf.jl
+++ b/src/parameterized/snedecorf.jl
@@ -16,18 +16,13 @@ end
 
 @inline function basemeasure(d::SnedecorF{(:ν1, :ν2)})
     ℓ = -logbeta(d.ν1 / 2, d.ν2 / 2)
-    weightedmeasure(ℓ, LebesgueMeasure())
+    weightedmeasure(ℓ, LebesgueBase())
 end
 
-xform(::SnedecorF) = asℝ₊
-
 Base.rand(rng::AbstractRNG, T::Type, d::SnedecorF) = rand(rng, proxy(d))
 
 proxy(d::SnedecorF{(:ν1, :ν2)}) = Dists.FDist(d.ν1, d.ν2)
 
-asparams(::Type{<:SnedecorF}, ::StaticSymbol{:ν1}) = asℝ₊
-asparams(::Type{<:SnedecorF}, ::StaticSymbol{:ν2}) = asℝ₊
-
 insupport(::SnedecorF, x) = x > 0
 
 # cdf(d::SnedecorF, x) = StatsFuns.fdistcdf(d.ν1, d.ν2, x)
diff --git a/src/parameterized/studentt.jl b/src/parameterized/studentt.jl
index 7ab7c96c..27a1c825 100644
--- a/src/parameterized/studentt.jl
+++ b/src/parameterized/studentt.jl
@@ -50,20 +50,22 @@ end
 
 @inline function basemeasure(d::StudentT{(:ν,)})
     ℓ = loggamma((d.ν + 1) / 2) - loggamma(d.ν / 2) - log(π * d.ν) / 2
-    weightedmeasure(ℓ, LebesgueMeasure())
+    weightedmeasure(ℓ, LebesgueBase())
 end
 
-xform(::StudentT) = asℝ
-
 Base.rand(rng::AbstractRNG, T::Type, μ::StudentT{(:ν,)}) = rand(rng, Dists.TDist(μ.ν))
 
-proxy(d::StudentT{(:ν,)}) = Dists.TDist(d.ν)
-
 @half StudentT
 
 HalfStudentT(ν, σ) = HalfStudentT((ν = ν, σ = σ))
 
-asparams(::Type{<:StudentT}, ::StaticSymbol{:ν}) = asℝ₊
-
 insupport(::StudentT, x) = true
 insupport(::StudentT) = Returns(true)
+
+proxy(d::StudentT{(:ν,)}) = Dists.TDist(d.ν)
+
+smf(d::StudentT, x) = smf(proxy(d), x)
+invsmf(d::StudentT, p) = invsmf(proxy(d), p)
+
+smf(d::StudentT{(:ν,)}, x) = cdf(proxy(d), x)
+invsmf(d::StudentT{(:ν,)}, p) = quantile(proxy(d), p)
diff --git a/src/parameterized/uniform.jl b/src/parameterized/uniform.jl
index fe19efe4..c14629e4 100644
--- a/src/parameterized/uniform.jl
+++ b/src/parameterized/uniform.jl
@@ -13,16 +13,12 @@ export Uniform
 insupport(::Uniform{()}) = in𝕀
 insupport(::Uniform{()}, x) = in𝕀(x)
 
-@inline basemeasure(::Uniform{()}) = LebesgueMeasure()
-
-proxy(::Uniform{()}) = Dists.Uniform()
+@inline basemeasure(::Uniform{()}) = LebesgueBase()
 
 density_def(::Uniform{()}, x) = 1.0
 
 logdensity_def(d::Uniform{()}, x) = 0.0
 
-xform(::Uniform{()}) = as𝕀
-
 Base.rand(rng::AbstractRNG, T::Type, μ::Uniform{()}) = rand(rng, T)
 
 ###############################################################################
@@ -35,3 +31,8 @@ Uniform(a, b) = Uniform((a = a, b = b))
 proxy(d::Uniform{(:a, :b)}) = affine((μ = d.a, σ = d.b - d.a), Uniform())
 @useproxy Uniform{(:a, :b)}
 Base.rand(rng::Random.AbstractRNG, ::Type{T}, μ::Uniform) where {T} = rand(rng, T, proxy(μ))
+
+smf(::Uniform{()}, x) = smf(StdUniform(), x)
+invsmf(::Uniform{()}, p) = invsmf(StdUniform(), p)
+
+proxy(::Uniform{()}) = Dists.Uniform()
diff --git a/src/smart-constructors.jl b/src/smart-constructors.jl
index ed2796cf..207a26a1 100644
--- a/src/smart-constructors.jl
+++ b/src/smart-constructors.jl
@@ -20,3 +20,5 @@ end
 function affine(f::AffineTransform{(:μ, :λ)}, parent::Lebesgue{RealNumbers})
     affine(AffineTransform((λ = f.λ,)), parent)
 end
+
+affine(f::AffineTransform, μ::LebesgueBase) = weightedmeasure(-logjac(f), μ)
diff --git a/src/transforms/affinetransform.jl b/src/transforms/affinetransform.jl
new file mode 100644
index 00000000..1f11805a
--- /dev/null
+++ b/src/transforms/affinetransform.jl
@@ -0,0 +1,29 @@
+struct OrthoTransform{I,N,T} <: TV.AbstractTransform
+    dim::I
+    par::NamedTuple{N,T}
+end
+
+TV.dimension(t::OrthoTransform) = t.dim
+
+TV.inverse_eltype(t::OrthoTransform, y::AbstractVector) = extended_eltype(y)
+
+function TV.inverse_at!(
+    x::AbstractVector,
+    index,
+    t::OrthoTransform{(:σ,)},
+    y::AbstractVector,
+) end
+
+function TV.transform_with(
+    flag::LogJacFlag,
+    t::OrthoTransform{(:σ,)},
+    x::AbstractVector,
+    index,
+)
+    ylen, xlen = size(d.σ)
+    T = extended_eltype(x)
+    ℓ = logjac_zero(flag, T)
+    y = d.σ * view(x, index, index + xlen - 1)
+    index += xlen
+    y, ℓ, index
+end
diff --git a/src/transforms/corrcholesky.jl b/src/transforms/corrcholesky.jl
index 9b3d39a0..47624f7a 100644
--- a/src/transforms/corrcholesky.jl
+++ b/src/transforms/corrcholesky.jl
@@ -30,9 +30,9 @@ function TV.transform_with(flag::TV.LogJacFlag, t::CorrCholesky, x::AbstractVect
     return Cholesky(U, 'U', 0), ℓ, index
 end
 
-TV.inverse_eltype(t::CorrCholesky, x::AbstractMatrix) = TV.extended_eltype(x)
+TV.inverse_eltype(t::CorrCholesky, x::AbstractMatrix) = eltype(x)
 
-TV.inverse_eltype(t::CorrCholesky, x::Cholesky) = TV.extended_eltype(x)
+TV.inverse_eltype(t::CorrCholesky, x::Cholesky) = eltype(x)
 
 function TV.inverse_at!(x::AbstractVector, index, t::CorrCholesky, L::LowerTriangular)
     return TV.inverse_at!(x, index, CorrCholeskyUpper(t.n), L')
diff --git a/src/transforms/corrcholeskylower.jl b/src/transforms/corrcholeskylower.jl
index e90fda06..7f457869 100644
--- a/src/transforms/corrcholeskylower.jl
+++ b/src/transforms/corrcholeskylower.jl
@@ -38,7 +38,7 @@ function TV.transform_with(
     return U', ℓ, index
 end
 
-TV.inverse_eltype(t::CorrCholeskyLower, L::LowerTriangular) = TV.extended_eltype(L)
+TV.inverse_eltype(t::CorrCholeskyLower, L::LowerTriangular) = eltype(L)
 
 function TV.inverse_at!(x::AbstractVector, index, t::CorrCholeskyLower, L::LowerTriangular)
     return TV.inverse_at!(x, index, CorrCholeskyUpper(t.n), L')
diff --git a/src/transforms/ordered.jl b/src/transforms/ordered.jl
index 7c927a66..72a641aa 100644
--- a/src/transforms/ordered.jl
+++ b/src/transforms/ordered.jl
@@ -49,7 +49,7 @@ function TV.transform_with(flag::TV.LogJacFlag, t::Ordered, x, index::T) where {
     return (y, ℓ, index)
 end
 
-TV.inverse_eltype(t::Ordered, y::AbstractVector) = TV.extended_eltype(y)
+TV.inverse_eltype(t::Ordered, y::AbstractVector) = eltype(y)
 
 Ordered(n::Int) = Ordered(asℝ, n)
 
diff --git a/src/utils.jl b/src/utils.jl
index dbd19f78..4348ccbf 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1,5 +1,4 @@
 using LinearAlgebra
-export mydot
 
 function solve(A::Union{AbstractMatrix,Factorization}, y::AbstractArray)
     (m, n) = size(A)
@@ -120,3 +119,15 @@ function getU(C::Cholesky)
 end
 
 const Triangular = Union{L,U} where {L<:LowerTriangular,U<:UpperTriangular}
+
+# Ideally we could just add one method smf(μ::$M, x::Dual{TAG}) where TAG
+# But then we have method ambiguities
+macro smfAD(M)
+    quote
+        function MeasureBase.smf(μ::$M, x::Dual{TAG}) where {TAG}
+            val = ForwardDiff.value(x)
+            Δ = ForwardDiff.partials(x)
+            Dual{TAG}(smf(μ, val), Δ * densityof(μ, val))
+        end
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 0bee927a..de93cebc 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,7 +3,7 @@ using StatsFuns
 using Base.Iterators: take
 using Random
 using LinearAlgebra
-using DynamicIterators: trace, TimeLift
+# using DynamicIterators: trace, TimeLift
 using TransformVariables: transform, as𝕀
 using FillArrays
 
@@ -73,28 +73,15 @@ test_measures = Any[
     MvNormal(σ = σ)
     MvNormal(λ = λ)
     Uniform()
-    Counting(Float64)
     Dirac(0.0) + Normal()
 ]
 
-testbroken_measures = Any[
-    Pushforward(as𝕀, Normal())
-    # InverseGamma(2) # Not defined yet
-    # MvNormal(I(3)) # Entirely broken for now
-    TrivialMeasure()
-]
-
 @testset "testvalue" begin
     for μ in test_measures
         @info "testing $μ"
         test_interface(μ)
     end
 
-    for μ in testbroken_measures
-        @info "testing $μ"
-        @test_broken test_measure(μ)
-    end
-
     @testset "testvalue(::Chain)" begin
         mc = Chain(x -> Normal(μ = x), Normal(μ = 0.0))
         r = testvalue(mc)
@@ -120,7 +107,7 @@ end
         @test rand(rng, Binomial(n=0, p=1.0)) == 0
         @test rand(rng, Binomial(n=10, p=1.0)) == 10
 
-        @test_broken logdensity_def(Binomial(n, p), CountingMeasure(ℤ[0:n]), x) ≈
+        @test_broken logdensity_def(Binomial(n, p), CountingBase(ℤ[0:n]), x) ≈
                      binomlogpdf(n, p, x)
     end
 
@@ -157,7 +144,7 @@ end
         sample2 = rand(MersenneTwister(123), NegativeBinomial(; r, logλ))
         @test sample1 == sample2
 
-        @test_broken logdensity_def(Binomial(n, p), CountingMeasure(ℤ[0:n]), x) ≈
+        @test_broken logdensity_def(Binomial(n, p), CountingBase(ℤ[0:n]), x) ≈
                      binomlogpdf(n, p, x)
     end
 
@@ -167,25 +154,24 @@ end
         @test sample1 == sample2
     end
 
-    # Fails because we need `asparams` for `::AffinePushfwd`
-    # @testset "Normal" begin
-    #     D = affine{(:μ,:σ), Normal}
-    #     par = transform(asparams(D), randn(2))
-    #     d = D(par)
-    #     @test params(d) == par
+    @testset "Normal" begin
+        D = Normal{(:μ, :σ)}
+        par = transform(asparams(D), randn(2))
+        d = D(par)
+        @test params(d) == par
 
-    #     μ = par.μ
-    #     σ = par.σ
-    #     σ² = σ^2
-    #     τ = 1/σ²
-    #     logσ = log(σ)
-    #     y = rand(d)
+        μ = par.μ
+        σ = par.σ
+        σ² = σ^2
+        τ = 1 / σ²
+        logσ = log(σ)
+        y = rand(d)
 
-    #     ℓ = logdensity_def(Normal(;μ,σ), y)
-    #     @test ℓ ≈ logdensity_def(Normal(;μ,σ²), y)
-    #     @test ℓ ≈ logdensity_def(Normal(;μ,τ), y)
-    #     @test ℓ ≈ logdensity_def(Normal(;μ,logσ), y)
-    # end
+        ℓ = logdensityof(Normal(; μ, σ), y)
+        @test ℓ ≈ logdensityof(Normal(; μ, σ²), y)
+        @test ℓ ≈ logdensityof(Normal(; μ, τ), y)
+        @test ℓ ≈ logdensityof(Normal(; μ, logσ), y)
+    end
 
     @testset "LKJCholesky" begin
         D = LKJCholesky{(:k, :η)}
@@ -288,6 +274,8 @@ end
     @test transform(t, []) == x
 end
 
+using DynamicIterators: trace, TimeLift
+
 # @testset "Univariate chain" begin
 #     ξ0 = 1.
 #     x = 1.2
@@ -313,7 +301,7 @@ end
 
 @testset "rootmeasure/logpdf" begin
     x = rand(Normal())
-    @test logdensityof(𝒹(Normal(), rootmeasure(Normal())), x) ≈ logdensityof(Normal(), x)
+    @test logdensity_rel(Normal(), rootmeasure(Normal()), x) ≈ logdensityof(Normal(), x)
 end
 
 @testset "Transforms" begin
@@ -430,7 +418,7 @@ end
     end
 
     @testset "Normal" begin
-        @test_broken repro(Normal, (:μ, :σ))
+        @test repro(Normal, (:μ, :σ))
     end
 
     @testset "Poisson" begin
@@ -480,45 +468,26 @@ end
     end
 end
 
-@testset "Density measures and Radon-Nikodym" begin
-    x = randn()
-    let d = ∫(𝒹(Cauchy(), Normal()), Normal())
-        @test logdensityof(𝒹(d, Cauchy()), x) ≈ 0 atol = 1e-12
-    end
-
-    let f = 𝒹(∫(x -> x^2, Normal()), Normal())
-        @test densityof(f, x) ≈ x^2
-    end
-
-    # let d = ∫exp(log𝒹(Cauchy(), Normal()), Normal())
-    #     @test logdensity_def(d, Cauchy(), x) ≈ 0 atol=1e-12
-    # end
-
-    let f = 𝒹(∫exp(x -> x^2, Normal()), Normal())
-        @test logdensityof(f, x) ≈ x^2
-    end
-end
-
 @testset "Half measures" begin
     @testset "HalfNormal" begin
         d = Normal(σ = 3)
         h = HalfNormal(3)
         x = rand(h)
-        @test densityof(𝒹(h, Lebesgue(ℝ)), x) ≈ 2 * densityof(𝒹(d, Lebesgue(ℝ)), x)
+        @test densityof(h, x) ≈ 2 * densityof(d, x)
     end
 
     @testset "HalfCauchy" begin
         d = Cauchy(σ = 3)
         h = HalfCauchy(3)
         x = rand(h)
-        @test densityof(𝒹(h, Lebesgue(ℝ)), x) ≈ 2 * densityof(𝒹(d, Lebesgue(ℝ)), x)
+        @test densityof(h, x) ≈ 2 * densityof(d, x)
     end
 
     @testset "HalfStudentT" begin
         d = StudentT(ν = 2, σ = 3)
         h = HalfStudentT(2, 3)
         x = rand(h)
-        @test densityof(𝒹(h, Lebesgue(ℝ)), x) ≈ 2 * densityof(𝒹(d, Lebesgue(ℝ)), x)
+        @test densityof(h, x) ≈ 2 * densityof(d, x)
     end
 end
 
@@ -667,7 +636,39 @@ end
     @test logdensityof(d, x) isa Real
 end
 
-@testset "Distributions.jl cdf" begin
-    @test cdf(Normal(0, 1), 0) == 0.5
-    @test cdf.((Normal(0, 1),), [0, 0]) == [0.5, 0.5]
+@testset "Normal smf" begin
+    @test smf(Normal(0, 1), 0) == 0.5
+    @test smf.((Normal(0, 1),), [0, 0]) == [0.5, 0.5]
+end
+
+affinepars_1d = [
+    (;)
+    (μ = randn(),)
+    (σ = randn(),)
+    (λ = randn(),)
+    (μ = randn(), σ = randn())
+    (μ = randn(), λ = randn())
+]
+
+smf_measures =
+    [
+        [Bernoulli(), Bernoulli(rand()), Bernoulli(logitp = randn())]
+        [Beta(rand(), rand())]
+        # [BetaBinomial(rand(5:20), 2 * rand(), 2 * rand())]
+        # [Binomial(rand(5:20), rand())]
+        Cauchy.(affinepars_1d)
+        Normal.(affinepars_1d)
+        StudentT.(merge.(Ref((ν = 1 + 2 * rand(),)), affinepars_1d))
+    ] |> vcat
+
+@testset "smf" begin
+    for μ in smf_measures
+        test_smf(μ)
+    end
+end
+
+@testset "more interface tests" begin
+    for μ in smf_measures
+        test_interface(μ)
+    end
 end