diff --git a/Project.toml b/Project.toml
index e5107d0..4ddda94 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,7 +7,6 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 Quaternions = "0.5.3, 0.6, 0.7"
diff --git a/src/Rotations.jl b/src/Rotations.jl
index 07229ab..965a9b4 100644
--- a/src/Rotations.jl
+++ b/src/Rotations.jl
@@ -8,8 +8,6 @@ using StaticArrays
 using Random
 using Quaternions
 
-import Statistics
-
 include("util.jl")
 include("core_types.jl")
 include("unitquaternion.jl")
diff --git a/src/mean.jl b/src/mean.jl
deleted file mode 100644
index 55e6bc1..0000000
--- a/src/mean.jl
+++ /dev/null
@@ -1,51 +0,0 @@
-#############################################################
-#
-#  add extra quaternion operations for Quaternions
-#  some should probably be merged back into Quaternions
-#
-#############################################################
-
-#=
-"""
-    mean(qvec, [method = 0])
-
-Function to find the mean of a set of rotations.
-
-Args:
-
-    qvec    - Vector of `Rotation`s
-
-    method  - if 0 (default), the returned rotation minimizes sum_i=1:n (||R(Qi) - R(Qbar)||)  with || the Frobenius norm, or
-              equivalently minimizes sum_i=1:n (sin(phi_i / 2))^2, where phi_i is the angle of Qe_i, with  Qbar = Q_i x Qe_i
-
-            - if 1, the returned rotation minimizes sum_i=1:n (phi_i ^2) where phi_i is the angle of Qe_i, with  Qbar = Q_i x Qe_i
-              this should be about the same as method == 0 if phi_i are all small (note: still not implemented)
-"""
-=#
-
-"""
-    mean(r::AbstractVector{Rotation})  ->  r̄
-
-Find the mean `r̄` of a set of rotations `r`. The returned rotation minimizes `sum_i=1:n (||r[i] - r̄||)`  with `||` the Frobenius norm, or
-equivalently minimizes `sum_i=1:n (sin(ϕ[i] / 2))^2`, where `ϕ[i] = rotation_angle(r[i] / r̄)`.
-"""
-function Statistics.mean(qvec::AbstractVector{QuatRotation{T}}, method::Integer = 0) where T
-    #if (method == 0)
-        M = zeros(4, 4)
-        for i in 1:length(qvec)
-            q = qvec[i]
-            Qi = @SVector [q.w, q.x, q.y, q.z]  # convert types to ensure we don't get quaternion multiplication
-            M .+= Qi * (Qi')
-        end
-        evec = eigfact(Symmetric(M), 4:4)
-        Qbar = QuatRotation(evec.vectors[1], evec.vectors[2], evec.vectors[3], evec.vectors[4]) # This will renormalize the quaternion...
-    #else
-    #    error("I haven't coded this")
-    #end
-
-    return Qbar
-end
-
-function Statistics.mean(vec::AbstractVector{R}) where R<:Rotation
-    R(mean(convert(Vector{QuatRotation{eltype(R)}}, vec)))
-end