diff --git a/Tudat/Mathematics/Statistics/CMakeLists.txt b/Tudat/Mathematics/Statistics/CMakeLists.txt
index b3a9c5daf5..4fb974b620 100755
--- a/Tudat/Mathematics/Statistics/CMakeLists.txt
+++ b/Tudat/Mathematics/Statistics/CMakeLists.txt
@@ -37,12 +37,14 @@
 set(STATISTICS_SOURCES
   "${SRCROOT}${MATHEMATICSDIR}/Statistics/basicStatistics.cpp"
   "${SRCROOT}${MATHEMATICSDIR}/Statistics/simpleLinearRegression.cpp"
+  "${SRCROOT}${MATHEMATICSDIR}/Statistics/randomSampling.cpp"
 )
 
 # Add header files.
 set(STATISTICS_HEADERS 
   "${SRCROOT}${MATHEMATICSDIR}/Statistics/basicStatistics.h"
   "${SRCROOT}${MATHEMATICSDIR}/Statistics/simpleLinearRegression.h"
+  "${SRCROOT}${MATHEMATICSDIR}/Statistics/randomSampling.h"
 )
 
 # Add static libraries.
@@ -57,3 +59,7 @@ target_link_libraries(test_SimpleLinearRegression tudat_statistics ${Boost_LIBRA
 add_executable(test_BasicStatistics "${SRCROOT}${MATHEMATICSDIR}/Statistics/UnitTests/unitTestBasicStatistics.cpp")
 setup_custom_test_program(test_BasicStatistics "${SRCROOT}${MATHEMATICSDIR}/Statistics")
 target_link_libraries(test_BasicStatistics tudat_statistics ${Boost_LIBRARIES})
+
+add_executable(test_RandomSampling "${SRCROOT}${MATHEMATICSDIR}/Statistics/UnitTests/unitTestRandomSampling.cpp")
+setup_custom_test_program(test_RandomSampling "${SRCROOT}${MATHEMATICSDIR}/Statistics")
+target_link_libraries(test_RandomSampling tudat_statistics ${Boost_LIBRARIES})
diff --git a/Tudat/Mathematics/Statistics/UnitTests/unitTestRandomSampling.cpp b/Tudat/Mathematics/Statistics/UnitTests/unitTestRandomSampling.cpp
new file mode 100644
index 0000000000..086add5cba
--- /dev/null
+++ b/Tudat/Mathematics/Statistics/UnitTests/unitTestRandomSampling.cpp
@@ -0,0 +1,181 @@
+/*    Copyright (c) 2010-2016, Delft University of Technology
+ *    All rights reserved.
+ *
+ *    Redistribution and use in source and binary forms, with or without modification, are
+ *    permitted provided that the following conditions are met:
+ *      - Redistributions of source code must retain the above copyright notice, this list of
+ *        conditions and the following disclaimer.
+ *      - Redistributions in binary form must reproduce the above copyright notice, this list of
+ *        conditions and the following disclaimer in the documentation and/or other materials
+ *        provided with the distribution.
+ *      - Neither the name of the Delft University of Technology nor the names of its contributors
+ *        may be used to endorse or promote products derived from this software without specific
+ *        prior written permission.
+ *
+ *    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS
+ *    OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+ *    MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ *    COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ *    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
+ *    GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
+ *    AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ *    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
+ *    OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ *    Changelog
+ *      YYMMDD    Author            Comment
+ *      160902    R. Hoogendoorn    File created.
+ *
+ *
+ *    References
+ *
+ *    Notes
+ *
+ */
+
+#define BOOST_TEST_MAIN
+
+#include <iostream> // cout sometimes needs this
+#include <vector>
+#include <limits>
+
+#include <Eigen/Core>
+
+
+#include <boost/bind.hpp>
+#include <boost/function.hpp>
+#include <boost/shared_ptr.hpp>
+#include <boost/make_shared.hpp>
+
+#include <boost/test/unit_test.hpp>
+#include "tudat/Basics/testMacros.h"
+
+#include "tudat/Mathematics/BasicMathematics/mathematicalConstants.h"
+#include "tudat/Mathematics/Statistics/randomSampling.h"
+#include <tudat/Mathematics/Statistics/basicStatistics.h>
+
+namespace tudat
+{
+namespace unit_tests
+{
+
+BOOST_AUTO_TEST_SUITE( test_Random_Sampler )
+
+//using tudat::mathematical_constants::PI;
+
+//BOOST_AUTO_TEST_CASE( test_Sobol_Sampler )
+//{
+//    int dimension = 3 ;
+//    int numberOfSamples = 1E6 ;
+
+//    Eigen::VectorXd lower(dimension);
+//    Eigen::VectorXd upper(dimension);
+//    lower << 0.0 , 1.0 , -2.0 ;
+//    upper << 1.0 , 3.0 , 4.0 ;
+
+//    Eigen::VectorXd width = upper - lower ;
+//    Eigen::VectorXd average = (upper + lower) / 2.0 ;
+
+//    std::vector<Eigen::VectorXd> sobolSamples = Thesis::Statistics::sobolSamplerXd(dimension,numberOfSamples,lower,upper);
+
+//    // Compute sample average
+//    Eigen::VectorXd sampleAverage(dimension);
+//    sampleAverage.setZero(dimension,1);
+//    for( unsigned int i = 0 ; i < sobolSamples.size() ; i++){
+//        sampleAverage += sobolSamples[i] ;
+//    }
+//    sampleAverage = sampleAverage / double(numberOfSamples) ;
+
+//    BOOST_CHECK_SMALL( std::fabs( average(0) - sampleAverage(0) ) , 2.0E-6 );
+//    BOOST_CHECK_SMALL( std::fabs( average(1) - sampleAverage(1) ) , 2.0E-6 );
+//    BOOST_CHECK_SMALL( std::fabs( average(2) - sampleAverage(2) ) , 2.0E-6 );
+
+//    sobolSamples = Thesis::Statistics::sobolSamplerXd(dimension,numberOfSamples);
+
+//    // Compute sample average
+//    sampleAverage.setZero(dimension,1);
+//    for( unsigned int i = 0 ; i < sobolSamples.size() ; i++){
+//        sampleAverage += sobolSamples[i] ;
+//    }
+//    sampleAverage = sampleAverage / double(numberOfSamples) ;
+
+//    BOOST_CHECK_SMALL( std::fabs( 0.5 - sampleAverage(0) ) , 2.0E-6 );
+//    BOOST_CHECK_SMALL( std::fabs( 0.5 - sampleAverage(1) ) , 2.0E-6 );
+//    BOOST_CHECK_SMALL( std::fabs( 0.5 - sampleAverage(2) ) , 2.0E-6 );
+//}
+
+BOOST_AUTO_TEST_CASE( test_randomVectorUniform )
+{
+    int dimension = 3 ;
+    int numberOfSamples = 1E6 ;
+    int seed = 511;
+
+    Eigen::VectorXd lower(dimension);
+    Eigen::VectorXd upper(dimension);
+    lower << 0.0 , 1.0 , -2.0 ;
+    upper << 1.0 , 3.0 , 4.0 ;
+
+    Eigen::VectorXd width = upper - lower ;
+    Eigen::VectorXd average = (upper + lower) / 2.0 ;
+
+    std::vector<Eigen::VectorXd> samples =
+            tudat::statistics::generateRandomVectorUniform( seed , numberOfSamples, lower, upper );
+
+    // Compute sample average
+    Eigen::VectorXd sampleAverage(dimension);
+    sampleAverage.setZero(dimension,1);
+    for( unsigned int i = 0 ; i < samples.size() ; i++){
+        sampleAverage += samples[i] ;
+    }
+    sampleAverage = sampleAverage / double(numberOfSamples) ;
+
+    BOOST_CHECK_SMALL( std::fabs( average(0) - sampleAverage(0) ) , 2.0E-3 );
+    BOOST_CHECK_SMALL( std::fabs( average(1) - sampleAverage(1) ) , 2.0E-3 );
+    BOOST_CHECK_SMALL( std::fabs( average(2) - sampleAverage(2) ) , 2.0E-3 );
+
+    samples = tudat::statistics::generateRandomVectorUniform( seed , numberOfSamples, dimension );
+
+    // Compute sample average
+    sampleAverage.setZero(dimension,1);
+    for( unsigned int i = 0 ; i < samples.size() ; i++){
+        sampleAverage += samples[i] ;
+    }
+    sampleAverage = sampleAverage / double(numberOfSamples) ;
+
+    BOOST_CHECK_SMALL( std::fabs( 0.5 - sampleAverage(0) ) , 2.0E-4 );
+    BOOST_CHECK_SMALL( std::fabs( 0.5 - sampleAverage(1) ) , 2.0E-4 );
+    BOOST_CHECK_SMALL( std::fabs( 0.5 - sampleAverage(2) ) , 2.0E-4 );
+}
+
+BOOST_AUTO_TEST_CASE( test_randomVectorGaussian )
+{
+    int dimension = 3 ;
+    int numberOfSamples = 1E6 ;
+    int seed = 511;
+
+    Eigen::VectorXd mean(dimension);
+    Eigen::VectorXd standardDev(dimension);
+    mean << 0.0 , 1.0 , -2.0 ;
+    standardDev << 1.0 , 3.0 , 4.0 ;
+
+    std::vector<Eigen::VectorXd> samples =
+            tudat::statistics::generateRandomVectorGaussian( seed , numberOfSamples, dimension, mean, standardDev );
+
+    // Compute sample average
+    Eigen::VectorXd sampleAverage(dimension);
+    sampleAverage.setZero(dimension,1);
+    for( unsigned int i = 0 ; i < samples.size() ; i++){
+        sampleAverage += samples[i] ;
+    }
+    sampleAverage = sampleAverage / double(numberOfSamples) ;
+
+    BOOST_CHECK_SMALL( std::fabs( mean(0) - sampleAverage(0) ) , 2.0E-3 );
+    BOOST_CHECK_SMALL( std::fabs( mean(1) - sampleAverage(1) ) , 2.0E-3 );
+    BOOST_CHECK_SMALL( std::fabs( mean(2) - sampleAverage(2) ) , 2.0E-3 );
+}
+
+
+BOOST_AUTO_TEST_SUITE_END( )
+
+} // namespace unit_tests
+} // namespace tudat
diff --git a/Tudat/Mathematics/Statistics/randomSampling.cpp b/Tudat/Mathematics/Statistics/randomSampling.cpp
new file mode 100644
index 0000000000..a094d343c6
--- /dev/null
+++ b/Tudat/Mathematics/Statistics/randomSampling.cpp
@@ -0,0 +1,191 @@
+/*    Copyright (c) 2010-2015, Delft University of Technology
+ *    All rights reserved.
+ *
+ *    Redistribution and use in source and binary forms, with or without modification, are
+ *    permitted provided that the following conditions are met:
+ *      - Redistributions of source code must retain the above copyright notice, this list of
+ *        conditions and the following disclaimer.
+ *      - Redistributions in binary form must reproduce the above copyright notice, this list of
+ *        conditions and the following disclaimer in the documentation and/or other materials
+ *        provided with the distribution.
+ *      - Neither the name of the Delft University of Technology nor the names of its contributors
+ *        may be used to endorse or promote products derived from this software without specific
+ *        prior written permission.
+ *
+ *    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS
+ *    OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+ *    MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ *    COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ *    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
+ *    GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
+ *    AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ *    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
+ *    OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ *    Changelog
+ *      YYMMDD    Author            Comment
+ *      160902    R. Hoogendoorn    File created.
+ *
+ *
+ *    References
+ *
+ *    Notes
+ *
+ */
+
+
+#include <tudat/Mathematics/Statistics/randomSampling.h>
+
+#include <boost/make_shared.hpp>
+#include <cmath>
+
+// Sobol sampler (Not in header, because it polutes.. )
+//#include <gsl/gsl_qrng.h>   // GSL Quasi Random generators -> Sobol
+
+// Pseudo random sampler
+#include <boost/random.hpp> // Random generator
+
+namespace tudat
+{
+
+namespace statistics
+{
+
+//! Generator random vector using pseudo random generator
+std::vector<Eigen::VectorXd> generateRandomVectorUniform(int seed, int numberOfSamples, int Dimension)
+{
+    // Setup Random generator
+    typedef boost::mt19937 RandomGeneratorType; // Mersenne Twister
+    RandomGeneratorType randomGenerator(seed);              // Create random generator
+
+    boost::uniform_real<> uniformDistribution(0.0 , 1.0); //
+    boost::variate_generator< RandomGeneratorType, boost::uniform_real<> >
+            Dice(randomGenerator, uniformDistribution); // define random generator
+
+    std::vector< Eigen::VectorXd > randomSamples(numberOfSamples) ;
+    Eigen::VectorXd randomSample( Dimension ) ;
+
+    // Sample
+    for(int i = 0 ; i < numberOfSamples ; i++ ) // Generate N samples
+    {
+        for(int j = 0 ; j < randomSample.rows() ; j++) // Generate vector of samples
+        {
+            randomSample(j) = Dice()  ;
+        }
+        randomSamples[i] = randomSample ;
+    }
+    return randomSamples;
+}
+
+//! Generator random vector using pseudo random generator
+std::vector<Eigen::VectorXd> generateRandomVectorUniform(int seed, int numberOfSamples,
+                                             Eigen::VectorXd lowerBound, Eigen::VectorXd upperBound){
+    // Compute properties
+    Eigen::VectorXd width = upperBound - lowerBound;
+    Eigen::VectorXd average = (upperBound + lowerBound)/2.0 ;
+
+    // Setup Random generator
+    typedef boost::mt19937 RandomGeneratorType; // Mersenne Twister
+    RandomGeneratorType randomGenerator(seed);              // Create random generator
+
+    boost::uniform_real<> uniformDistribution(0.0 , 1.0); //
+    boost::variate_generator< RandomGeneratorType, boost::uniform_real<> >
+            Dice(randomGenerator, uniformDistribution); // define random generator
+
+    std::vector< Eigen::VectorXd > randomSamples(numberOfSamples) ;
+    Eigen::VectorXd randomSample( lowerBound.rows() ) ;
+
+    // Sample
+    for(int i = 0 ; i < numberOfSamples ; i++ ) // Generate N samples
+    {
+        for(int j = 0 ; j < randomSample.rows() ; j++) // Generate vector of samples
+        {
+            randomSample(j) = Dice() - 0.5 ;
+        }
+        randomSamples[i] = randomSample.cwiseProduct(width) + average ;
+    }
+    return randomSamples;
+}
+
+////! Generator random vector using Sobol sampler
+//std::vector< Eigen::VectorXd > sobolSamplerXd(const int Dimension, int numberOfSamples,
+//                                              Eigen::VectorXd lowerBound, Eigen::VectorXd upperBound){
+//    // Compute properties
+//    Eigen::VectorXd width = upperBound - lowerBound;
+//    Eigen::VectorXd average = (upperBound + lowerBound)/2.0 ;
+
+//    std::vector< Eigen::VectorXd > sobolSamples( numberOfSamples );
+
+//    Eigen::VectorXd randomSample( Dimension ) ;
+//    double randomSampleArray[ Dimension ];
+
+//    gsl_qrng * q = gsl_qrng_alloc (gsl_qrng_sobol, Dimension );
+
+//    for(int j = 0 ; j < numberOfSamples ; j++){ // Loop over samples
+//        gsl_qrng_get (q, randomSampleArray); // Generate sobol [0,1]
+
+//        for(int i = 0 ; i < Dimension ; i++){ // Fill vector
+//            randomSample(i) = randomSampleArray[i] - 0.5 ;
+//        }
+//        sobolSamples[j] = randomSample.cwiseProduct(width) + average ; // Save vector
+//    }
+
+//    gsl_qrng_free (q); // don't know why?
+
+//    return sobolSamples;
+//}
+
+////! Generator random vector using Sobol sampler
+//std::vector< Eigen::VectorXd > sobolSamplerXd(const int Dimension, int numberOfSamples){
+//    std::vector< Eigen::VectorXd > sobolSamples( numberOfSamples );
+
+//    Eigen::VectorXd randomSample( Dimension ) ;
+//    double randomSampleArray[ Dimension ];
+
+//    gsl_qrng * q = gsl_qrng_alloc (gsl_qrng_sobol, Dimension );
+
+//    for(int j = 0 ; j < numberOfSamples ; j++){ // Loop over samples
+//        gsl_qrng_get (q, randomSampleArray); // Generate sobol [0,1]
+
+//        for(int i = 0 ; i < Dimension ; i++){ // Fill vector
+//            randomSample(i) = randomSampleArray[i] ;
+//        }
+//        sobolSamples[j] = randomSample ; // Save vector
+//    }
+
+//    gsl_qrng_free (q); // don't know why?
+
+//    return sobolSamples;
+//}
+
+//! Generator random vector using pseudo random generator with gaussian distribution (without correlation)
+std::vector< Eigen::VectorXd > generateRandomVectorGaussian(
+        int seed, int numberOfSamples, int Dimension, Eigen::VectorXd mean, Eigen::VectorXd standardDeviation )
+{
+    // Setup Random generator
+    typedef boost::mt19937 RandomGeneratorType; // Mersenne Twister
+    RandomGeneratorType randomGenerator(seed);              // Create random generator
+
+    boost::normal_distribution<> standardNormalDistribution(0.0 , 1.0); // standard normal
+    boost::variate_generator< RandomGeneratorType, boost::normal_distribution<> >
+            Dice(randomGenerator, standardNormalDistribution); // define random generator
+
+    std::vector< Eigen::VectorXd > randomSamples(numberOfSamples) ;
+    Eigen::VectorXd randomSample( Dimension ) ;
+
+    // Sample
+    for(int i = 0 ; i < numberOfSamples ; i++ ) // Generate N samples
+    {
+        for(int j = 0 ; j < randomSample.rows() ; j++) // Generate vector of samples
+        {
+            randomSample(j) = Dice()*standardDeviation(j)  + mean(j)  ;
+        }
+        randomSamples[i] = randomSample ;
+    }
+    return randomSamples;
+}
+
+} // Close Namespace statistics
+
+} // Close Namespace tudat
+
diff --git a/Tudat/Mathematics/Statistics/randomSampling.h b/Tudat/Mathematics/Statistics/randomSampling.h
new file mode 100644
index 0000000000..e9c9eb6121
--- /dev/null
+++ b/Tudat/Mathematics/Statistics/randomSampling.h
@@ -0,0 +1,39 @@
+#ifndef TUDAT_RANDOM_SAMPLING_H
+#define TUDAT_RANDOM_SAMPLING_H
+
+#include <iostream> // cout sometimes needs this
+#include <map>
+
+#include <Eigen/Core>
+
+#include <boost/shared_ptr.hpp>
+
+namespace tudat
+{
+
+namespace statistics
+{
+
+//! Generator random vector using pseudo random generator
+std::vector< Eigen::VectorXd > generateRandomVectorUniform(int seed, int numberOfSamples, int Dimension);
+
+//! Generator random vector using pseudo random generator
+std::vector< Eigen::VectorXd > generateRandomVectorUniform(int seed, int numberOfSamples,
+                                             Eigen::VectorXd lowerBound, Eigen::VectorXd upperBound);
+
+////! Generator random vector using Sobol sampler
+//std::vector< Eigen::VectorXd > sobolSamplerXd(const int Dimension, int numberOfSamples,
+//                                              Eigen::VectorXd lowerBound, Eigen::VectorXd upperBound);
+
+////! Generator random vector using Sobol sampler
+//std::vector< Eigen::VectorXd > sobolSamplerXd(const int Dimension, int numberOfSamples);
+
+//! Generator random vector using pseudo random generator with gaussian distribution (without correlation)
+std::vector< Eigen::VectorXd > generateRandomVectorGaussian(int seed, int numberOfSamples, int Dimension,
+                                    Eigen::VectorXd mean , Eigen::VectorXd standardDeviation );
+
+} // Close Namespace statistics
+
+} // Close Namespace tudat
+
+#endif // TUDAT_RANDOM_SAMPLING_H