added kernel density wrapper

added general kld solution
fixed minor bugs
added tests

math/distribution/KernelDensity.h

@@ -0,0 +1,35 @@
+#ifndef KERNELDENSITY_H
+#define KERNELDENSITY_H
+
+#include <cmath>
+#include <random>
+#include <functional>
+
+#include <eigen3/Eigen/Dense>
+
+#include "../../Assertions.h"
+#include "../Random.h"
+
+
+namespace Distribution {
+
+    template <typename T, typename Sample> class KernelDensity{
+
+    private:
+        const std::function<T(Sample)> probabilityFunction;
+
+
+    public:
+        KernelDensity(const std::function<T(Sample)> probabilityFunction) : probabilityFunction(probabilityFunction){
+
+        }
+
+        T getProbability(Sample sample){
+            return probabilityFunction(sample);
+        }
+
+    };
+
+}
+
+#endif // KERNELDENSITY_H

math/distribution/NormalN.h

@@ -15,8 +15,8 @@ namespace Distribution {
 
     private:
 
-        const Eigen::VectorXd mu;
-        const Eigen::MatrixXd sigma;
+        Eigen::VectorXd mu;
+        Eigen::MatrixXd sigma;
 
         const double _a;
         const Eigen::MatrixXd _sigmaInv;
@@ -61,6 +61,14 @@ namespace Distribution {
             return this->_sigmaInv;
         }
 
+        void setSigma(Eigen::MatrixXd sigma){
+            this->sigma = sigma;
+        }
+
+        void setMu(Eigen::VectorXd mu){
+            this->mu = mu;
+        }
+
         /** return a NormalN based on given data */
         static NormalDistributionN getNormalNFromSamples(const Eigen::MatrixXd& data) {
 

math/divergence/KullbackLeibler.h

@@ -9,10 +9,77 @@
 
 namespace Divergence {
 
+    enum LOGMODE{
+        BASE2,
+        BASE10,
+        NATURALIS
+    };
+
     template <typename Scalar> class KullbackLeibler {
 
         public:
 
+            /** Calculate the Kullback Leibler Distance from a set of sample densities
+             *  Info: https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence
+             * @param P is the vector containing the densities of a set of samples
+             * @param Q is a vector containg the densities of the same samples set then P
+            */
+            static inline Scalar getGeneralFromSamples(Eigen::VectorXd P, Eigen::VectorXd Q, LOGMODE mode){
+
+                // Assure P and Q have the same size and are finite in all values
+                Assert::equal(P.size(), Q.size(), "The sample vectors P and Q do not have the same size.");
+                if(!P.allFinite() || !Q.allFinite()){
+                    Assert::doThrow("The sample vectors P and Q are not finite. Check for NaN or Inf.");
+                }
+
+                Scalar dist = 0.0;
+                Scalar PQratio = 0.0;
+
+                //normalize to 1
+                P /= P.sum();
+                Q /= Q.sum();
+
+                Assert::isNear((double)P.sum(), 1.0, 0.01,"Normalization failed.. this shouldn't happen");
+                Assert::isNear((double)Q.sum(), 1.0, 0.01, "Normalization failed.. this shouldn't happen");
+
+                //sum up the logarithmic difference between P and Q, also called kullback leibler...
+                for(int i = 0; i < P.size(); ++i){
+
+                    if((P[i] == 0.0) && (Q[i] == 0.0)){
+                            dist += 0.0;
+                    } else {
+
+                        // calc PQratio
+                        if(Q[i] == 0.0){
+                            Assert::doThrow("Division by zero is not allowed ;). TODO: What if the densities are to small?");
+                            //PQratio = P[i] / 0.00001;
+                        } else {
+                            PQratio = P[i] / Q[i];
+                        }
+
+                        //use log for dist
+                        if (mode == NATURALIS){
+                            dist += P[i] * log(PQratio);
+                        }
+
+                        else if (mode == BASE2){
+                            dist += P[i] * log2(PQratio);
+                        }
+
+                        else if (mode == BASE10){
+                            dist += P[i] * log10(PQratio);
+                        }
+
+                    }
+                }
+
+                return dist;
+            }
+
+            static inline Scalar getGeneralFromSamplesSymmetric(Eigen::VectorXd P, Eigen::VectorXd Q, LOGMODE mode){
+                return getGeneralFromSamples(P, Q, mode) + getGeneralFromSamples(Q, P, mode);
+            }
+
             /** Calculate the Kullback Leibler Distance for a univariate Gaussian distribution
             *   Info: https://tgmstat.wordpress.com/2013/07/10/kullback-leibler-divergence/
             */