diff --git a/Project.toml b/Project.toml
index 8ddfcc5a..5fc04a8d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -2,6 +2,9 @@ name = "FixedPointNumbers"
 uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
 version = "0.8.1"
 
+[deps]
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
 [compat]
 julia = "1"
 
diff --git a/src/FixedPointNumbers.jl b/src/FixedPointNumbers.jl
index 3c3bd338..81be767f 100644
--- a/src/FixedPointNumbers.jl
+++ b/src/FixedPointNumbers.jl
@@ -8,6 +8,8 @@ import Base: ==, <, <=, -, +, *, /, ~, isapprox,
              div, fld, rem, mod, mod1, fld1, min, max, minmax,
              rand, length
 
+import Statistics   # for _mean_promote
+
 using Base.Checked: checked_add, checked_sub, checked_div
 
 using Base: @pure
@@ -225,6 +227,9 @@ Base.mul_prod(x::FixedPoint, y::FixedPoint) = Treduce(x) * Treduce(y)
 Base.reduce_empty(::typeof(Base.mul_prod), ::Type{F}) where {F<:FixedPoint} = one(Treduce)
 Base.reduce_first(::typeof(Base.mul_prod), x::FixedPoint)  = Treduce(x)
 
+if isdefined(Statistics, :_mean_promote)
+    Statistics._mean_promote(x::Real, y::FixedPoint) = Treduce(y)
+end
 
 """
     sd, ad = scaledual(s::Number, a)
diff --git a/test/fixed.jl b/test/fixed.jl
index dd0764a8..0feec22f 100644
--- a/test/fixed.jl
+++ b/test/fixed.jl
@@ -1,4 +1,4 @@
-using FixedPointNumbers, Test
+using FixedPointNumbers, Statistics, Test
 using FixedPointNumbers: bitwidth
 
 function test_op(fun::F, ::Type{T}, fx, fy, fxf, fyf, tol) where {F,T}
@@ -242,6 +242,17 @@ end
     @test prod(a, dims=1) == [acmp]
 end
 
+@testset "reductions, Statistics" begin
+    a = Q1f6[0.75, 0.5]
+    af = FixedPointNumbers.Treduce.(a)
+    @test mean(a) === mean(af)
+    @test std(a)  === std(af)
+    @test var(a)  === var(af)
+    m = mean(a)
+    @test stdm(a, m) === stdm(af, m)
+    @test varm(a, m) === varm(af, m)
+end
+
 @testset "convert result type" begin
     x = Fixed{Int8,7}(0.75)
     for T in (Float16, Float32, Float64, BigFloat)
diff --git a/test/normed.jl b/test/normed.jl
index 56831d55..aa1cbf6e 100644
--- a/test/normed.jl
+++ b/test/normed.jl
@@ -1,4 +1,4 @@
-using FixedPointNumbers, Test
+using FixedPointNumbers, Statistics, Test
 using FixedPointNumbers: bitwidth
 
 @testset "domain of f" begin
@@ -449,6 +449,17 @@ end
     @test prod(a, dims=1) == [acmp]
 end
 
+@testset "reductions, Statistics" begin
+    a = N0f8[reinterpret(N0f8, 0x80), reinterpret(N0f8, 0x40)]
+    af = FixedPointNumbers.Treduce.(a)
+    @test mean(a) === mean(af)
+    @test std(a)  === std(af)
+    @test var(a)  === var(af)
+    m = mean(a)
+    @test stdm(a, m) === stdm(af, m)
+    @test varm(a, m) === varm(af, m)
+end
+
 @testset "rand" begin
     for T in (Normed{UInt8,8}, Normed{UInt8,6},
               Normed{UInt16,16}, Normed{UInt16,14},