diff --git a/VirtualRobot/CMakeLists.txt b/VirtualRobot/CMakeLists.txt
index 1bf2d54b71163e15b6843bcb3767387584f9a407..7af2f9508daa92b4f2dc8b40608c1b9a09a01cf6 100644
--- a/VirtualRobot/CMakeLists.txt
+++ b/VirtualRobot/CMakeLists.txt
@@ -334,6 +334,7 @@ math/Plane.cpp
 math/Primitive.cpp
 math/Triangle.cpp
 math/WeightedAverage.cpp
+math/statistics/measures.cpp
 Util/xml/tinyxml2.cpp
 Util/json/eigen_conversion.cpp
 )
@@ -507,6 +508,7 @@ math/TransformedFunctionR1R3.h
 math/TransformedFunctionR2R3.h
 math/Triangle.h
 math/WeightedAverage.h
+math/statistics/measures.h
 Util/xml/tinyxml2.h
 Util/json/json.hpp
 Util/json/eigen_conversion.h
diff --git a/VirtualRobot/math/statistics/measures.cpp b/VirtualRobot/math/statistics/measures.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..59fa4486d9161e488cae2735b8a40c9da9a916e1
--- /dev/null
+++ b/VirtualRobot/math/statistics/measures.cpp
@@ -0,0 +1,110 @@
+#include "measures.h"
+
+#include <algorithm>
+#include <cmath>
+#include <stdexcept>
+
+
+namespace math { namespace stat
+{
+
+void sort(std::vector<float>& values)
+{
+    std::sort(values.begin(), values.end());
+}
+
+std::vector<float> sorted(const std::vector<float>& values)
+{
+    std::vector<float> s = values;
+    std::sort(s.begin(), s.end());
+    return s;
+}
+
+float min(const std::vector<float>& values, bool isSorted)
+{
+    return isSorted ? values.front() : *std::min_element(values.begin(), values.end());
+}
+
+float max(const std::vector<float>& values, bool isSorted)
+{
+    return isSorted ? values.back() : *std::max_element(values.begin(), values.end());
+}
+
+float mean(const std::vector<float>& values)
+{
+    if (values.empty() == 0)
+        throw std::range_error("Given values are empty.");
+    
+    float sum = 0;
+    for (float v : values)
+    {
+        sum += v;
+    }
+    return sum / values.size();
+}
+
+float stddev(const std::vector<float>& values)
+{
+    return stddev(values, mean(values));
+}
+
+float stddev(const std::vector<float>& values, float mean)
+{
+    float sum = 0;
+    for (float v : values)
+    {
+        float diff = v - mean;
+        sum += diff * diff;
+    }
+    float variance = sum / (values.size() - 1);
+    return std::sqrt(variance);
+}
+
+float quantile(const std::vector<float>& _values, float p, bool isSorted)
+{
+    const std::vector<float>& values = isSorted ? _values : sorted(_values);
+    
+    float location = p * values.size();
+    
+    std::size_t floor = static_cast<std::size_t>(std::floor(location));
+    std::size_t ceil = static_cast<std::size_t>(std::ceil(location));
+    
+    if (floor == ceil)
+    {
+        return values.at(floor);
+    }
+    else
+    {
+        float t = location - floor;
+        return (1 - t) * values.at(floor) + t * values.at(ceil);
+    }
+}
+
+float lowerQuartile(const std::vector<float>& values, bool isSorted)
+{
+    return quantile(values, .25, isSorted);
+}
+
+float median(const std::vector<float>& values, bool isSorted)
+{
+    return quantile(values, .5, isSorted);
+}
+
+float upperQuartile(const std::vector<float>& values, bool isSorted)
+{
+    return quantile(values, .75, isSorted);
+}
+
+float interquartileRange(const std::vector<float>& _values, bool isSorted)
+{
+    const std::vector<float>& values = isSorted ? _values : sorted(_values);
+    return interquartileRange(lowerQuartile(values, true), upperQuartile(values, true));
+}
+
+float interquartileRange(float lowerQuartile, float upperQuartile)
+{
+    return upperQuartile - lowerQuartile;
+}
+
+
+}}
diff --git a/VirtualRobot/math/statistics/measures.h b/VirtualRobot/math/statistics/measures.h
new file mode 100644
index 0000000000000000000000000000000000000000..dd18caa2b7a4d1497eac75ea82dd6d43004cb2b2
--- /dev/null
+++ b/VirtualRobot/math/statistics/measures.h
@@ -0,0 +1,31 @@
+#pragma once
+
+#include <vector>
+
+
+namespace math { namespace stat
+{
+
+    void sort(std::vector<float>& values);
+    std::vector<float> sorted(const std::vector<float>& values);
+
+    
+    float min(const std::vector<float>& values, bool isSorted=false);
+    float max(const std::vector<float>& values, bool isSorted=false);
+    float mean(const std::vector<float>& values);
+    
+    float stddev(const std::vector<float>& values);
+    float stddev(const std::vector<float>& values, float mean);
+    
+    float quantile(const std::vector<float>& values, float p, bool isSorted=false);
+    
+    float lowerQuartile(const std::vector<float>& values, bool isSorted=false);
+    float median(const std::vector<float>& values, bool isSorted=false);
+    float upperQuartile(const std::vector<float>& values, bool isSorted=false);
+    
+    float interquartileRange(const std::vector<float>& values, bool isSorted=false);
+    float interquartileRange(float lowerQuartile, float upperQuartile);
+    
+
+
+}}