diff --git a/applications/test/Distribution/DistributionTest.C b/applications/test/Distribution/DistributionTest.C
index cdbbc4e469..687d5746ac 100644
--- a/applications/test/Distribution/DistributionTest.C
+++ b/applications/test/Distribution/DistributionTest.C
@@ -28,12 +28,24 @@ Application
 Description
     Test the Distribution class
 
+    Plot normal distribution test in gnuplot using:
+
+    @verbatim
+    normalDistribution(mean, sigma, x) = \
+        sqrt(1.0/(2.0*pi*sigma**2))*exp(-(x - mean)**2.0/(2.0*sigma**2))
+
+    plot normalDistribution(8.5, 2.5, x)
+    @endverbatim
+
 \*---------------------------------------------------------------------------*/
 
 #include "vector.H"
 #include "labelVector.H"
 #include "tensor.H"
 #include "Distribution.H"
+#include "Random.H"
+#include "dimensionedTypes.H"
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -41,36 +53,173 @@ using namespace Foam;
 
 int main(int argc, char *argv[])
 {
-    Distribution<scalar> dS(scalar(3.0));
-    Distribution<vector> dV(vector(0.1, 0.24, 0.18));
-    Distribution<labelVector> dLV(labelVector(2,3,4));
+    Random R(918273);
 
-    dS.add(1.0);
-    dS.add(2.0);
-    dS.add(12.0);
-    dS.add(1.3);
+    {
+        // scalar
+        label randomDistributionTestSize = 50000000;
 
-    Info<< dS << nl << dS.raw() << endl;
+        Distribution<scalar> dS(scalar(5e-2));
 
-    vector vA(1.2, 1.3, 1.1);
-    vector vB(1.3, 1.5, 1.6);
-    vector vC(0.5, 5.3, 1.1);
+        Info<< nl << "Distribution<scalar>" << nl
+            << "Sampling "
+            << randomDistributionTestSize
+            << " times from GaussNormal distribution."
+            << endl;
 
-    dV.add(vA);
-    dV.add(vB);
-    dV.add(vC);
+        for (label i = 0; i < randomDistributionTestSize; i++)
+        {
+            dS.add(2.5*R.GaussNormal() + 8.5);
+        }
 
-    Info<< dV << nl << dV.raw() << endl;
+        Info<< "Mean " << dS.mean() << nl
+            << "Median " << dS.median()
+            << endl;
 
-    labelVector lVA(6, 8, 11);
-    labelVector lVB(-12, -3, 6);
-    labelVector lVC(-4, -2, 5);
+        dS.write("Distribution_scalar_test", dS.normalised());
+    }
 
-    dLV.add(lVA);
-    dLV.add(lVB);
-    dLV.add(lVC);
+    {
+        // vector
+        Distribution<vector> dV(vector(0.1, 0.05, 0.15));
 
-    Info<< dLV << nl << dLV.raw() << endl;
+        label randomDistributionTestSize = 1000000;
+
+        Info<< nl << "Distribution<vector>" << nl
+            << "Sampling "
+            << randomDistributionTestSize
+            << " times from uniform and GaussNormal distribution."
+            << endl;
+
+        for (label i = 0; i < randomDistributionTestSize; i++)
+        {
+            dV.add(R.vector01());
+
+            // Adding separate GaussNormal components with component
+            // weights
+
+            dV.add
+            (
+                vector
+                (
+                    R.GaussNormal()*3.0 + 1.5,
+                    R.GaussNormal()*0.25 + 4.0,
+                    R.GaussNormal()*3.0 - 1.5
+                ),
+                vector(1.0, 2.0, 5.0)
+            );
+        }
+
+        Info<< "Mean " << dV.mean() << nl
+            << "Median " << dV.median()
+            << endl;
+
+        dV.write("Distribution_vector_test", dV.normalised());
+    }
+
+    {
+        // labelVector
+        Distribution<labelVector> dLV(labelVector::one*10);
+
+        label randomDistributionTestSize = 2000000;
+
+        Info<< nl << "Distribution<labelVector>" << nl
+            << "Sampling "
+            << randomDistributionTestSize
+            << " times from uniform distribution."
+            << endl;
+
+        for (label i = 0; i < randomDistributionTestSize; i++)
+        {
+            dLV.add
+            (
+                labelVector
+                (
+                    R.integer(-1000, 1000),
+                    R.integer(-5000, 5000),
+                    R.integer(-2000, 7000)
+                )
+            );
+        }
+
+        Info<< "Mean " << dLV.mean() << nl
+            << "Median " << dLV.median()
+            << endl;
+
+        dLV.write("Distribution_labelVector_test", dLV.normalised());
+    }
+
+    {
+        // tensor
+        Distribution<tensor> dT(tensor::one*1e-2);
+
+        label randomDistributionTestSize = 2000000;
+
+        Info<< nl << "Distribution<tensor>" << nl
+            << "Sampling "
+            << randomDistributionTestSize
+            << " times from uniform distribution."
+            << endl;
+
+        for (label i = 0; i < randomDistributionTestSize; i++)
+        {
+            dT.add(R.tensor01());
+        }
+
+        Info<< "Mean " << dT.mean() << nl
+            << "Median " << dT.median()
+            << endl;
+
+        dT.write("Distribution_tensor_test", dT.normalised());
+    }
+
+    {
+        // symmTensor
+        Distribution<symmTensor> dSyT(symmTensor::one*1e-2);
+
+        label randomDistributionTestSize = 2000000;
+
+        Info<< nl << "Distribution<symmTensor>" << nl
+            << "Sampling "
+            << randomDistributionTestSize
+            << " times from uniform distribution."
+            << endl;
+
+        for (label i = 0; i < randomDistributionTestSize; i++)
+        {
+            dSyT.add(R.symmTensor01());
+        }
+
+        Info<< "Mean " << dSyT.mean() << nl
+            << "Median " << dSyT.median()
+            << endl;
+
+        dSyT.write("Distribution_symmTensor_test", dSyT.normalised());
+    }
+
+    {
+        // sphericalTensor
+        Distribution<sphericalTensor> dSpT(sphericalTensor::one*1e-2);
+
+        label randomDistributionTestSize = 50000000;
+
+        Info<< nl << "Distribution<sphericalTensor>" << nl
+            << "Sampling "
+            << randomDistributionTestSize
+            << " times from uniform distribution."
+            << endl;
+
+        for (label i = 0; i < randomDistributionTestSize; i++)
+        {
+            dSpT.add(R.sphericalTensor01());
+        }
+
+        Info<< "Mean " << dSpT.mean() << nl
+            << "Median " << dSpT.median()
+            << endl;
+
+        dSpT.write("Distribution_sphericalTensor_test", dSpT.normalised());
+    }
 
     Info<< nl << "End" << nl << endl;
 
diff --git a/src/OpenFOAM/containers/HashTables/Map/Distribution/Distribution.C b/src/OpenFOAM/containers/HashTables/Map/Distribution/Distribution.C
index 3a0defe22c..cee8014f09 100644
--- a/src/OpenFOAM/containers/HashTables/Map/Distribution/Distribution.C
+++ b/src/OpenFOAM/containers/HashTables/Map/Distribution/Distribution.C
@@ -25,16 +25,14 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "Distribution.H"
-
-// * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * * //
-
+#include "OFstream.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Type>
 Foam::Distribution<Type>::Distribution()
 :
-    List< Map<label> >(pTraits<Type>::nComponents),
+    List< Map<scalar> >(pTraits<Type>::nComponents),
     binWidth_(pTraits<Type>::one)
 {}
 
@@ -42,7 +40,7 @@ Foam::Distribution<Type>::Distribution()
 template<class Type>
 Foam::Distribution<Type>::Distribution(const Type& binWidth)
 :
-    List< Map<label> >(pTraits<Type>::nComponents),
+    List< Map<scalar> >(pTraits<Type>::nComponents),
     binWidth_(binWidth)
 {}
 
@@ -53,7 +51,7 @@ Foam::Distribution<Type>::Distribution(const Type& binWidth)
 //     const cmptType& binWidth
 // )
 // :
-//     List< Map<label> >(pTraits<Type>::nComponents),
+//     List< Map<scalar> >(pTraits<Type>::nComponents),
 //     binWidth_(binWidth*pTraits<Type>::one)
 // {}
 
@@ -61,7 +59,7 @@ Foam::Distribution<Type>::Distribution(const Type& binWidth)
 template<class Type>
 Foam::Distribution<Type>::Distribution(const Distribution<Type>& d)
 :
-    List< Map<label> >(static_cast< const List< Map<label> >& >(d)),
+    List< Map<scalar> >(static_cast< const List< Map<scalar> >& >(d)),
     binWidth_(d.binWidth())
 {}
 
@@ -80,13 +78,145 @@ Foam::Distribution<Type>::~Distribution()
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 template<class Type>
-void Foam::Distribution<Type>::add(const Type& valueToAdd)
+Foam::scalar Foam::Distribution<Type>::totalWeight(direction cmpt) const
 {
-    for (direction cmpt=0; cmpt < pTraits<Type>::nComponents; cmpt++)
-    {
-        Map<label>& cmptDistribution = (*this)[cmpt];
+    const Map<scalar>& cmptDistribution = (*this)[cmpt];
 
-        Map<label>::iterator iter(cmptDistribution.begin());
+    scalar sumOfWeights = 0.0;
+
+    forAllConstIter(Map<scalar>, cmptDistribution, iter)
+    {
+        sumOfWeights += iter();
+    }
+
+    return sumOfWeights;
+}
+
+
+template<class Type>
+inline Type Foam::Distribution<Type>::mean() const
+{
+    Type meanValue(pTraits<Type>::zero);
+
+    for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+    {
+        const Map<scalar>& cmptDistribution = (*this)[cmpt];
+
+        scalar totalCmptWeight = totalWeight(cmpt);
+
+        List<label> keys = cmptDistribution.sortedToc();
+
+        forAll(keys,k)
+        {
+            label key = keys[k];
+
+            setComponent(meanValue, cmpt) +=
+                (0.5 + scalar(key))
+               *component(binWidth_, cmpt)
+               *cmptDistribution[key]
+               /totalCmptWeight;
+        }
+    }
+
+    return meanValue;
+}
+
+
+template<class Type>
+inline Type Foam::Distribution<Type>::median()
+{
+    Type medianValue(pTraits<Type>::zero);
+
+    List< List < Pair<scalar> > > normDistribution = normalised();
+
+    for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+    {
+        List<Pair<scalar> >& normDist = normDistribution[cmpt];
+
+        scalar cumulative = 0.0;
+
+        if (normDist.size())
+        {
+            if (normDist.size() == 1)
+            {
+                setComponent(medianValue, cmpt) = normDist[0].first();
+            }
+            else if
+            (
+                normDist.size() > 1
+             && normDist[0].second()*component(binWidth_, cmpt) > 0.5
+            )
+            {
+                scalar xk = normDist[1].first();
+
+                scalar xkm1 = normDist[0].first();
+
+                scalar Sk =
+                    (normDist[0].second() + normDist[1].second())
+                   *component(binWidth_, cmpt);
+
+                scalar Skm1 = normDist[0].second()*component(binWidth_, cmpt);
+
+                setComponent(medianValue, cmpt) =
+                    (0.5 - Skm1)*(xk - xkm1)/(Sk - Skm1) + xkm1;
+            }
+            else
+            {
+                label lastNonZeroIndex = 0;
+
+                forAll(normDist,nD)
+                {
+                    if
+                    (
+                        cumulative
+                      + (normDist[nD].second()*component(binWidth_, cmpt))
+                      > 0.5
+                    )
+                    {
+                        scalar xk = normDist[nD].first();
+
+                        scalar xkm1 = normDist[lastNonZeroIndex].first();
+
+                        scalar Sk =
+                            cumulative
+                          + (normDist[nD].second()*component(binWidth_, cmpt));
+
+                        scalar Skm1 = cumulative;
+
+                        setComponent(medianValue, cmpt) =
+                            (0.5 - Skm1)*(xk - xkm1)/(Sk - Skm1) + xkm1;
+
+                        break;
+                    }
+                    else if (normDist[nD].second() > 0.0)
+                    {
+                        cumulative +=
+                            normDist[nD].second()*component(binWidth_, cmpt);
+
+                        lastNonZeroIndex = nD;
+                    }
+                }
+            }
+        }
+
+    }
+
+    return medianValue;
+}
+
+
+template<class Type>
+void Foam::Distribution<Type>::add
+(
+    const Type& valueToAdd,
+    const Type& weight
+)
+{
+    for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+    {
+        Map<scalar>& cmptDistribution = (*this)[cmpt];
+
+        Map<scalar>::iterator iter(cmptDistribution.begin());
 
         label n =
             label(component(valueToAdd, cmpt)/component(binWidth_, cmpt))
@@ -96,23 +226,23 @@ void Foam::Distribution<Type>::add(const Type& valueToAdd)
 
         if (iter == cmptDistribution.end())
         {
-            cmptDistribution.insert(n,1);
+            cmptDistribution.insert(n, component(weight, cmpt));
         }
         else
         {
-            cmptDistribution[n]++;
+            cmptDistribution[n] += component(weight, cmpt);
         }
 
-        if (cmptDistribution[n] < 0)
-        {
-            FatalErrorIn("Distribution::add(const scalar valueToAdd)")
-                << "Accumulated Distribution value has become negative: "
-                << "bin = " << (0.5 + scalar(n))*component(binWidth_, cmpt)
-                << ", value = " << cmptDistribution[n]
-                << ". This is most likely to be because too many samples "
-                << "have been added to a bin and the label has 'rolled round'"
-                << abort(FatalError);
-        }
+        // if (cmptDistribution[n] < 0)
+        // {
+        //     FatalErrorIn("Distribution::add(const scalar valueToAdd)")
+        //         << "Accumulated Distribution value has become negative: "
+        //         << "bin = " << (0.5 + scalar(n))*component(binWidth_, cmpt)
+        //         << ", value = " << cmptDistribution[n]
+        //         << ". This is most likely to be because too many samples "
+        //         << "have been added to a bin and the weight has 'rolled round'"
+        //         << abort(FatalError);
+        // }
     }
 }
 
@@ -120,11 +250,11 @@ void Foam::Distribution<Type>::add(const Type& valueToAdd)
 template<class Type>
 void Foam::Distribution<Type>::insertMissingKeys()
 {
-    for (direction cmpt=0; cmpt < pTraits<Type>::nComponents; cmpt++)
+    for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
     {
-        Map<label>& cmptDistribution = (*this)[cmpt];
+        Map<scalar>& cmptDistribution = (*this)[cmpt];
 
-        Map<label>::iterator iter(cmptDistribution.begin());
+        Map<scalar>::iterator iter(cmptDistribution.begin());
 
         List<label> keys = cmptDistribution.sortedToc();
 
@@ -146,19 +276,57 @@ void Foam::Distribution<Type>::insertMissingKeys()
 
 template<class Type>
 Foam::List< Foam::List< Foam::Pair<Foam::scalar> > >Foam::
-Distribution<Type>::raw()
+Distribution<Type>::normalised()
 {
-    List< List < Pair<scalar> > > rawDistributions(pTraits<Type>::nComponents);
+    List< List < Pair<scalar> > > normDistribution(pTraits<Type>::nComponents);
 
     insertMissingKeys();
 
-    for (direction cmpt=0; cmpt < pTraits<Type>::nComponents; cmpt++)
+    for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
     {
-        Map<label>& cmptDistribution = (*this)[cmpt];
+        Map<scalar>& cmptDistribution = (*this)[cmpt];
+
+        scalar totalCmptWeight = totalWeight(cmpt);
 
         List<label> keys = cmptDistribution.sortedToc();
 
-        List<Pair<scalar> >& rawDist = rawDistributions[cmpt];
+        List<Pair<scalar> >& normDist = normDistribution[cmpt];
+
+        normDist.setSize(keys.size());
+
+        forAll(keys,k)
+        {
+            label key = keys[k];
+
+            normDist[k].first() =
+                (0.5 + scalar(key))*component(binWidth_, cmpt);
+
+            normDist[k].second() =
+                cmptDistribution[key]
+               /totalCmptWeight
+               /component(binWidth_, cmpt);
+        }
+    }
+
+    return normDistribution;
+}
+
+
+template<class Type>
+Foam::List< Foam::List< Foam::Pair<Foam::scalar> > >Foam::
+Distribution<Type>::raw()
+{
+    List< List < Pair<scalar> > > rawDistribution(pTraits<Type>::nComponents);
+
+    insertMissingKeys();
+
+    for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+    {
+        Map<scalar>& cmptDistribution = (*this)[cmpt];
+
+        List<label> keys = cmptDistribution.sortedToc();
+
+        List<Pair<scalar> >& rawDist = rawDistribution[cmpt];
 
         rawDist.setSize(keys.size());
 
@@ -168,11 +336,48 @@ Distribution<Type>::raw()
 
             rawDist[k].first() = (0.5 + scalar(key))*component(binWidth_, cmpt);
 
-            rawDist[k].second() = scalar(cmptDistribution[key]);
+            rawDist[k].second() = cmptDistribution[key];
         }
     }
 
-    return rawDistributions;
+    return rawDistribution;
+}
+
+
+template<class Type>
+void Foam::Distribution<Type>::write
+(
+    const fileName& filePrefix,
+    const List< List<Pair<scalar> > >& pairs
+) const
+{
+    if (pairs.size() != pTraits<Type>::nComponents)
+    {
+        FatalErrorIn
+        (
+            "Distribution::write"
+            "("
+                "const fileName& filePrefix,"
+                "const List< List<Pair<scalar> > >& pairs"
+            ")"
+        )
+            << "List of pairs (" << pairs.size()
+            << ") is not the same size as the number of components ("
+            << pTraits<Type>::nComponents << ")." << nl
+            << abort(FatalError);
+    }
+
+    for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+    {
+        const List<Pair<scalar> >& cmptPairs = pairs[cmpt];
+
+        OFstream os(filePrefix + '_' + pTraits<Type>::componentNames[cmpt]);
+
+        forAll(cmptPairs, i)
+        {
+            os  << cmptPairs[i].first() << ' ' << cmptPairs[i].second() << nl;
+        }
+    }
 }
 
 
@@ -195,7 +400,7 @@ void Foam::Distribution<Type>::operator=
             << abort(FatalError);
     }
 
-    List< Map<label> >::operator=(rhs);
+    List< Map<scalar> >::operator=(rhs);
 
     binWidth_ = rhs.binWidth();
 }
@@ -211,7 +416,7 @@ Foam::Ostream& Foam::operator<<
 )
 {
     os  << d.binWidth_
-        << static_cast<const List< Map<label> >& >(d);
+        << static_cast<const List< Map<scalar> >& >(d);
 
     // Check state of Ostream
     os.check
diff --git a/src/OpenFOAM/containers/HashTables/Map/Distribution/Distribution.H b/src/OpenFOAM/containers/HashTables/Map/Distribution/Distribution.H
index 0f60ee5511..377ffae830 100644
--- a/src/OpenFOAM/containers/HashTables/Map/Distribution/Distribution.H
+++ b/src/OpenFOAM/containers/HashTables/Map/Distribution/Distribution.H
@@ -60,37 +60,19 @@ Ostream& operator<<(Ostream&, const Distribution<Type>&);
 template<class Type>
 class Distribution
 :
-    public List< Map<label> >
+    public List< Map<scalar> >
 {
     // Private data
 
         Type binWidth_;
 
 
-    // Private Member Functions
-
-        // //- Disallow default bitwise copy construct
-        // Distribution(const Distribution<Type>&);
-
-        // //- Disallow default bitwise assignment
-        // void operator=(const Distribution<Type>&);
-
-
 public:
 
     //- Component type
     typedef typename pTraits<Type>::cmptType cmptType;
 
 
-    // Static data members
-
-        // static void write
-        // (
-        //     const fileName& file,
-        //     const List<Pair<scalar> >& pairs
-        // );
-
-
     // Constructors
 
         //- Construct null
@@ -99,7 +81,7 @@ public:
         //- Construct from separate binWidth for each component
         Distribution(const Type& binWidth);
 
-        // //- Construct from single binWidth for each component
+        //- Construct from single binWidth for each component
         // Distribution(const cmptType& binWidth);
 
         //- Construct as copy
@@ -111,17 +93,24 @@ public:
 
     // Member Functions
 
-        // Type totalEntries() const;
+        scalar totalWeight(direction cmpt) const;
 
-        // Type mean() const;
+        Type mean() const;
 
-        // Type median();
+        // From http://mathworld.wolfram.com/StatisticalMedian.html
+        // The statistical median is the value of the Distribution
+        // variable where the cumulative Distribution = 0.5.
+        Type median();
 
-        void add(const Type& valueToAdd);
+        void add
+        (
+            const Type& valueToAdd,
+            const Type& weight = pTraits<Type>::one
+        );
 
         void insertMissingKeys();
 
-        // List<Pair<scalar> > normalised();
+        List< List<Pair<scalar> > > normalised();
 
         // List<Pair<scalar> > normalisedMinusMean();
 
@@ -133,15 +122,19 @@ public:
 
         inline const Type& binWidth() const;
 
+    // Write
+
+        void write
+        (
+            const fileName& filePrefix,
+            const List< List<Pair<scalar> > >& pairs
+        ) const;
+
+
     // Member Operators
 
         void operator=(const Distribution<Type>&);
 
-
-    // Friend Functions
-
-    // Friend Operators
-
     // IOstream Operators
 
         friend Ostream& operator<< <Type>