diff --git a/distributions/univariate/continuous/ContinuousDistribution.cpp b/distributions/univariate/continuous/ContinuousDistribution.cpp
index 74a2046c..617cb4c6 100644
--- a/distributions/univariate/continuous/ContinuousDistribution.cpp
+++ b/distributions/univariate/continuous/ContinuousDistribution.cpp
@@ -10,6 +10,15 @@ void ContinuousDistribution::ProbabilityDensityFunction(const std::vector<double
         y[i] = f(x[i]);
 }
 
+void ContinuousDistribution::LogProbabilityDensityFunction(const std::vector<double> &x, std::vector<double> &y) const
+{
+    size_t size = x.size();
+    if (size > y.size())
+        return;
+    for (size_t i = 0; i != size; ++i)
+        y[i] = logf(x[i]);
+}
+
 double ContinuousDistribution::quantileImpl(double p) const
 {
     double guess = 0.0;
diff --git a/distributions/univariate/continuous/ContinuousDistribution.h b/distributions/univariate/continuous/ContinuousDistribution.h
index c38af6a4..3f9857b6 100644
--- a/distributions/univariate/continuous/ContinuousDistribution.h
+++ b/distributions/univariate/continuous/ContinuousDistribution.h
@@ -36,6 +36,14 @@ class RANDLIBSHARED_EXPORT ContinuousDistribution : public virtual UnivariateDis
      */
     void ProbabilityDensityFunction(const std::vector<double> &x, std::vector<double> &y) const;
 
+    /**
+     * @brief LogProbabilityDensityFunction
+     * fill vector y by logf(x)
+     * @param x
+     * @param y
+     */
+    void LogProbabilityDensityFunction(const std::vector<double> &x, std::vector<double> &y) const;
+
     double Mode() const override;
 
 protected:
diff --git a/distributions/univariate/discrete/DiscreteDistribution.cpp b/distributions/univariate/discrete/DiscreteDistribution.cpp
index d5ce0637..34c79116 100644
--- a/distributions/univariate/discrete/DiscreteDistribution.cpp
+++ b/distributions/univariate/discrete/DiscreteDistribution.cpp
@@ -10,6 +10,15 @@ void DiscreteDistribution::ProbabilityMassFunction(const std::vector<int> &x, st
         y[i] = P(x[i]);
 }
 
+void DiscreteDistribution::LogProbabilityMassFunction(const std::vector<int> &x, std::vector<double> &y) const
+{
+    size_t size = x.size();
+    if (size > y.size())
+        return;
+    for (size_t i = 0; i != size; ++i)
+        y[i] = logP(x[i]);
+}
+
 int DiscreteDistribution::Mode() const
 {
     /// Works only for unimodal and monotone from starting point to the mode distributions
diff --git a/distributions/univariate/discrete/DiscreteDistribution.h b/distributions/univariate/discrete/DiscreteDistribution.h
index 5f5a9658..7d07eb53 100644
--- a/distributions/univariate/discrete/DiscreteDistribution.h
+++ b/distributions/univariate/discrete/DiscreteDistribution.h
@@ -36,6 +36,14 @@ class RANDLIBSHARED_EXPORT DiscreteDistribution : public virtual UnivariateDistr
      */
     void ProbabilityMassFunction(const std::vector<int> &x, std::vector<double> &y) const;
 
+    /**
+     * @brief LogProbabilityMassFunction
+     * fill vector y with logP(x)
+     * @param x
+     * @param y
+     */
+    void LogProbabilityMassFunction(const std::vector<int> &x, std::vector<double> &y) const;
+
     int Mode() const override;
 
 private: