diff --git a/Project.toml b/Project.toml
index 0612503698..f90a728cee 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Manifolds"
 uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
 authors = ["Seth Axen <seth.axen@gmail.com>", "Mateusz Baran <mateuszbaran89@gmail.com>", "Ronny Bergmann <manopt@ronnybergmann.net>", "Antoine Levitt <antoine.levitt@gmail.com>"]
-version = "0.8.74"
+version = "0.8.75"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
diff --git a/docs/Project.toml b/docs/Project.toml
index ed58b3d015..a253c731c4 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -3,6 +3,7 @@ BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
@@ -21,6 +22,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 BoundaryValueDiffEq = "2"
 CondaPkg = "0.2"
 DiffEqCallbacks = "2"
+Distributions = "0.22.6, 0.23, 0.24, 0.25"
 Documenter = "0.27"
 FiniteDifferences = "0.12"
 Graphs = "1.4"
diff --git a/docs/make.jl b/docs/make.jl
index bc38f2c19c..7417755f1c 100755
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -91,6 +91,7 @@ makedocs(
             "🚀 Get Started with `Manifolds.jl`" => "tutorials/getstarted.md",
             "work in charts" => "tutorials/working-in-charts.md",
             "perform Hand gesture analysis" => "tutorials/hand-gestures.md",
+            "integrate on manifolds and handle probability densities" => "tutorials/integration.md",
         ],
         "Manifolds" => [
             "Basic manifolds" => [
@@ -148,6 +149,7 @@ makedocs(
             "Atlases and charts" => "features/atlases.md",
             "Differentiation" => "features/differentiation.md",
             "Distributions" => "features/distributions.md",
+            "Integration" => "features/integration.md",
             "Statistics" => "features/statistics.md",
             "Testing" => "features/testing.md",
             "Utilities" => "features/utilities.md",
diff --git a/docs/src/features/integration.md b/docs/src/features/integration.md
new file mode 100644
index 0000000000..5249ea0f3c
--- /dev/null
+++ b/docs/src/features/integration.md
@@ -0,0 +1,6 @@
+# Integration
+
+```@docs
+manifold_volume(::AbstractManifold)
+volume_density(::AbstractManifold, ::Any, ::Any)
+```
diff --git a/docs/src/manifolds/probabilitysimplex.md b/docs/src/manifolds/probabilitysimplex.md
index 2242802562..19237f6491 100644
--- a/docs/src/manifolds/probabilitysimplex.md
+++ b/docs/src/manifolds/probabilitysimplex.md
@@ -8,6 +8,15 @@ Private=false
 Public=true
 ```
 
+```@autodocs
+Modules = [Manifolds]
+Pages = ["manifolds/ProbabilitySimplexEuclideanMetric.jl"]
+Order = [:type, :function]
+Private=false
+Public=true
+```
+
+
 ## Real probability amplitudes
 
 An isometric embedding of interior of [`ProbabilitySimplex`](@ref) in positive orthant of the
diff --git a/docs/src/misc/notation.md b/docs/src/misc/notation.md
index c0020f5bb1..67fcf7259c 100644
--- a/docs/src/misc/notation.md
+++ b/docs/src/misc/notation.md
@@ -53,4 +53,6 @@ Within the documented functions, the utf8 symbols are used whenever possible, as
 | ``B`` | a vector bundle | |
 | ``\mathcal T_{q\gets p}X`` | vector transport | | of the vector ``X`` from ``T_p\mathcal M`` to ``T_q\mathcal M``
 | ``\mathcal T_{p,Y}X`` | vector transport in direction ``Y`` | | of the vector ``X`` from ``T_p\mathcal M`` to ``T_q\mathcal M``, where ``q`` is deretmined by ``Y``, for example using the exponential map or some retraction.
+| ``\operatorname{Vol}(\mathcal M)`` | volume of manifold ``\mathcal M`` | |
+| ``\theta_p(X)`` | volume density for vector ``X`` tangent at point ``p`` | |
 | ``0_k`` | the ``k\times k`` zero matrix. | |
diff --git a/src/Manifolds.jl b/src/Manifolds.jl
index 647e2f2717..cb0001b925 100644
--- a/src/Manifolds.jl
+++ b/src/Manifolds.jl
@@ -319,6 +319,7 @@ using Random
 using RecursiveArrayTools: ArrayPartition
 using Requires
 using SimpleWeightedGraphs: AbstractSimpleWeightedGraph, get_weight
+using SpecialFunctions
 using StaticArrays
 using Statistics
 using StatsBase
@@ -487,6 +488,37 @@ function Base.in(X, TpM::TangentSpaceAtPoint; kwargs...)
     return is_vector(base_manifold(TpM), TpM.point, X, false; kwargs...)
 end
 
+@doc raw"""
+    manifold_volume(M::AbstractManifold)
+
+Volume of manifold `M` defined through integration of Riemannian volume element in a chart.
+Note that for many manifolds there is no universal agreement over the exact ranges over
+which the integration should happen. For details see [^BoyaSudarshanTilma2003].
+
+[^BoyaSudarshanTilma2003]:
+    > L. J. Boya, E. C. G. Sudarshan, and T. Tilma, “Volumes of Compact Manifolds,” Reports
+    > on Mathematical Physics, vol. 52, no. 3, pp. 401–422, Dec. 2003,
+    > doi: [10.1016/S0034-4877(03)80038-1](https://doi.org/10.1016/S0034-4877(03)80038-1)
+"""
+manifold_volume(::AbstractManifold)
+
+@doc raw"""
+    volume_density(M::AbstractManifold, p, X)
+
+Volume density function of manifold `M`, i.e. determinant of the differential of exponential map
+`exp(M, p, X)`. Determinant can be understood as computed in a basis, from the matrix
+of the linear operator said differential corresponds to. Details are available in Section 4.1
+of [^ChevallierLiLuDunson2022].
+
+Note that volume density is well-defined only for `X` for which `exp(M, p, X)` is injective.
+
+[^ChevallierLiLuDunson2022]:
+    > E. Chevallier, D. Li, Y. Lu, and D. B. Dunson, “Exponential-wrapped distributions on
+    > symmetric spaces.” arXiv, Oct. 09, 2022.
+    > doi: [10.48550/arXiv.2009.01983](https://doi.org/10.48550/arXiv.2009.01983).
+"""
+volume_density(::AbstractManifold, p, X)
+
 # functions populated with methods by extensions
 
 function solve_chart_log_bvp end
@@ -789,6 +821,7 @@ export ×,
     log!,
     log_local_metric_density,
     manifold_dimension,
+    manifold_volume,
     metric,
     mean,
     mean!,
@@ -852,6 +885,7 @@ export ×,
     vee!,
     vertical_component,
     vertical_component!,
+    volume_density,
     zero_vector,
     zero_vector!
 # Lie group types & functions
diff --git a/src/groups/general_unitary_groups.jl b/src/groups/general_unitary_groups.jl
index 04ec75eb41..c56311b86c 100644
--- a/src/groups/general_unitary_groups.jl
+++ b/src/groups/general_unitary_groups.jl
@@ -279,6 +279,10 @@ function log_lie!(
     return project!(G, X, Identity(G), X)
 end
 
+function manifold_volume(M::GeneralUnitaryMultiplicationGroup)
+    return manifold_volume(M.manifold)
+end
+
 function Random.rand!(G::GeneralUnitaryMultiplicationGroup, pX; kwargs...)
     rand!(G.manifold, pX; kwargs...)
     return pX
@@ -328,3 +332,7 @@ function translate_diff!(
 )
     return copyto!(G, Y, inv(G, p) * X * p)
 end
+
+function volume_density(M::GeneralUnitaryMultiplicationGroup, p, X)
+    return volume_density(M.manifold, p, X)
+end
diff --git a/src/manifolds/Circle.jl b/src/manifolds/Circle.jl
index 7b35513dbd..b6fa91c327 100644
--- a/src/manifolds/Circle.jl
+++ b/src/manifolds/Circle.jl
@@ -331,6 +331,13 @@ i.e. $\dim(𝕊^1) = 1$.
 """
 manifold_dimension(::Circle) = 1
 
+@doc raw"""
+    manifold_volume(M::Circle)
+
+Return the volume of the [`Circle`](@ref) `M`, i.e. ``2π``.
+"""
+manifold_volume(::Circle) = 2 * π
+
 @doc raw"""
     mean(M::Circle{ℝ}, x::AbstractVector[, w::AbstractWeights])
 
@@ -525,5 +532,12 @@ function parallel_transport_to!(M::Circle{ℂ}, Y, p, X, q)
     return Y
 end
 
+"""
+    volume_density(::Circle, p, X)
+
+Return volume density of [`Circle`](@ref), i.e. 1.
+"""
+volume_density(::Circle, p, X) = one(eltype(X))
+
 zero_vector(::Circle, p::T) where {T<:Number} = zero(p)
 zero_vector!(::Circle, X, p) = fill!(X, 0)
diff --git a/src/manifolds/Euclidean.jl b/src/manifolds/Euclidean.jl
index ce2eef2afe..cd3b663aa1 100644
--- a/src/manifolds/Euclidean.jl
+++ b/src/manifolds/Euclidean.jl
@@ -454,6 +454,13 @@ function manifold_dimension(M::Euclidean{N,𝔽}) where {N,𝔽}
 end
 manifold_dimension(::Euclidean{Tuple{},𝔽}) where {𝔽} = real_dimension(𝔽)
 
+"""
+    manifold_volume(::Euclidean)
+
+Return volume of the [`Euclidean`](@ref) manifold, i.e. infinity.
+"""
+manifold_volume(::Euclidean) = Inf
+
 Statistics.mean(::Euclidean{Tuple{}}, x::AbstractVector{<:Number}; kwargs...) = mean(x)
 function Statistics.mean(
     ::Euclidean{Tuple{}},
@@ -709,6 +716,15 @@ function Statistics.var(::Euclidean, x::AbstractVector{<:Number}, m::Number; kwa
     return sum(var(x; mean=m, kwargs...))
 end
 
+@doc raw"""
+    volume_density(M::Euclidean, p, X)
+
+Return volume density function of [`Euclidean`](@ref) manifold `M`, i.e. 1.
+"""
+function volume_density(::Euclidean, p, X)
+    return one(eltype(X))
+end
+
 """
     zero_vector(M::Euclidean, x)
 
diff --git a/src/manifolds/GeneralUnitaryMatrices.jl b/src/manifolds/GeneralUnitaryMatrices.jl
index 494f3a485a..d2abb58a19 100644
--- a/src/manifolds/GeneralUnitaryMatrices.jl
+++ b/src/manifolds/GeneralUnitaryMatrices.jl
@@ -679,6 +679,122 @@ Return the dimension of the manifold of special unitary matrices.
 """
 manifold_dimension(::GeneralUnitaryMatrices{n,ℂ,DeterminantOneMatrices}) where {n} = n^2 - 1
 
+@doc raw"""
+    manifold_volume(::GeneralUnitaryMatrices{n,ℝ,AbsoluteDeterminantOneMatrices}) where {n}
+
+Volume of the manifold of real orthogonal matrices of absolute determinant one. The
+formula reads [^BoyaSudarshanTilma2003]:
+
+```math
+\begin{cases}
+\frac{2^{k}(2\pi)^{k^2}}{\prod_{s=1}^{k-1} (2s)!} & \text{ if } n = 2k \\
+\frac{2^{k+1}(2\pi)^{k(k+1)}}{\prod_{s=1}^{k-1} (2s+1)!} & \text{ if } n = 2k+1
+\end{cases}
+```
+
+[^BoyaSudarshanTilma2003]:
+    > L. J. Boya, E. C. G. Sudarshan, and T. Tilma, “Volumes of Compact Manifolds,” Reports
+    > on Mathematical Physics, vol. 52, no. 3, pp. 401–422, Dec. 2003,
+    > doi: [10.1016/S0034-4877(03)80038-1](https://doi.org/10.1016/S0034-4877(03)80038-1)
+"""
+function manifold_volume(
+    ::GeneralUnitaryMatrices{n,ℝ,AbsoluteDeterminantOneMatrices},
+) where {n}
+    return 2 * manifold_volume(GeneralUnitaryMatrices{n,ℝ,DeterminantOneMatrices}())
+end
+@doc raw"""
+    manifold_volume(::GeneralUnitaryMatrices{n,ℝ,DeterminantOneMatrices}) where {n}
+
+Volume of the manifold of real orthogonal matrices of determinant one. The
+formula reads [^BoyaSudarshanTilma2003]:
+
+```math
+\begin{cases}
+2 & \text{ if } n = 0 \\
+\frac{2^{k-1/2}(2\pi)^{k^2}}{\prod_{s=1}^{k-1} (2s)!} & \text{ if } n = 2k+2 \\
+\frac{2^{k+1/2}(2\pi)^{k(k+1)}}{\prod_{s=1}^{k-1} (2s+1)!} & \text{ if } n = 2k+1
+\end{cases}
+```
+
+It differs from the paper by a factor of `sqrt(2)` due to a different choice of
+normalization.
+
+[^BoyaSudarshanTilma2003]:
+    > L. J. Boya, E. C. G. Sudarshan, and T. Tilma, “Volumes of Compact Manifolds,” Reports
+    > on Mathematical Physics, vol. 52, no. 3, pp. 401–422, Dec. 2003,
+    > doi: [10.1016/S0034-4877(03)80038-1](https://doi.org/10.1016/S0034-4877(03)80038-1)
+"""
+function manifold_volume(::GeneralUnitaryMatrices{n,ℝ,DeterminantOneMatrices}) where {n}
+    vol = 1.0
+    if n % 2 == 0
+        k = div(n, 2)
+        vol *= 2^(k - 1) * (2π)^(k^2)
+        for s in 1:(k - 1)
+            vol /= factorial(2 * s)
+        end
+    else
+        k = div(n - 1, 2)
+        vol *= 2^k * (2π)^(k * (k + 1))
+        for s in 1:(k - 1)
+            vol /= factorial(2 * s + 1)
+        end
+    end
+    if n > 1
+        vol *= sqrt(2)
+    end
+    return vol
+end
+@doc raw"""
+    manifold_volume(::GeneralUnitaryMatrices{n,ℂ,AbsoluteDeterminantOneMatrices}) where {n}
+
+Volume of the manifold of complex general unitary matrices of absolute determinant one. The
+formula reads [^BoyaSudarshanTilma2003]:
+
+```math
+\sqrt{n 2^{n+1}} π^{n(n+1)/2} \prod_{k=1}^{n-1}\frac{1}{k!}
+```
+
+[^BoyaSudarshanTilma2003]:
+    > L. J. Boya, E. C. G. Sudarshan, and T. Tilma, “Volumes of Compact Manifolds,” Reports
+    > on Mathematical Physics, vol. 52, no. 3, pp. 401–422, Dec. 2003,
+    > doi: [10.1016/S0034-4877(03)80038-1](https://doi.org/10.1016/S0034-4877(03)80038-1)
+"""
+function manifold_volume(
+    ::GeneralUnitaryMatrices{n,ℂ,AbsoluteDeterminantOneMatrices},
+) where {n}
+    vol = sqrt(n * 2^(n + 1)) * π^(((n + 1) * n) // 2)
+    kf = 1
+    for k in 1:(n - 1)
+        kf *= k
+        vol /= kf
+    end
+    return vol
+end
+@doc raw"""
+    manifold_volume(::GeneralUnitaryMatrices{n,ℂ,DeterminantOneMatrices}) where {n}
+
+Volume of the manifold of complex general unitary matrices of determinant one. The formula
+reads [^BoyaSudarshanTilma2003]:
+
+```math
+\sqrt{n 2^{n-1}} π^{(n-1)(n+2)/2} \prod_{k=1}^{n-1}\frac{1}{k!}
+```
+
+[^BoyaSudarshanTilma2003]:
+    > L. J. Boya, E. C. G. Sudarshan, and T. Tilma, “Volumes of Compact Manifolds,” Reports
+    > on Mathematical Physics, vol. 52, no. 3, pp. 401–422, Dec. 2003,
+    > doi: [10.1016/S0034-4877(03)80038-1](https://doi.org/10.1016/S0034-4877(03)80038-1)
+"""
+function manifold_volume(::GeneralUnitaryMatrices{n,ℂ,DeterminantOneMatrices}) where {n}
+    vol = sqrt(n * 2^(n - 1)) * π^(((n - 1) * (n + 2)) // 2)
+    kf = 1
+    for k in 1:(n - 1)
+        kf *= k
+        vol /= kf
+    end
+    return vol
+end
+
 """
     mean(
         M::Rotations,
@@ -805,3 +921,68 @@ function riemann_tensor!(::GeneralUnitaryMatrices, Xresult, p, X, Y, Z)
     Xresult .= 1 // 4 .* (Z * Xtmp .- Xtmp * Z)
     return Xresult
 end
+
+@doc raw"""
+    volume_density(M::GeneralUnitaryMatrices{n,ℝ}, p, X) where {n}
+
+Compute volume density function of a sphere, i.e. determinant of the differential of
+exponential map `exp(M, p, X)`. It is derived from Eq. (4.1) and Corollary 4.4
+in [^ChevallierLiLuDunson2022]. See also Theorem 4.1 in [^FalorsideHaanDavidsonForré2019],
+(note that it uses a different convention).
+
+[^ChevallierLiLuDunson2022]:
+    > E. Chevallier, D. Li, Y. Lu, and D. B. Dunson, “Exponential-wrapped distributions on
+    > symmetric spaces.” arXiv, Oct. 09, 2022.
+    > doi: [10.48550/arXiv.2009.01983](https://doi.org/10.48550/arXiv.2009.01983).
+
+[^FalorsideHaanDavidsonForré2019]:
+    > L. Falorsi, P. de Haan, T. R. Davidson, and P. Forré, “Reparameterizing Distributions
+    > on Lie Groups,” arXiv:1903.02958 [cs, math, stat], Mar. 2019
+    > doi: [10.48550/arXiv.1903.02958](https://doi.org/10.48550/arXiv.1903.02958)
+"""
+function volume_density(M::GeneralUnitaryMatrices{n,ℝ}, p, X) where {n}
+    dens = one(eltype(X))
+    B = get_basis(M, p, DefaultOrthonormalBasis())
+    Ys = get_vectors(M, p, B)
+    Z = similar(X)
+    op_coeffs = similar(X, manifold_dimension(M), manifold_dimension(M))
+    for k in 1:manifold_dimension(M)
+        Y = Ys[k]
+        Z .= X * Y .- Y * X
+        get_coordinates!(M, view(op_coeffs, :, k), p, Z, DefaultOrthonormalBasis())
+    end
+    for ev in eigvals(op_coeffs)
+        if abs(ev) > eps(eltype(X))
+            cm = (1 - exp(-ev)) / ev
+            dens *= real(cm)
+        end
+    end
+
+    return dens
+end
+
+@doc raw"""
+    volume_density(M::GeneralUnitaryMatrices{3,ℝ}, p, X)
+
+Compute the volume density on O(3)/SO(3). The formula reads [^FalorsideHaanDavidsonForré2019]:
+```math
+\frac{1-1\cos(\sqrt{2}\lVert X \rVert)}{\lVert X \rVert^2}.
+```
+"""
+function volume_density(M::GeneralUnitaryMatrices{3,ℝ}, p, X)
+    nX = norm(M, p, X)
+    if nX > eps(eltype(X))
+        return (1 - 1 * cos(sqrt(2) * nX)) / nX^2
+    else
+        return one(nX)
+    end
+end
+
+@doc raw"""
+    volume_density(M::GeneralUnitaryMatrices{2,ℝ}, p, X)
+
+Volume density on O(2)/SO(2) is equal to 1.
+"""
+function volume_density(::GeneralUnitaryMatrices{2,ℝ}, p, X)
+    return one(eltype(X))
+end
diff --git a/src/manifolds/Hyperbolic.jl b/src/manifolds/Hyperbolic.jl
index d563c07614..9e18b5ec22 100644
--- a/src/manifolds/Hyperbolic.jl
+++ b/src/manifolds/Hyperbolic.jl
@@ -299,6 +299,13 @@ Return the dimension of the hyperbolic space manifold $\mathcal H^n$, i.e. $\dim
 """
 manifold_dimension(::Hyperbolic{N}) where {N} = N
 
+@doc raw"""
+    manifold_dimension(M::Hyperbolic)
+
+Return the volume of the hyperbolic space manifold $\mathcal H^n$, i.e. infinity.
+"""
+manifold_volume(::Hyperbolic) = Inf
+
 """
     mean(
         M::Hyperbolic,
diff --git a/src/manifolds/HyperbolicHyperboloid.jl b/src/manifolds/HyperbolicHyperboloid.jl
index a751477d48..681fc561db 100644
--- a/src/manifolds/HyperbolicHyperboloid.jl
+++ b/src/manifolds/HyperbolicHyperboloid.jl
@@ -432,3 +432,25 @@ function parallel_transport_to!(::Hyperbolic, Y, p, X, q)
         X .+ minkowski_metric(q, X) ./ (1 - minkowski_metric(p, q)) .* (p + q),
     )
 end
+
+@doc raw"""
+    volume_density(M::Hyperbolic, p, X)
+
+Compute volume density function of the hyperbolic manifold. The formula reads
+``(\sinh(\lVert X\rVert)/\lVert X\rVert)^(n-1)`` where `n` is the dimension of `M`.
+It is derived from Eq. (4.1) in [^ChevallierLiLuDunson2022].
+
+[^ChevallierLiLuDunson2022]:
+    > E. Chevallier, D. Li, Y. Lu, and D. B. Dunson, “Exponential-wrapped distributions on
+    > symmetric spaces.” arXiv, Oct. 09, 2022.
+    > doi: [10.48550/arXiv.2009.01983](https://doi.org/10.48550/arXiv.2009.01983).
+"""
+function volume_density(M::Hyperbolic, p, X)
+    Xnorm = norm(X)
+    if Xnorm == 0
+        return one(eltype(X))
+    else
+        n = manifold_dimension(M) - 1
+        return (sinh(Xnorm) / Xnorm)^n
+    end
+end
diff --git a/src/manifolds/PositiveNumbers.jl b/src/manifolds/PositiveNumbers.jl
index cf9ad5b89b..6196e2b4ec 100644
--- a/src/manifolds/PositiveNumbers.jl
+++ b/src/manifolds/PositiveNumbers.jl
@@ -234,6 +234,15 @@ i.e. of the 1-dimensional hyperbolic space,
 """
 manifold_dimension(::PositiveNumbers) = 1
 
+@doc raw"""
+    manifold_volume(M::PositiveNumbers)
+
+Return volume of [`PositiveNumbers`](@ref) `M`, i.e. `Inf`.
+"""
+function manifold_volume(::PositiveNumbers)
+    return Inf
+end
+
 mid_point(M::PositiveNumbers, p1, p2) = exp(M, p1, log(M, p1, p2) / 2)
 
 @inline LinearAlgebra.norm(::PositiveNumbers, p, X) = sum(abs.(X ./ p))
@@ -311,5 +320,17 @@ function vector_transport_direction(
     return vector_transport_to(M, p, X, q, m)
 end
 
+@doc raw"""
+    volume_density(M::PositiveNumbers, p, X)
+
+Compute volume density function of [`PositiveNumbers`](@ref). The formula reads
+```math
+\theta_p(X) = \exp(X / p)
+```
+"""
+function volume_density(::PositiveNumbers, p, X)
+    return exp(X / p)
+end
+
 zero_vector(::PositiveNumbers, p::Real) = zero(p)
 zero_vector!(::PositiveNumbers, X, p) = fill!(X, 0)
diff --git a/src/manifolds/PowerManifold.jl b/src/manifolds/PowerManifold.jl
index bc5c3d6a9b..43cf98aef4 100644
--- a/src/manifolds/PowerManifold.jl
+++ b/src/manifolds/PowerManifold.jl
@@ -177,6 +177,15 @@ Base.@propagate_inbounds function Base.getindex(
     return get_component(M, p, I...)
 end
 
+@doc raw"""
+    manifold_volume(M::PowerManifold)
+
+Return the manifold volume of an [`PowerManifold`](@ref) `M`.
+"""
+function manifold_volume(M::PowerManifold{𝔽,<:AbstractManifold,TSize}) where {𝔽,TSize}
+    return manifold_volume(M.manifold)^prod(size_to_tuple(TSize))
+end
+
 function Random.rand(rng::AbstractRNG, d::PowerFVectorDistribution)
     fv = zero_vector(d.type, d.point)
     Distributions._rand!(rng, d, fv)
@@ -294,6 +303,23 @@ function vector_bundle_transport(fiber::VectorSpaceType, M::PowerManifold)
     return ParallelTransport()
 end
 
+@doc raw"""
+    volume_density(M::PowerManifold, p, X)
+
+Return volume density on the [`PowerManifold`](@ref) `M`, i.e. product of constituent
+volume densities.
+"""
+function volume_density(M::PowerManifold, p, X)
+    density = one(float(eltype(X)))
+    rep_size = representation_size(M.manifold)
+    for i in get_iterator(M)
+        p_i = _read(M, rep_size, p, i)
+        X_i = _read(M, rep_size, X, i)
+        density *= volume_density(M.manifold, p_i, X_i)
+    end
+    return density
+end
+
 @inline function _write(
     ::PowerManifoldMultidimensional,
     rep_size::Tuple,
diff --git a/src/manifolds/ProbabilitySimplex.jl b/src/manifolds/ProbabilitySimplex.jl
index 1580401580..5a1d2c5ec6 100644
--- a/src/manifolds/ProbabilitySimplex.jl
+++ b/src/manifolds/ProbabilitySimplex.jl
@@ -327,6 +327,15 @@ Returns the manifold dimension of the probability simplex in $ℝ^{n+1}$, i.e.
 """
 manifold_dimension(::ProbabilitySimplex{n}) where {n} = n
 
+@doc raw"""
+    manifold_volume(::ProbabilitySimplex{n}) where {n}
+
+Return the volume of the [`ProbabilitySimplex`](@ref), i.e. volume of the `n`-dimensional
+[`Sphere`](@ref) divided by ``2^{n+1}``, corresponding to the volume of its positive
+orthant.
+"""
+manifold_volume(::ProbabilitySimplex{n}) where {n} = manifold_volume(Sphere(n)) / 2^(n + 1)
+
 @doc raw"""
     mean(
         M::ProbabilitySimplex,
@@ -516,6 +525,17 @@ function Base.show(io::IO, ::ProbabilitySimplex{n,boundary}) where {n,boundary}
     return print(io, "ProbabilitySimplex($(n); boundary=:$boundary)")
 end
 
+@doc raw"""
+    volume_density(M::ProbabilitySimplex{N}, p, X) where {N}
+
+Compute the volume density at point `p` on [`ProbabilitySimplex`](@ref) `M` for tangent
+vector `X`. It is computed using isometry with positive orthant of a sphere.
+"""
+function volume_density(M::ProbabilitySimplex{N}, p, X) where {N}
+    pe = simplex_to_amplitude(M, p)
+    return volume_density(Sphere(N), pe, simplex_to_amplitude_diff(M, p, X))
+end
+
 @doc raw"""
     zero_vector(M::ProbabilitySimplex, p)
 
diff --git a/src/manifolds/ProbabilitySimplexEuclideanMetric.jl b/src/manifolds/ProbabilitySimplexEuclideanMetric.jl
index 3665bd7656..791237a6c6 100644
--- a/src/manifolds/ProbabilitySimplexEuclideanMetric.jl
+++ b/src/manifolds/ProbabilitySimplexEuclideanMetric.jl
@@ -9,6 +9,18 @@ function exp!(
     return (q .= p .+ t .* X)
 end
 
+@doc raw"""
+    manifold_volume(::MetricManifold{ℝ,<:ProbabilitySimplex{n},<:EuclideanMetric})) where {n}
+
+Return the volume of the [`ProbabilitySimplex`](@ref) with the Euclidean metric.
+The formula reads ``\frac{\sqrt{n+1}}{n!}``
+"""
+function manifold_volume(
+    ::MetricManifold{ℝ,<:ProbabilitySimplex{n},<:EuclideanMetric},
+) where {n}
+    return sqrt(n + 1) / factorial(n)
+end
+
 @doc raw"""
     rand(::MetricManifold{ℝ,<:ProbabilitySimplex,<:EuclideanMetric}; vector_at=nothing, σ::Real=1.0)
 
@@ -44,3 +56,13 @@ function Random.rand!(
     end
     return pX
 end
+
+@doc raw"""
+    volume_density(::MetricManifold{ℝ,<:ProbabilitySimplex,<:EuclideanMetric}, p, X)
+
+Compute the volume density at point `p` on [`ProbabilitySimplex`](@ref) `M` for tangent
+vector `X`. It is equal to 1.
+"""
+function volume_density(::MetricManifold{ℝ,<:ProbabilitySimplex,<:EuclideanMetric}, p, X)
+    return one(eltype(X))
+end
diff --git a/src/manifolds/ProductManifold.jl b/src/manifolds/ProductManifold.jl
index eec31d32d9..197bf9088e 100644
--- a/src/manifolds/ProductManifold.jl
+++ b/src/manifolds/ProductManifold.jl
@@ -967,6 +967,14 @@ manifold dimensions the product is made of.
 """
 manifold_dimension(M::ProductManifold) = mapreduce(manifold_dimension, +, M.manifolds)
 
+"""
+    manifold_dimension(M::ProductManifold)
+
+Return the volume of [`ProductManifold`](@ref) `M`, i.e. product of volumes of the
+manifolds `M` is constructed from.
+"""
+manifold_volume(M::ProductManifold) = mapreduce(manifold_volume, *, M.manifolds)
+
 function mid_point!(M::ProductManifold, q, p1, p2)
     map(
         mid_point!,
@@ -1639,6 +1647,22 @@ function vector_transport_to!(M::ProductManifold, Y, p, X, q, m::ParallelTranspo
     return Y
 end
 
+@doc raw"""
+    volume_density(M::ProductManifold, p, X)
+
+Return volume density on the [`ProductManifold`](@ref) `M`, i.e. product of constituent
+volume densities.
+"""
+function volume_density(M::ProductManifold, p, X)
+    dens = map(
+        volume_density,
+        M.manifolds,
+        submanifold_components(M, p),
+        submanifold_components(M, X),
+    )
+    return prod(dens)
+end
+
 function zero_vector!(M::ProductManifold, X, p)
     map(
         zero_vector!,
diff --git a/src/manifolds/ProjectiveSpace.jl b/src/manifolds/ProjectiveSpace.jl
index 71d7dbf62f..3a5bc4e6cd 100644
--- a/src/manifolds/ProjectiveSpace.jl
+++ b/src/manifolds/ProjectiveSpace.jl
@@ -344,6 +344,28 @@ function manifold_dimension(M::AbstractProjectiveSpace{𝔽}) where {𝔽}
     return manifold_dimension(get_embedding(M)) - real_dimension(𝔽)
 end
 
+@doc raw"""
+    manifold_volume(M::AbstractProjectiveSpace{ℝ})
+
+Volume of the ``n``-dimensional [`AbstractProjectiveSpace`](@ref) `M`. The formula reads:
+
+````math
+\frac{\pi^{(n+1)/2}}{Γ((n+1)/2)},
+````
+
+where ``Γ`` denotes the [Gamma function](https://en.wikipedia.org/wiki/Gamma_function).
+For details see [^BoyaSudarshanTilma2003].
+
+[^BoyaSudarshanTilma2003]:
+    > L. J. Boya, E. C. G. Sudarshan, and T. Tilma, “Volumes of Compact Manifolds,” Reports
+    > on Mathematical Physics, vol. 52, no. 3, pp. 401–422, Dec. 2003,
+    > doi: [10.1016/S0034-4877(03)80038-1](https://doi.org/10.1016/S0034-4877(03)80038-1)
+"""
+function manifold_volume(M::AbstractProjectiveSpace{ℝ})
+    n = manifold_dimension(M) + 1
+    return pi^(n / 2) / gamma(n / 2)
+end
+
 """
     mean(
         M::AbstractProjectiveSpace,
diff --git a/src/manifolds/Sphere.jl b/src/manifolds/Sphere.jl
index a40e0e5503..7d91c0432b 100644
--- a/src/manifolds/Sphere.jl
+++ b/src/manifolds/Sphere.jl
@@ -363,6 +363,22 @@ dimension of the embedding -1.
 """
 manifold_dimension(M::AbstractSphere) = manifold_dimension(get_embedding(M)) - 1
 
+@doc raw"""
+    manifold_volume(M::AbstractSphere{ℝ})
+
+Volume of the ``n``-dimensional [`Sphere`](@ref) `M`. The formula reads
+
+````math
+\operatorname{Vol}(𝕊^{n}) = \frac{2\pi^{(n+1)/2}}{Γ((n+1)/2)},
+````
+
+where ``Γ`` denotes the [Gamma function](https://en.wikipedia.org/wiki/Gamma_function).
+"""
+function manifold_volume(M::AbstractSphere{ℝ})
+    n = manifold_dimension(M) + 1
+    return 2 * pi^(n / 2) / gamma(n / 2)
+end
+
 """
     mean(
         S::AbstractSphere,
@@ -518,6 +534,24 @@ function riemann_tensor!(M::AbstractSphere{ℝ}, Xresult, p, X, Y, Z)
     return Xresult
 end
 
+@doc raw"""
+    volume_density(M::AbstractSphere{ℝ}, p, X)
+
+Compute volume density function of a sphere, i.e. determinant of the differential of
+exponential map `exp(M, p, X)`. The formula reads ``(\sin(\lVert X\rVert)/\lVert X\rVert)^(n-1)``
+where `n` is the dimension of `M`. It is derived from Eq. (4.1) in [^ChevallierLiLuDunson2022].
+
+[^ChevallierLiLuDunson2022]:
+    > E. Chevallier, D. Li, Y. Lu, and D. B. Dunson, “Exponential-wrapped distributions on
+    > symmetric spaces.” arXiv, Oct. 09, 2022.
+    > doi: [10.48550/arXiv.2009.01983](https://doi.org/10.48550/arXiv.2009.01983).
+"""
+function volume_density(M::AbstractSphere{ℝ}, p, X)
+    Xnorm = norm(X)
+    n = manifold_dimension(M) - 1
+    return usinc(Xnorm)^n
+end
+
 """
     StereographicAtlas()
 
diff --git a/src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl b/src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl
index 6bce8bd389..c6987940bb 100644
--- a/src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl
+++ b/src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl
@@ -368,6 +368,13 @@ function log!(::SymmetricPositiveDefinite{N}, X, p, q) where {N}
     return mul!(X, pUe, Se * transpose(pUe))
 end
 
+"""
+    manifold_volume(::SymmetricPositiveDefinite)
+
+Return volume of the [`SymmetricPositiveDefinite`](@ref) manifold, i.e. infinity.
+"""
+manifold_volume(::SymmetricPositiveDefinite) = Inf
+
 @doc raw"""
     parallel_transport_to(M::SymmetricPositiveDefinite, p, X, q)
     parallel_transport_to(M::MetricManifold{SymmetricPositiveDefinite,AffineInvariantMetric}, p, X, y)
@@ -423,7 +430,8 @@ The formula reads[^Rentmeesters2011] ``R(X,Y)Z=p^{1/2}R(X_I, Y_I)Z_Ip^{1/2}``, w
 [^Rentmeesters2011]:
     > Q. Rentmeesters, “A gradient method for geodesic data fitting on some symmetric
     > Riemannian manifolds,” in 2011 50th IEEE Conference on Decision and Control and
-    > European Control Conference, Dec. 2011, pp. 7141–7146. doi: [10.1109/CDC.2011.6161280](https://doi.org/10.1109/CDC.2011.6161280).
+    > European Control Conference, Dec. 2011, pp. 7141–7146.
+    > doi: [10.1109/CDC.2011.6161280](https://doi.org/10.1109/CDC.2011.6161280).
 """
 riemann_tensor(::SymmetricPositiveDefinite, p, X, Y, Z)
 
@@ -437,3 +445,30 @@ function riemann_tensor!(::SymmetricPositiveDefinite, Xresult, p, X, Y, Z)
     Xresult .= ps * (1 // 4 .* (ZI * Xtmp .- Xtmp * ZI)) * ps
     return Xresult
 end
+
+"""
+    volume_density(::SymmetricPositiveDefinite{n}, p, X) where {n}
+
+Compute the volume density of the [`SymmetricPositiveDefinite`](@ref) manifold at `p`
+in direction `X`. See [^ChevallierKalungaAngulo2017], Section 6.2 for details.
+Note that metric in Manifolds.jl has a different scaling factor than the reference.
+
+[^ChevallierKalungaAngulo2017]:
+    > E. Chevallier, E. Kalunga, and J. Angulo, “Kernel Density Estimation on Spaces of
+    > Gaussian Distributions and Symmetric Positive Definite Matrices,” SIAM J. Imaging Sci.,
+    > vol. 10, no. 1, pp. 191–215, Jan. 2017,
+    > doi: [10.1137/15M1053566](https://doi.org/10.1137/15M1053566).
+"""
+function volume_density(::SymmetricPositiveDefinite{n}, p, X) where {n}
+    eig = eigvals(X)
+    dens = 1.0
+    for i in 1:length(eig)
+        for j in (i + 1):length(eig)
+            absdiff = abs(eig[i] - eig[j])
+            if absdiff > eps(absdiff)
+                dens *= sinh(absdiff) / absdiff
+            end
+        end
+    end
+    return dens
+end
diff --git a/test/groups/general_unitary_groups.jl b/test/groups/general_unitary_groups.jl
index a301f2fb02..18c3463067 100644
--- a/test/groups/general_unitary_groups.jl
+++ b/test/groups/general_unitary_groups.jl
@@ -23,6 +23,57 @@ include("group_utils.jl")
             @test isapprox(On, e, X2, X)
             @test log_lie(On, e) == zeros(n, n)
         end
+
+        @test manifold_volume(Orthogonal(1)) ≈ 2
+        @test manifold_volume(Orthogonal(2)) ≈ 4 * π * sqrt(2)
+        @test manifold_volume(Orthogonal(3)) ≈ 16 * π^2 * sqrt(2)
+        @test manifold_volume(Orthogonal(4)) ≈ 2 * (2 * π)^4 * sqrt(2)
+        @test manifold_volume(Orthogonal(5)) ≈ 8 * (2 * π)^6 / 6 * sqrt(2)
+    end
+
+    @testset "Special Orthogonal Group" begin
+        @test manifold_volume(SpecialOrthogonal(1)) ≈ 1
+        @test manifold_volume(SpecialOrthogonal(2)) ≈ 2 * π * sqrt(2)
+        @test manifold_volume(SpecialOrthogonal(3)) ≈ 8 * π^2 * sqrt(2)
+        @test manifold_volume(SpecialOrthogonal(4)) ≈ (2 * π)^4 * sqrt(2)
+        @test manifold_volume(SpecialOrthogonal(5)) ≈ 4 * (2 * π)^6 / 6 * sqrt(2)
+
+        M = SpecialOrthogonal(2)
+        p = [
+            0.22632098602578843 0.9740527764368391
+            -0.9740527764368391 0.22632098602578843
+        ]
+        X = [0.0 -0.7071067811865475; 0.7071067811865475 0.0]
+        @test volume_density(M, p, X) ≈ 1.0
+
+        M = SpecialOrthogonal(3)
+        p = [
+            -0.5908399013383766 -0.6241917041179139 0.5111681988316876
+            -0.7261666986267721 0.13535732881097293 -0.6740625485388226
+            0.35155388888753836 -0.7694563730631729 -0.5332417398896261
+        ]
+        X = [
+            0.0 -0.30777760628130063 0.5499897386953444
+            0.30777760628130063 0.0 -0.32059980100053004
+            -0.5499897386953444 0.32059980100053004 0.0
+        ]
+        @test volume_density(M, p, X) ≈ 0.8440563052346255
+        @test volume_density(M, p, zero(X)) ≈ 1.0
+
+        M = SpecialOrthogonal(4)
+        p = [
+            -0.09091199873970474 -0.5676546886791307 -0.006808638869334249 0.8182034009599919
+            -0.8001176365300662 0.3161567169523502 -0.4938592872334223 0.12633171594159726
+            -0.5890394255366699 -0.2597679221590146 0.7267279425385695 -0.23962403743004465
+            -0.0676707570677516 -0.7143764493344514 -0.4774129704812182 -0.5071132150619608
+        ]
+        X = [
+            0.0 0.2554704296965055 0.26356215573144676 -0.4070678736115306
+            -0.2554704296965055 0.0 -0.04594199053786204 -0.10586374034761421
+            -0.26356215573144676 0.04594199053786204 0.0 0.43156436122007846
+            0.4070678736115306 0.10586374034761421 -0.43156436122007846 0.0
+        ]
+        @test volume_density(M, p, X) ≈ 0.710713830700454
     end
 
     @testset "Unitary Group" begin
@@ -30,6 +81,11 @@ include("group_utils.jl")
         @test repr(U2) == "Unitary(2)"
         @test injectivity_radius(U2) == π
 
+        @test manifold_volume(Unitary(1)) ≈ 2 * π
+        @test manifold_volume(Unitary(2)) ≈ 4 * π^3
+        @test manifold_volume(Unitary(3)) ≈ sqrt(3) * 2 * π^6
+        @test manifold_volume(Unitary(4)) ≈ sqrt(2) * 8 * π^10 / 12
+
         for n in [1, 2, 3]
             Un = Unitary(n)
             X = zeros(ComplexF64, n, n)
@@ -115,6 +171,11 @@ include("group_utils.jl")
         ) # base point wrong
         e = Identity(MultiplicationOperation())
         @test_throws DomainError is_vector(SU2, e, Xe, true, true) # Xe not skew hermitian
+
+        @test manifold_volume(SpecialUnitary(1)) ≈ 1
+        @test manifold_volume(SpecialUnitary(2)) ≈ 2 * π^2
+        @test manifold_volume(SpecialUnitary(3)) ≈ sqrt(3) * π^5
+        @test manifold_volume(SpecialUnitary(4)) ≈ sqrt(2) * 4 * π^9 / 12
     end
 
     @testset "SO(4) and O(4) exp/log edge cases" begin
diff --git a/test/manifolds/circle.jl b/test/manifolds/circle.jl
index 2a6b548965..0187eb732b 100644
--- a/test/manifolds/circle.jl
+++ b/test/manifolds/circle.jl
@@ -115,6 +115,10 @@ using Manifolds: TFVector, CoTFVector
             2.0,
             ManifoldDiff.βdifferential_shortest_geodesic_startpoint,
         ) === 2.0
+
+        # volume
+        @test manifold_volume(M) ≈ 2 * π
+        @test volume_density(M, 0.0, 2.0) == 1.0
     end
     TEST_STATIC_SIZED && @testset "Real Circle and static sized arrays" begin
         X = @MArray fill(0.0)
@@ -236,6 +240,10 @@ using Manifolds: TFVector, CoTFVector
         @test mean(Mc, angles, [1.0, 1.0, 1.0]) ≈ exp(-π * im / 2)
         @test_throws ErrorException mean(Mc, [-1.0 + 0im, 1.0 + 0im])
         @test_throws ErrorException mean(Mc, [-1.0 + 0im, 1.0 + 0im], [1.0, 1.0])
+
+        # volume
+        @test manifold_volume(Mc) ≈ 2 * π
+        @test volume_density(Mc, 1.0 + 0.0im, 2im) == 1.0
     end
     types = [Complex{Float64}]
     TEST_FLOAT32 && push!(types, Complex{Float32})
diff --git a/test/manifolds/euclidean.jl b/test/manifolds/euclidean.jl
index 98e8fe6049..eb071d4403 100644
--- a/test/manifolds/euclidean.jl
+++ b/test/manifolds/euclidean.jl
@@ -365,4 +365,9 @@ using Manifolds: induced_basis
             ManifoldDiff.βdifferential_shortest_geodesic_startpoint,
         ) === 2.0
     end
+
+    @testset "Volume" begin
+        @test manifold_volume(Euclidean(2)) == Inf
+        @test volume_density(E, p, X) == 1.0
+    end
 end
diff --git a/test/manifolds/hyperbolic.jl b/test/manifolds/hyperbolic.jl
index 638daa683f..734e091ccd 100644
--- a/test/manifolds/hyperbolic.jl
+++ b/test/manifolds/hyperbolic.jl
@@ -339,4 +339,12 @@ include("../utils.jl")
         @test dpX[2][1] == -1.0
         @test dpX[2][2].X == 2 * normalize(X)
     end
+    @testset "Manifold volume" begin
+        M = Hyperbolic(2)
+        @test manifold_volume(M) == Inf
+        p = [1.0, 1.0, sqrt(3)]
+        X = [1.0, 2.0, sqrt(3)]
+        @test volume_density(M, p, X) ≈ 2.980406103535168
+        @test volume_density(M, p, [0.0, 0.0, 0.0]) ≈ 1.0
+    end
 end
diff --git a/test/manifolds/positive_numbers.jl b/test/manifolds/positive_numbers.jl
index f88269f7df..c86215d834 100644
--- a/test/manifolds/positive_numbers.jl
+++ b/test/manifolds/positive_numbers.jl
@@ -77,4 +77,12 @@ include("../utils.jl")
             )
         end
     end
+
+    @testset "Manifold volume" begin
+        M5 = PositiveVectors(3)
+        @test isinf(manifold_volume(M))
+        @test isinf(manifold_volume(M5))
+        @test volume_density(M, 0.5, 2.0) ≈ exp(4.0)
+        @test volume_density(M5, [0.5, 1.0, 2.0], [1.0, -1.0, 2.0]) ≈ exp(2.0)
+    end
 end
diff --git a/test/manifolds/power_manifold.jl b/test/manifolds/power_manifold.jl
index db37291fc5..541f04830e 100644
--- a/test/manifolds/power_manifold.jl
+++ b/test/manifolds/power_manifold.jl
@@ -476,4 +476,11 @@ end
             }
         end
     end
+
+    @testset "Manifold volume" begin
+        @test manifold_volume(Ms1) ≈ manifold_volume(Ms)^5
+        p = repeat([1.0, 0.0, 0.0], 1, 5)
+        X = repeat([0.0, 1.0, 0.0], 1, 5)
+        @test volume_density(Ms1, p, X) ≈ volume_density(Ms, p[:, 1], X[:, 1])^5
+    end
 end
diff --git a/test/manifolds/probability_simplex.jl b/test/manifolds/probability_simplex.jl
index 9cfc96b82f..e34c24127a 100644
--- a/test/manifolds/probability_simplex.jl
+++ b/test/manifolds/probability_simplex.jl
@@ -158,4 +158,11 @@ include("../utils.jl")
         @test riemann_tensor(M, p, X, Y, Z) ≈
               [-0.0034821428571428577, -0.005625, 0.009107142857142857]
     end
+
+    @testset "Volume density" begin
+        @test manifold_volume(M) ≈ pi / 2
+        @test volume_density(M, p, Y) ≈ 0.986956111346216
+        @test manifold_volume(M_euc) ≈ sqrt(3) / 2
+        @test volume_density(M_euc, p, Y) ≈ 1.0
+    end
 end
diff --git a/test/manifolds/product_manifold.jl b/test/manifolds/product_manifold.jl
index c9eb3fddc1..237239250b 100644
--- a/test/manifolds/product_manifold.jl
+++ b/test/manifolds/product_manifold.jl
@@ -755,4 +755,27 @@ using RecursiveArrayTools: ArrayPartition
         @test X2 ==
               ArrayPartition([0.3535533905932738, -0.35355339059327373, 0.0], [1.0, 1.5])
     end
+
+    @testset "Manifold volume" begin
+        MS2 = Sphere(2)
+        MS3 = Sphere(3)
+        PM = ProductManifold(MS2, MS3)
+        @test manifold_volume(PM) ≈ manifold_volume(MS2) * manifold_volume(MS3)
+        p1 = [-0.9171596991960276, 0.39792260844341604, -0.02181017790481868]
+        p2 = [
+            -0.5427653626654726,
+            5.420303965772687e-5,
+            -0.8302022885580579,
+            -0.12716099333369416,
+        ]
+        X1 = [-0.35333565579879633, -0.7896159441709865, 0.45204526334685574]
+        X2 = [
+            -0.33940201562492356,
+            0.8092470417550779,
+            0.18290591742514573,
+            0.2548785571950708,
+        ]
+        @test volume_density(PM, ArrayPartition(p1, p2), ArrayPartition(X1, X2)) ≈
+              volume_density(MS2, p1, X1) * volume_density(MS3, p2, X2)
+    end
 end
diff --git a/test/manifolds/projective_space.jl b/test/manifolds/projective_space.jl
index ad381a7839..412922419a 100644
--- a/test/manifolds/projective_space.jl
+++ b/test/manifolds/projective_space.jl
@@ -324,4 +324,11 @@ include("../utils.jl")
             end
         end
     end
+
+    @testset "Volume" begin
+        @test manifold_volume(ProjectiveSpace(0)) ≈ 1
+        @test manifold_volume(ProjectiveSpace(1)) ≈ π
+        @test manifold_volume(ProjectiveSpace(2)) ≈ 2 * π
+        @test manifold_volume(ProjectiveSpace(3)) ≈ π * π
+    end
 end
diff --git a/test/manifolds/rotations.jl b/test/manifolds/rotations.jl
index 25c0f62264..5a6ca9f38b 100644
--- a/test/manifolds/rotations.jl
+++ b/test/manifolds/rotations.jl
@@ -1,7 +1,7 @@
 include("../utils.jl")
 
 @testset "Rotations" begin
-    M = Manifolds.Rotations(2)
+    M = Rotations(2)
     @test repr(M) == "Rotations(2)"
     @test representation_size(M) == (2, 2)
     @test is_flat(M)
@@ -16,15 +16,14 @@ include("../utils.jl")
     types = [Matrix{Float64}, SMatrix{2,2,Float64,4}]
     TEST_FLOAT32 && push!(types, Matrix{Float32})
     TEST_STATIC_SIZED && push!(types, MMatrix{2,2,Float64,4})
-    retraction_methods = [Manifolds.PolarRetraction(), Manifolds.QRRetraction()]
+    retraction_methods = [PolarRetraction(), QRRetraction()]
 
-    inverse_retraction_methods =
-        [Manifolds.PolarInverseRetraction(), Manifolds.QRInverseRetraction()]
+    inverse_retraction_methods = [PolarInverseRetraction(), QRInverseRetraction()]
 
     basis_types = (DefaultOrthonormalBasis(), ProjectedOrthonormalBasis(:svd))
 
     @testset "vee/hat" begin
-        M = Manifolds.Rotations(2)
+        M = Rotations(2)
         Xf = [1.23]
         p = Matrix{Float64}(I, 2, 2)
         X = Manifolds.hat(M, p, Xf)
@@ -80,12 +79,12 @@ include("../utils.jl")
         exp!(M, q, q, X)
         @test norm(q - q2) ≈ 0
 
-        X14_polar = inverse_retract(M, pts[1], pts[4], Manifolds.PolarInverseRetraction())
-        p4_polar = retract(M, pts[1], X14_polar, Manifolds.PolarRetraction())
+        X14_polar = inverse_retract(M, pts[1], pts[4], PolarInverseRetraction())
+        p4_polar = retract(M, pts[1], X14_polar, PolarRetraction())
         @test isapprox(M, pts[4], p4_polar)
 
-        X14_qr = inverse_retract(M, pts[1], pts[4], Manifolds.QRInverseRetraction())
-        p4_qr = retract(M, pts[1], X14_qr, Manifolds.QRRetraction())
+        X14_qr = inverse_retract(M, pts[1], pts[4], QRInverseRetraction())
+        p4_qr = retract(M, pts[1], X14_qr, QRRetraction())
         @test isapprox(M, pts[4], p4_qr)
     end
 
@@ -107,7 +106,7 @@ include("../utils.jl")
     Random.seed!(42)
     for n in (3, 4, 5)
         @testset "Rotations: SO($n)" begin
-            SOn = Manifolds.Rotations(n)
+            SOn = Rotations(n)
             @test !is_flat(SOn)
             ptd = Manifolds.normal_rotation_distribution(SOn, Matrix(1.0I, n, n), 1.0)
             tvd = Manifolds.normal_tvector_distribution(SOn, Matrix(1.0I, n, n), 1.0)
@@ -168,7 +167,7 @@ include("../utils.jl")
         end
     end
     @testset "Test AbstractManifold Point and Tangent Vector checks" begin
-        M = Manifolds.Rotations(2)
+        M = Rotations(2)
         for p in [1, [2.0 0.0; 0.0 1.0], [1.0 0.5; 0.0 1.0]]
             @test_throws DomainError is_point(M, p, true)
             @test !is_point(M, p)
@@ -185,12 +184,12 @@ include("../utils.jl")
         @test is_vector(M, p, X, true)
     end
     @testset "Project point" begin
-        M = Manifolds.Rotations(2)
+        M = Rotations(2)
         p = Matrix{Float64}(I, 2, 2)
         p1 = project(M, p)
         @test is_point(M, p1, true)
 
-        M = Manifolds.Rotations(3)
+        M = Rotations(3)
         p = collect(reshape(1.0:9.0, (3, 3)))
         p2 = project(M, p)
         @test is_point(M, p2, true)
@@ -200,7 +199,7 @@ include("../utils.jl")
         @test is_point(M, x3, true)
     end
     @testset "Convert from Lie algebra representation of tangents to Riemannian submanifold representation" begin
-        M = Manifolds.Rotations(3)
+        M = Rotations(3)
         p = project(M, collect(reshape(1.0:9.0, (3, 3))))
         x = [[0, -1, 3] [1, 0, 2] [-3, -2, 0]]
         @test is_vector(M, p, x, true)
@@ -236,4 +235,33 @@ include("../utils.jl")
         @test injectivity_radius(M, PolarRetraction()) == 0.0
         @test injectivity_radius(M, p, PolarRetraction()) == 0.0
     end
+
+    @testset "riemann_tensor" begin
+        M = Rotations(3)
+        p = [
+            -0.5908399013383766 -0.6241917041179139 0.5111681988316876
+            -0.7261666986267721 0.13535732881097293 -0.6740625485388226
+            0.35155388888753836 -0.7694563730631729 -0.5332417398896261
+        ]
+        X = [
+            0.0 -0.30777760628130063 0.5499897386953444
+            0.30777760628130063 0.0 -0.32059980100053004
+            -0.5499897386953444 0.32059980100053004 0.0
+        ]
+        Y = [
+            0.0 -0.4821890003925358 -0.3513148535122392
+            0.4821890003925358 0.0 0.37956770358148356
+            0.3513148535122392 -0.37956770358148356 0.0
+        ]
+        Z = [
+            0.0 0.3980141785048982 0.09735377380829331
+            -0.3980141785048982 0.0 -0.576287216962475
+            -0.09735377380829331 0.576287216962475 0.0
+        ]
+        @test riemann_tensor(M, p, X, Y, Z) ≈ [
+            0.0 0.04818900625787811 -0.050996416671166
+            -0.04818900625787811 0.0 0.024666891276861697
+            0.050996416671166 -0.024666891276861697 0.0
+        ]
+    end
 end
diff --git a/test/manifolds/sphere.jl b/test/manifolds/sphere.jl
index 41fac3e55d..61cbe0d013 100644
--- a/test/manifolds/sphere.jl
+++ b/test/manifolds/sphere.jl
@@ -264,4 +264,14 @@ using ManifoldsBase: TFVector
         @test dpX[2][1] == 1.0
         @test dpX[2][2].X == normalize(X)
     end
+
+    @testset "Volume density" begin
+        M = Sphere(2)
+        p = [1.0, 0.0, 0.0]
+        @test manifold_volume(Sphere(0)) ≈ 2
+        @test manifold_volume(Sphere(1)) ≈ 2 * π
+        @test manifold_volume(Sphere(2)) ≈ 4 * π
+        @test manifold_volume(Sphere(3)) ≈ 2 * π * π
+        @test volume_density(M, p, [0.0, 0.5, 0.5]) ≈ 0.9187253698655684
+    end
 end
diff --git a/test/manifolds/symmetric_positive_definite.jl b/test/manifolds/symmetric_positive_definite.jl
index 75c9694d16..27cb3f7d3d 100644
--- a/test/manifolds/symmetric_positive_definite.jl
+++ b/test/manifolds/symmetric_positive_definite.jl
@@ -276,4 +276,20 @@ include("../utils.jl")
         ]
         @test is_point(M, p1)
     end
+
+    @testset "Volume density" begin
+        M = SymmetricPositiveDefinite(3)
+        @test manifold_volume(M) == Inf
+        p = [
+            1.680908185710701 -0.030208936760309613 -0.284402826783584
+            -0.030208936760309613 1.7199740873691465 -0.0066638025832747305
+            -0.284402826783584 -0.0066638025832747305 1.6842441059901379
+        ]
+        X = [
+            -1.3447374559982306 0.2591587910302816 0.9739140395169474
+            0.2591587910302816 0.3649975914025044 -0.2888865063584093
+            0.9739140395169474 -0.2888865063584093 -0.9564259306801289
+        ]
+        @test volume_density(M, p, X) ≈ 5.141867280770719
+    end
 end
diff --git a/tutorials/Project.toml b/tutorials/Project.toml
index 1d0d4c7d2f..ffcf36e7be 100644
--- a/tutorials/Project.toml
+++ b/tutorials/Project.toml
@@ -1,4 +1,5 @@
 [deps]
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
@@ -11,10 +12,12 @@ MultivariateStats = "6f286f6a-111f-5878-ab1e-185364afe411"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
 CSV = "0.10"
 DataFrames = "1"
+Distributions = "0.22.6, 0.23, 0.24, 0.25"
 IJulia = "1"
 Manifolds = "0.8.46"
 ManifoldsBase = "0.14.5"
diff --git a/tutorials/integration.qmd b/tutorials/integration.qmd
new file mode 100644
index 0000000000..4f51f6631d
--- /dev/null
+++ b/tutorials/integration.qmd
@@ -0,0 +1,198 @@
+---
+title: Integration
+---
+
+```{julia}
+#| echo: false
+#| code-fold: true
+#| output: false
+using Pkg;
+cd(@__DIR__)
+Pkg.activate("."); # for reproducibility use the local tutorial environment.
+Pkg.develop(PackageSpec(; path=(@__DIR__) * "/../"))
+using Markdown
+```
+
+This part of documentation covers integration of scalar functions defined on manifolds $f \colon \mathcal{M} \to \mathbb{R}$:
+
+```math
+\int_{\mathcal M} f(p) \mathrm{d}p
+```
+
+The basic concepts are derived from geometric measure theory.
+In principle, there are many ways in which a manifold can be equipped with a measure that can be later used to define an integral.
+One of the most popular ways is based on pushing the Lebesgue measure on a tangent space through the exponential map.
+Any other suitable atlas could be used, not just the one defined by normal coordinates, though each one requires different volume density corrections due to the Jacobian determinant of the pushforward.
+`Manifolds.jl` provides the function [`volume_density`](@ref) that calculates that quantity, denoted $\theta_p(X)$.
+See for example [^leBrigantPuechmorel2019], Definition 11, for a precise description using Jacobi fields.
+
+While many sources define volume density as a function of two points, `Manifolds.jl` decided to use the more general point-tangent vector formulation. The two-points variant can be implemented as
+```{julia}
+using Manifolds
+volume_density_two_points(M::AbstractManifold, p, q) = volume_density(M, p, log(M, p, q))
+```
+
+The simplest way to of integrating a function on a compact manifold is through a [📖 Monte Carlo integrator](https://en.wikipedia.org/wiki/Monte_Carlo_integration).
+A simple variant can be implemented as follows (assuming uniform distribution of `rand`):
+
+```{julia}
+using LinearAlgebra, Distributions, SpecialFunctions
+function simple_mc_integrate(M::AbstractManifold, f; N::Int = 1000)
+    V = manifold_volume(M)
+    sum = 0.0
+    q = rand(M)
+    for i in 1:N
+        sum += f(M, q)
+        rand!(M, q)
+    end
+    return V * sum/N
+end
+```
+
+We used the function [`manifold_volume`](@ref) to get the volume of the set over which the integration is performed, as described in the linked Wikipedia article.
+
+## Distributions
+
+We will now try to verify that volume density correction correctly changes probability density of an exponential-wrapped normal distribution.
+`pdf_tangent_space` (defined in the next code block) represents probability density of a normally distributed random variable $X_T$ in the tangent space $T_p \mathcal{M}$.
+Its probability density (with respect to the Lebesgue measure of the tangent space) is $f_{X_T}\colon T_p \mathcal{M} \to \mathbb{R}$.
+
+`pdf_manifold` (defined below) refers to the probability density of the distribution $X_M$ from the tangent space $T_p \mathcal{M}$ wrapped using exponential map on the manifold.
+The formula for probability density with respect to pushforward measure of the Lebesgue measure in the tangent space reads
+```math
+f_{X_M}(q) = \sum_{X \in T_p\mathcal{M}, \exp_p(X)=q} \frac{f_{X_T}(X)}{\theta_p(X)}
+```
+[`volume_density`](@ref) function calculates the correction $\theta_p(X)$.
+
+
+```{julia}
+function pdf_tangent_space(M::AbstractManifold, p)
+    return pdf(MvNormal(zeros(manifold_dimension(M)), 0.2*I), p)
+end
+
+function pdf_manifold(M::AbstractManifold, q)
+    p = [1.0, 0.0, 0.0]
+    X = log(M, p, q)
+    Xc = get_coordinates(M, p, X, DefaultOrthonormalBasis())
+    vd = abs(volume_density(M, p, X))
+    if vd > eps()
+        return pdf_tangent_space(M, Xc) / vd
+    else
+        return 0.0
+    end
+end
+
+println(simple_mc_integrate(Sphere(2), pdf_manifold; N=1000000))
+```
+
+The function `simple_mc_integrate`, defined in the previous section, is used to verify that the density integrates to 1 over the manifold.
+
+Note that our `pdf_manifold` implements a simplified version of $f_{X_M}$ which assumes that the probability mass of `pdf_tangent_space` outside of (local) injectivity radius at $p$ is negligible.
+In such case there is only one non-zero summand in the formula for $f_{X_M}(q)$, namely $X=\log_p(q)$.
+Otherwise we would have to consider other vectors $Y\in T_p \mathcal{M}$ such that $\exp_p(Y) = q$ in that sum.
+
+Remarkably, exponential-wrapped distributions possess three important qualities [^ChevallierLiLuDunson2022]:
+
+* Densities of $X_M$ are explicit. There is no normalization constant that needs to be computed like in truncated distributions.
+* Sampling from $X_M$ is easy. It suffices to get a sample from $X_T$ and pass it to the exponential map.
+* If mean of $X_T$ is 0, then there is a simple correspondence between moments of $X_M$ and $X_T$, for example $p$ is the mean of $X_M$.
+
+## Kernel density estimation
+
+We can also make a Pelletier's isotropic kernel density estimator. Given points $p_1, p_2, \dots, p_n$ on $d$-dimensional manifold $\mathcal M$
+the density at point $q$ is defined as
+```math
+f(q) = \frac{1}{n h^d} \sum_{i=1}^n \frac{1}{\theta_q(\log_q(p_i))}K\left( \frac{d(q, p_i)}{h} \right),
+```
+
+where $h$ is the bandwidth, a small positive number less than the injectivity radius of $\mathcal M$ and $K\colon\mathbb{R}\to\mathbb{R}$ is a kernel function.
+Note that Pelletier's estimator can only use radially-symmetric kernels.
+The radially symmetric multivariate Epanechnikov kernel used in the example below is described in [^LangrenéWarin2019].
+
+```{julia}
+struct PelletierKDE{TM<:AbstractManifold,TPts<:AbstractVector}
+    M::TM
+    bandwidth::Float64
+    pts::TPts
+end
+
+(kde::PelletierKDE)(::AbstractManifold, p) = kde(p)
+function (kde::PelletierKDE)(p)
+    n = length(kde.pts)
+    d = manifold_dimension(kde.M)
+    sum_kde = 0.0
+    function epanechnikov_kernel(x)
+        if x < 1
+            return gamma(2+d/2) * (1-x^2)/(π^(d/2))
+        else
+            return 0.0
+        end
+    end
+    for i in 1:n
+        X = log(kde.M, p, kde.pts[i])
+        Xn = norm(kde.M, p, X)
+        sum_kde += epanechnikov_kernel(Xn / kde.bandwidth) / volume_density(kde.M, p, X)
+    end
+    sum_kde /= n * kde.bandwidth^d
+    return sum_kde
+end
+
+M = Sphere(2)
+pts = rand(M, 8)
+kde = PelletierKDE(M, 0.7, pts)
+println(simple_mc_integrate(Sphere(2), kde; N=1000000))
+println(kde(rand(M)))
+```
+
+
+## Technical notes
+
+This section contains a few technical notes that are relevant to the problem of integration on manifolds but can be freely skipped on the first read of the tutorial.
+
+### Conflicting statements about volume of a manifold
+
+[`manifold_volume`](@ref) and [`volume_density`](@ref) are closely related to each other, though very few sources explore this connection, and some even claiming a certain level of arbitrariness in defining `manifold_volume`.
+Volume is sometimes considered arbitrary because Riemannian metrics on some spaces like the manifold of rotations are defined with arbitrary constants.
+However, once a constant is picked (and it must be picked before any useful computation can be performed), all geometric operations must follow in a consistent way: inner products, exponential and logarithmic maps, volume densities, etc.
+`Manifolds.jl` consistently picks such constants and provides a unified framework, though it sometimes results in picking a different constant than what is the most popular in some sub-communities.
+
+### Haar measures
+
+On Lie groups the situation regarding integration is more complicated.
+Invariance under left or right group action is a desired property that leads one to consider Haar measures [^Tornier2020].
+It is, however, unclear what are the practical benefits of considering Haar measures over the Lebesgue measure of the underlying manifold, which often turns out to be invariant anyway.
+
+### Integration in charts
+
+Integration through charts is an approach currently not supported by `Manifolds.jl`.
+One has to define a suitable set of disjoint charts covering the entire manifold and use a method for multivariate Euclidean integration.
+Note that ranges of parameters have to be adjusted for each manifold and scaling based on the metric needs to be applied.
+See [^BoyaSudarshanTilma2003] for some considerations on symmetric spaces.
+
+
+## References
+
+[^Tornier2020]:
+    > S. Tornier, “Haar Measures.” arXiv, Jun. 19, 2020.
+    > doi: [10.48550/arXiv.2006.10956](https://doi.org/10.48550/arXiv.2006.10956).
+
+[^LangrenéWarin2019]:
+    > N. Langrené and X. Warin, “Fast and Stable Multivariate Kernel Density Estimation by
+    > Fast Sum Updating,” Journal of Computational and Graphical Statistics, vol. 28, no. 3,
+    > pp. 596–608, Jul. 2019,
+    > doi: [10.1080/10618600.2018.1549052](https://doi.org/10.1080/10618600.2018.1549052).
+
+[^leBrigantPuechmorel2019]:
+    > A. le Brigant and S. Puechmorel, “Approximation of Densities on Riemannian Manifolds,”
+    > Entropy, vol. 21, no. 1, Art. no. 1, Jan. 2019,
+    > doi: [10.3390/e21010043](https://doi.org/10.3390/e21010043).
+
+[^ChevallierLiLuDunson2022]:
+    > E. Chevallier, D. Li, Y. Lu, and D. B. Dunson, “Exponential-wrapped distributions on
+    > symmetric spaces.” arXiv, Oct. 09, 2022.
+    > doi: [10.48550/arXiv.2009.01983](https://doi.org/10.48550/arXiv.2009.01983).
+
+[^BoyaSudarshanTilma2003]:
+    > L. J. Boya, E. C. G. Sudarshan, and T. Tilma, “Volumes of Compact Manifolds,” Reports
+    > on Mathematical Physics, vol. 52, no. 3, pp. 401–422, Dec. 2003,
+    > doi: [10.1016/S0034-4877(03)80038-1](https://doi.org/10.1016/S0034-4877(03)80038-1)