diff --git a/lib/src/Uncertainty/Distribution/CMakeLists.txt b/lib/src/Uncertainty/Distribution/CMakeLists.txt
index 4b50220edb..d105a448f3 100644
--- a/lib/src/Uncertainty/Distribution/CMakeLists.txt
+++ b/lib/src/Uncertainty/Distribution/CMakeLists.txt
@@ -29,6 +29,7 @@ ot_add_source_file (ClaytonCopulaFactory.cxx)
 ot_add_source_file (BlockIndependentCopula.cxx)
 ot_add_source_file (JointDistribution.cxx)
 ot_add_source_file (CompositeDistribution.cxx)
+ot_add_source_file (CombinationsDistribution.cxx)
 ot_add_source_file (ConditionalDistribution.cxx)
 ot_add_source_file (CumulativeDistributionNetwork.cxx)
 ot_add_source_file (Dirac.cxx)
@@ -209,6 +210,7 @@ ot_install_header_file (ClaytonCopulaFactory.hxx)
 ot_install_header_file (BlockIndependentCopula.hxx)
 ot_install_header_file (JointDistribution.hxx)
 ot_install_header_file (CompositeDistribution.hxx)
+ot_install_header_file (CombinationsDistribution.hxx)
 ot_install_header_file (ConditionalDistribution.hxx)
 ot_install_header_file (CumulativeDistributionNetwork.hxx)
 ot_install_header_file (DiracFactory.hxx)
diff --git a/lib/src/Uncertainty/Distribution/CombinationsDistribution.cxx b/lib/src/Uncertainty/Distribution/CombinationsDistribution.cxx
new file mode 100644
index 0000000000..6bb34e18b4
--- /dev/null
+++ b/lib/src/Uncertainty/Distribution/CombinationsDistribution.cxx
@@ -0,0 +1,468 @@
+//                                               -*- C++ -*-
+/**
+ *  @brief The CombinationsDistribution distribution
+ *
+ *  Copyright 2005-2024 Airbus-EDF-IMACS-ONERA-Phimeca
+ *
+ *  This library is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU Lesser General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This library is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU Lesser General Public License for more details.
+ *
+ *  You should have received a copy of the GNU Lesser General Public License
+ *  along with this library.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+#include <cmath>
+#include "openturns/Collection.hxx"
+#include "openturns/CombinationsDistribution.hxx"
+#include "openturns/UserDefined.hxx"
+#include "openturns/Combinations.hxx"
+#include "openturns/RandomGenerator.hxx"
+#include "openturns/SpecFunc.hxx"
+#include "openturns/Exception.hxx"
+#include "openturns/PersistentObjectFactory.hxx"
+
+BEGIN_NAMESPACE_OPENTURNS
+
+typedef Collection<UnsignedInteger>     UnsignedIntegerCollection;
+
+CLASSNAMEINIT(CombinationsDistribution)
+
+static const Factory<CombinationsDistribution> Factory_CombinationsDistribution;
+
+/* Default constructor */
+CombinationsDistribution::CombinationsDistribution()
+  : DiscreteDistribution()
+  , k_(0)
+  , n_(0)
+{
+  setName("CombinationsDistribution");
+  setKN(1, 1);
+}
+
+/* Parameters constructor */
+CombinationsDistribution::CombinationsDistribution(const UnsignedInteger k,
+    const UnsignedInteger n)
+  : DiscreteDistribution()
+  , k_(0)
+  , n_(0)
+{
+  setName("CombinationsDistribution");
+  setKN(k, n);
+}
+
+/* Comparison operator */
+Bool CombinationsDistribution::operator ==(const CombinationsDistribution & other) const
+{
+  if (this == &other) return true;
+  return (k_ == other.k_) && (n_ == other.n_);
+}
+
+Bool CombinationsDistribution::equals(const DistributionImplementation & other) const
+{
+  const CombinationsDistribution* p_other = dynamic_cast<const CombinationsDistribution*>(&other);
+  return p_other && (*this == *p_other);
+}
+
+/* String converter */
+String CombinationsDistribution::__repr__() const
+{
+  OSS oss;
+  oss << "class=" << CombinationsDistribution::GetClassName()
+      << " name=" << getName()
+      << " dimension=" << getDimension()
+      << " k=" << k_
+      << " n=" << n_;
+  return oss;
+}
+
+String CombinationsDistribution::__str__(const String & ) const
+{
+  OSS oss;
+  oss << getClassName() << "(k = " << k_ << ", n = " << n_ << ")";
+  return oss;
+}
+
+/* Virtual constructor */
+CombinationsDistribution * CombinationsDistribution::clone() const
+{
+  return new CombinationsDistribution(*this);
+}
+
+/* Compute the numerical range of the distribution given the parameters values */
+void CombinationsDistribution::computeRange()
+{
+  Point lowerBound(k_);
+  Point upperBound(k_);
+  for (UnsignedInteger i = 0; i < k_; ++i)
+    {
+      lowerBound[i] = i;
+      upperBound[i] = n_ - k_ + i;
+    }
+  const Interval::BoolCollection finiteLowerBound(k_, true);
+  const Interval::BoolCollection finiteUpperBound(k_, true);
+  setRange(Interval(lowerBound, upperBound, finiteLowerBound, finiteUpperBound));
+}
+
+/* Get one realization of the distribution
+   See https://cs.stackexchange.com/questions/104930/efficient-n-choose-k-random-sampling
+*/
+Point CombinationsDistribution::getRealization() const
+{
+  Indices flags(n_);
+  Indices integralRealization;
+  Bool done = false;
+  // Build the complementary subset if its cardinal is smaller than k_
+  const UnsignedInteger actualK = (k_ > n_ / 2 ? n_ - k_ : k_);
+  while (integralRealization.getSize() < actualK)
+    {
+      const UnsignedInteger i = RandomGenerator::IntegerGenerate(n_);
+      if (flags[i] == 0)
+        {
+          integralRealization.add(i);
+          flags[i] = 1;
+        }
+    } // while
+  // Di I built the complementary set?
+  if (actualK != k_)
+    integralRealization = integralRealization.complement(n_);
+  // The realization must be sorted in ascending order
+  std::sort(integralRealization.begin(), integralRealization.end());
+  Point realization(k_);
+  for (UnsignedInteger i = 0; i < k_; ++i)
+    realization[i] = integralRealization[i];
+  return realization;
+}
+
+/* Get the PDF of the distribution */
+Scalar CombinationsDistribution::computeLogPDF(const Point & point) const
+{
+  if (point.getDimension() != k_) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=" << k_ << ", here dimension=" << point.getDimension();
+  Indices x(k_);
+  for (UnsignedInteger i = 0; i < k_; ++i)
+  {
+    const Scalar k = point[i];
+    if ((k < i - supportEpsilon_) || (k > n_ - k_ + i + supportEpsilon_)) return SpecFunc::LowestScalar;
+    const UnsignedInteger ik = static_cast< UnsignedInteger > (round(k));
+    if (std::abs(k - ik) > supportEpsilon_) return SpecFunc::LowestScalar;
+    x[i] = ik;
+  }
+  if (!x.isStrictlyIncreasing()) return SpecFunc::LowestScalar;
+  return logPDFValue_;
+}
+
+Scalar CombinationsDistribution::computePDF(const Point & point) const
+{
+  const Scalar logPDF = computeLogPDF(point);
+  if (logPDF == SpecFunc::LowestScalar) return 0.0;
+  return std::exp(logPDF);
+}
+
+Scalar CombinationsDistribution::exploreTree(const Scalar j,
+                                             const UnsignedInteger lower,
+                                             const UnsignedInteger upper,
+                                             const UnsignedInteger count,
+                                             const Point & xReduced) const
+{
+  const Scalar normalizationProba = (k_ - j) / (n_ - j);
+  // Upper branch of the tree
+  const UnsignedInteger lower1 = lower;
+  const UnsignedInteger upper1 = upper;
+  // The value lower1 == dimension is a guard, telling us that the lower bound is 0 and not a component of xReduced
+  const Scalar a1 = (lower1 < xReduced.getDimension()) ? xReduced[lower1] : -1.0;
+  const Scalar b1 = xReduced[upper1];
+  const UnsignedInteger count1 = count;
+  const Scalar f1 = (b1 - a1 - count1) / (count1 + 1);
+  // Lower branch of the tree
+  const UnsignedInteger lower2 = upper;
+  const UnsignedInteger upper2 = j;
+  const Scalar a2 = xReduced[lower2];
+  const Scalar b2 = xReduced[upper2];
+  const UnsignedInteger count2 = 1;
+  const Scalar f2 = b2 - a2;
+  if (j == k_ - 1)
+    return (f1 + f2) * normalizationProba;
+  Scalar value = 0.0;
+  // This test allows one to cut upper parts of the tree
+  if (f1 > 0.0)
+    value += f1 * exploreTree(j + 1, lower1, upper1, count1 + 1, xReduced);
+  // This test allows one to cut lower parts of the tree
+  if (f2 > 0.0)
+    value += f2 * exploreTree(j + 1, lower2, upper2, count2, xReduced);
+  return value * normalizationProba;
+}
+
+/* Get the CDF of the distribution */
+Scalar CombinationsDistribution::computeCDF(const Point & point) const
+{
+  if (point.getDimension() != k_) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=" << k_ << ", here dimension=" << point.getDimension();
+
+  Point sortedPoint(k_);
+  Scalar oldX = n_ - 1.0;
+  if (sortedPoint[k_ - 1] < -supportEpsilon_) return 0.0;
+  for (SignedInteger i = k_ - 1; i >= 0; --i)
+  {
+    const Scalar x = std::min(oldX, floor(point[i] + supportEpsilon_));
+    if (x < -supportEpsilon_) return 0.0;
+    sortedPoint[i] = x;
+    oldX = x - 1.0;
+  }
+  // Explore the tree describing the domain of summation
+  // k = 1
+  // lower = dimension_ (guard value to tell that the current interval is [0, xReduced[0]])
+  // upper = 0
+  // count = 1
+  return ((sortedPoint[0] + 1.0) * k_) / n_ * exploreTree(1, k_, 0, 1, sortedPoint);
+}
+
+/* Compute the scalar quantile of the 1D CombinationsDistribution distribution */
+Scalar CombinationsDistribution::computeScalarQuantile(const Scalar prob,
+    const Bool tail) const
+{
+  const UnsignedInteger i = static_cast< UnsignedInteger >(ceil(prob * (n_ - 1.0)));
+  return (tail ? n_ - 1.0 - i : i);
+} // computeScalarQuantile
+
+/* Compute the quantile of the CombinationsDistribution distribution */
+Point CombinationsDistribution::computeQuantile(const Scalar prob,
+    const Bool tail,
+    Scalar & marginalProb) const
+{
+  const Scalar p = (tail ? 1.0 - prob : prob);
+  if (p <= 0.0) return Point(k_, 0.0);
+  if (p >= 1.0) return Point(k_, n_);
+  UnsignedInteger iMin = 0;
+  UnsignedInteger iMax = n_;
+  while (iMax > iMin + 1)
+  {
+    const UnsignedInteger iMiddle = (iMax + iMin) / 2;
+    const Scalar cdfMiddle = computeCDF(Point(k_, iMiddle));
+    if (cdfMiddle < p)
+    {
+      iMin = iMiddle;
+    }
+    else
+    {
+      iMax = iMiddle;
+    }
+  } // while
+  marginalProb = computeScalarQuantile(prob, tail);
+  return Point(k_, iMax);
+} // computeQuantile
+
+/* Get the i-th marginal distribution */
+Distribution CombinationsDistribution::getMarginal(const UnsignedInteger i) const
+{
+  if (i >= k_) throw InvalidArgumentException(HERE) << "The index of a marginal distribution must be in the range [0, dim-1]";
+  Sample support(n_ - k_ + 1, 1);
+  // First, compute the probabilities on a log scale
+  Point probabilities(n_ - k_ + 1, -SpecFunc::LogBinomialCoefficient(k_, n_));
+  for (UnsignedInteger x = i; x <= n_ - k_ + i; ++x)
+    {
+      support(x - i, 0) = x;
+      probabilities[x - i] += SpecFunc::LogBinomialCoefficient(x, i) + SpecFunc::LogBinomialCoefficient(n_ - 1 - x, k_ - 1 - i);
+    }
+  // Then, go back to the [0, 1] interval
+  for (UnsignedInteger j = 0; j <= n_ - k_; ++j)
+    probabilities[j] = SpecFunc::Clip01(std::exp(probabilities[j]));
+  UserDefined marginal(support, probabilities);
+  marginal.setDescription(Description(1, getDescription()[i]));
+  return marginal;
+}
+
+/* Get the distribution of the marginal distribution corresponding to indices dimensions */
+Distribution CombinationsDistribution::getMarginal(const Indices & indices) const
+{
+  const UnsignedInteger dimension = getDimension();
+  if (!indices.check(dimension)) throw InvalidArgumentException(HERE) << "The indices of a marginal distribution must be in the range [0, dim-1] and must be different";
+  // Special case for dimension 1
+  if (dimension == 1) return clone();
+  // Special case for indices of length 1
+  if (indices.getSize() == 1) return getMarginal(indices[0]);
+  // General case
+  const UnsignedInteger outputDimension = indices.getSize();
+  const Sample support(getSupport().getMarginal(indices));
+  Point probabilities(support.getSize());
+  for (UnsignedInteger j = 0; j < support.getSize(); ++j)
+    {
+      const UnsignedInteger j0 = indices[0];
+      const UnsignedInteger x0 = static_cast<UnsignedInteger>(support(j, 0));
+      const UnsignedInteger je = indices[outputDimension - 1];
+      const UnsignedInteger xe = static_cast<UnsignedInteger>(support(j, outputDimension - 1));
+      probabilities[j] += SpecFunc::LogBinomialCoefficient(x0, j0) + SpecFunc::LogBinomialCoefficient(n_ - xe - 1, k_ - je - 1);
+      for (UnsignedInteger i = 1; i < outputDimension - 1; ++i)
+        probabilities[j] += SpecFunc::LogBinomialCoefficient(xe, je);
+      probabilities[j] = SpecFunc::Clip01(std::exp(probabilities[j]));
+    }
+  UserDefined marginal(support, probabilities);
+  marginal.setDescription(getDescription().select(indices));
+  return marginal;
+}
+
+/* Get the support of a discrete distribution that intersect a given interval */
+Sample CombinationsDistribution::getSupport(const Interval & interval) const
+{
+  if (interval.getDimension() != getDimension()) throw InvalidArgumentException(HERE) << "Error: the given interval has a dimension that does not match the distribution dimension.";
+  // Convert int values into float
+  const IndicesCollection intResult(Combinations(k_, n_).generate());
+  const UnsignedInteger size = intResult.getSize();
+  if (size == 0) return Sample();
+  const Interval inter(interval.intersect(range_));
+  // Common case: get the full support
+  if (inter == range_)
+    {
+      Sample result(size, dimension_);
+      for (UnsignedInteger i = 0; i < size; ++i)
+        for (UnsignedInteger j = 0; j < dimension_; ++j)
+            result(i, j) = intResult(i, j);
+      return result;
+    }
+  Sample result(0, dimension_);
+  for (UnsignedInteger i = 0; i < size; ++i)
+    {
+      Point point(dimension_);
+      for (UnsignedInteger j = 0; j < dimension_; ++j)
+        point[j] = intResult(i, j);
+      if (inter.contains(point))
+        result.add(point);
+    }
+  return result;
+}
+
+/* Compute the mean of the distribution */
+void CombinationsDistribution::computeMean() const
+{
+  mean_ = Point(dimension_);
+  for (UnsignedInteger i = 0; i < dimension_; ++i)
+    mean_[i] = getMarginal(i).getMean()[0];
+  isAlreadyComputedMean_ = true;
+}
+
+/* Compute the covariance of the distribution */
+void CombinationsDistribution::computeCovariance() const
+{
+  covariance_ = CovarianceMatrix(dimension_);
+  for (UnsignedInteger j = 0; j < dimension_; ++j)
+    {
+      for (UnsignedInteger i = 0; i < j; ++i)
+        {
+          const CovarianceMatrix covIJ(getMarginal({i, j}).getCovariance());
+          covariance_(i, j) = covIJ(0, 1);
+          covariance_(i, i) = covIJ(0, 0);
+        }
+    }
+  isAlreadyComputedCovariance_ = true;
+}
+
+/* Parameters value and description accessor */
+CombinationsDistribution::PointWithDescriptionCollection CombinationsDistribution::getParametersCollection() const
+{
+  const UnsignedInteger dimension = getDimension();
+  PointWithDescriptionCollection parameters((dimension == 1 ? 1 : dimension + 1));
+  for (UnsignedInteger i = 0; i < dimension; ++i)
+  {
+    PointWithDescription point(1);
+    point[0] = n_;
+    Description description(1);
+    description[0] = "n";
+    point.setDescription(description);
+    point.setName(getDescription()[i]);
+    parameters[i] = point;
+  }
+  if (dimension > 1)
+  {
+    PointWithDescription point(2);
+    Description description = {"k", "n"};
+    point[0] = k_;
+    point[1] = n_;
+    point.setDescription(description);
+    point.setName("dependence");
+    parameters[dimension] = point;
+  }
+  return parameters;
+}
+
+/* K accessor */
+void CombinationsDistribution::setK(const UnsignedInteger k)
+{
+  if (k == 0) throw InvalidArgumentException(HERE) << "Error: k must be > 0.";
+  if (k > n_) throw InvalidArgumentException(HERE) << "Error: k must be less or equal to n, here k=" << k << " and n=" << n_;
+  if (k != k_)
+  {
+    k_ = k;
+    logPDFValue_ = SpecFunc::LogGamma(n_ - k_ + 1) - SpecFunc::LogGamma(n_ + 1) + SpecFunc::LogGamma(k_ + 1);
+    setDimension(k);
+    isAlreadyComputedMean_ = false;
+    isAlreadyComputedCovariance_ = false;
+    isAlreadyCreatedGeneratingFunction_ = false;
+  }
+}
+
+/* K accessor */
+UnsignedInteger CombinationsDistribution::getK() const
+{
+  return k_;
+}
+
+/* N accessor */
+void CombinationsDistribution::setN(const UnsignedInteger n)
+{
+  if (n == 0) throw InvalidArgumentException(HERE) << "Error: n must be > 0.";
+  if (n < k_) throw InvalidArgumentException(HERE) << "Error: n must be greater or equal to k, here n=" << n << " and k=" << k_;
+  if (n != n_)
+  {
+    n_ = n;
+    logPDFValue_ = SpecFunc::LogGamma(n_ - k_ + 1) - SpecFunc::LogGamma(n_ + 1) + SpecFunc::LogGamma(k_ + 1);
+    isAlreadyComputedMean_ = false;
+    isAlreadyComputedCovariance_ = false;
+    computeRange();
+  }
+}
+
+UnsignedInteger CombinationsDistribution::getN() const
+{
+  return n_;
+}
+
+/* K/N accessor */
+void CombinationsDistribution::setKN(const UnsignedInteger k,
+                                      const UnsignedInteger n)
+{
+  if (k == 0) throw InvalidArgumentException(HERE) << "Error: k must be > 0.";
+  if (n == 0) throw InvalidArgumentException(HERE) << "Error: n must be > 0.";
+  if (k > n) throw InvalidArgumentException(HERE) << "Error: k must be less or equal to n, here k=" << k << " and n=" << n;
+  k_ = k;
+  setDimension(k);
+  n_ = n;
+  logPDFValue_ = SpecFunc::LogGamma(n_ - k_ + 1) - SpecFunc::LogGamma(n_ + 1) + SpecFunc::LogGamma(k_ + 1);
+  isAlreadyComputedMean_ = false;
+  isAlreadyComputedCovariance_ = false;
+  computeRange();
+}
+
+/* Method save() stores the object through the StorageManager */
+void CombinationsDistribution::save(Advocate & adv) const
+{
+  DiscreteDistribution::save(adv);
+  adv.saveAttribute( "k_", k_ );
+  adv.saveAttribute( "n_", n_ );
+  adv.saveAttribute( "logPDFValue_", logPDFValue_ );
+}
+
+/* Method load() reloads the object from the StorageManager */
+void CombinationsDistribution::load(Advocate & adv)
+{
+  DiscreteDistribution::load(adv);
+  adv.loadAttribute( "k_", k_ );
+  adv.loadAttribute( "n_", n_ );
+  adv.loadAttribute( "logPDFValue_", logPDFValue_ );
+  computeRange();
+}
+
+END_NAMESPACE_OPENTURNS
diff --git a/lib/src/Uncertainty/Distribution/openturns/CombinationsDistribution.hxx b/lib/src/Uncertainty/Distribution/openturns/CombinationsDistribution.hxx
new file mode 100644
index 0000000000..5fa093a55c
--- /dev/null
+++ b/lib/src/Uncertainty/Distribution/openturns/CombinationsDistribution.hxx
@@ -0,0 +1,161 @@
+//                                               -*- C++ -*-
+/**
+ *  @brief The CombinationsDistribution distribution
+ *
+ *  Copyright 2005-2024 Airbus-EDF-IMACS-ONERA-Phimeca
+ *
+ *  This library is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU Lesser General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This library is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU Lesser General Public License for more details.
+ *
+ *  You should have received a copy of the GNU Lesser General Public License
+ *  along with this library.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+#ifndef OPENTURNS_COMBINATIONSDISTRIBUTION_HXX
+#define OPENTURNS_COMBINATIONSDISTRIBUTION_HXX
+
+#include "openturns/OTprivate.hxx"
+#include "openturns/DiscreteDistribution.hxx"
+#include "openturns/Distribution.hxx"
+
+BEGIN_NAMESPACE_OPENTURNS
+
+/**
+ * @class CombinationsDistribution
+ *
+ * The CombinationsDistribution distribution.
+ */
+class OT_API CombinationsDistribution
+  : public DiscreteDistribution
+{
+  CLASSNAME
+public:
+
+  /** Default constructor */
+  CombinationsDistribution();
+
+
+  /** Parameters constructor */
+  CombinationsDistribution(const UnsignedInteger k,
+                            const UnsignedInteger n);
+
+
+  /** Comparison operator */
+  using DiscreteDistribution::operator ==;
+  Bool operator ==(const CombinationsDistribution & other) const;
+protected:
+  Bool equals(const DistributionImplementation & other) const override;
+public:
+
+  /** String converter */
+  String __repr__() const override;
+  String __str__(const String & offset = "") const override;
+
+
+
+  /* Interface inherited from Distribution */
+
+  /** Virtual constructor */
+  CombinationsDistribution * clone() const override;
+
+  /** Get one realization of the distribution */
+  Point getRealization() const override;
+
+  /** Get the PDF of the distribution */
+  using DiscreteDistribution::computePDF;
+  Scalar computePDF(const Point & point) const override;
+  using DiscreteDistribution::computeLogPDF;
+  Scalar computeLogPDF(const Point & point) const override;
+
+  /** Get the CDF of the distribution */
+  using DiscreteDistribution::computeCDF;
+  Scalar computeCDF(const Point & point) const override;
+
+  /** Get the quantile of the distribution */
+  using DiscreteDistribution::computeQuantile;
+  Point computeQuantile(const Scalar prob,
+                        const Bool tail,
+                        Scalar & marginalProb) const override;
+
+  /** Get the i-th marginal distribution */
+  using DiscreteDistribution::getMarginal;
+  Distribution getMarginal(const UnsignedInteger i) const override;
+
+  /** Get the distribution of the marginal distribution corresponding to indices dimensions */
+  Distribution getMarginal(const Indices & indices) const override;
+
+  /** Get the support of a discrete distribution that intersect a given interval */
+  using DistributionImplementation::getSupport;
+  Sample getSupport(const Interval & interval) const override;
+
+  /** Parameters value and description accessor */
+  PointWithDescriptionCollection getParametersCollection() const override;
+
+  /* Interface specific to CombinationsDistribution */
+
+  /** K accessor */
+  void setK(const UnsignedInteger k);
+  UnsignedInteger getK() const;
+
+  /** N accessor */
+  void setN(const UnsignedInteger n);
+  UnsignedInteger getN() const;
+
+private:
+  /** K/N accessor */
+  void setKN(const UnsignedInteger k,
+             const UnsignedInteger n);
+
+public:
+  /** Method save() stores the object through the StorageManager */
+  void save(Advocate & adv) const override;
+
+  /** Method load() reloads the object from the StorageManager */
+  void load(Advocate & adv) override;
+
+protected:
+
+
+private:
+
+  /** Compute the numerical range of the distribution given the parameters values */
+  void computeRange() override;
+
+  /** Quantile computation for dimension=1 */
+  Scalar computeScalarQuantile(const Scalar prob, const Bool tail = false) const override;
+
+  /** Compute the mean of the distribution */
+  void computeMean() const override;
+
+  /** Compute the covariance of the distribution */
+  void computeCovariance() const override;
+
+  /** Helper to compute the CDF */
+  Scalar exploreTree(const Scalar j,
+                     const UnsignedInteger lower,
+                     const UnsignedInteger upper,
+                     const UnsignedInteger count,
+                     const Point & xReduced) const;
+
+  /** Size of the permutations */
+  UnsignedInteger k_;
+
+  /** Size of the base set */
+  UnsignedInteger n_;
+
+  /** Log PDF value */
+  Scalar logPDFValue_;
+
+}; /* class CombinationsDistribution */
+
+
+END_NAMESPACE_OPENTURNS
+
+#endif /* OPENTURNS_KPERMUTATIONSDISTRIBUTION_HXX */
diff --git a/lib/src/Uncertainty/Distribution/openturns/CombinationsDistributions.hxx b/lib/src/Uncertainty/Distribution/openturns/CombinationsDistributions.hxx
new file mode 100644
index 0000000000..2555eae45a
--- /dev/null
+++ b/lib/src/Uncertainty/Distribution/openturns/CombinationsDistributions.hxx
@@ -0,0 +1,153 @@
+//                                               -*- C++ -*-
+/**
+ *  @brief The CombinationsDistribution distribution
+ *
+ *  Copyright 2005-2024 Airbus-EDF-IMACS-ONERA-Phimeca
+ *
+ *  This library is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU Lesser General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This library is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU Lesser General Public License for more details.
+ *
+ *  You should have received a copy of the GNU Lesser General Public License
+ *  along with this library.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+#ifndef OPENTURNS_COMBINATIONSDISTRIBUTION_HXX
+#define OPENTURNS_COMBINATIONSDISTRIBUTION_HXX
+
+#include "openturns/OTprivate.hxx"
+#include "openturns/DiscreteDistribution.hxx"
+
+BEGIN_NAMESPACE_OPENTURNS
+
+/**
+ * @class CombinationsDistribution
+ *
+ * The CombinationsDistribution distribution.
+ */
+class OT_API CombinationsDistribution
+  : public DiscreteDistribution
+{
+  CLASSNAME
+public:
+
+  /** Default constructor */
+  CombinationsDistribution();
+
+
+  /** Parameters constructor */
+  CombinationsDistribution(const UnsignedInteger k,
+                            const UnsignedInteger n);
+
+
+  /** Comparison operator */
+  using DiscreteDistribution::operator ==;
+  Bool operator ==(const CombinationsDistribution & other) const;
+protected:
+  Bool equals(const DistributionImplementation & other) const override;
+public:
+
+  /** String converter */
+  String __repr__() const override;
+  String __str__(const String & offset = "") const override;
+
+
+
+  /* Interface inherited from Distribution */
+
+  /** Virtual constructor */
+  CombinationsDistribution * clone() const override;
+
+  /** Get one realization of the distribution */
+  Point getRealization() const override;
+
+  /** Get the PDF of the distribution */
+  using DiscreteDistribution::computePDF;
+  Scalar computePDF(const Point & point) const override;
+  using DiscreteDistribution::computeLogPDF;
+  Scalar computeLogPDF(const Point & point) const override;
+
+  /** Get the CDF of the distribution */
+  using DiscreteDistribution::computeCDF;
+  Scalar computeCDF(const Point & point) const override;
+
+  /** Get the quantile of the distribution */
+  using DiscreteDistribution::computeQuantile;
+  Point computeQuantile(const Scalar prob,
+                        const Bool tail,
+                        Scalar & marginalProb) const override;
+
+  /** Get the i-th marginal distribution */
+  using DiscreteDistribution::getMarginal;
+  Distribution getMarginal(const UnsignedInteger i) const override;
+
+  /** Get the distribution of the marginal distribution corresponding to indices dimensions */
+  Distribution getMarginal(const Indices & indices) const override;
+
+  /** Get the support of a discrete distribution that intersect a given interval */
+  using DistributionImplementation::getSupport;
+  Sample getSupport(const Interval & interval) const override;
+
+  /** Parameters value and description accessor */
+  PointWithDescriptionCollection getParametersCollection() const override;
+
+  /* Interface specific to CombinationsDistribution */
+
+  /** K accessor */
+  void setK(const UnsignedInteger k);
+  UnsignedInteger getK() const;
+
+  /** N accessor */
+  void setN(const UnsignedInteger n);
+  UnsignedInteger getN() const;
+
+private:
+  /** K/N accessor */
+  void setKN(const UnsignedInteger k,
+             const UnsignedInteger n);
+
+public:
+  /** Method save() stores the object through the StorageManager */
+  void save(Advocate & adv) const override;
+
+  /** Method load() reloads the object from the StorageManager */
+  void load(Advocate & adv) override;
+
+protected:
+
+
+private:
+
+  /** Compute the numerical range of the distribution given the parameters values */
+  void computeRange() override;
+
+  /** Quantile computation for dimension=1 */
+  Scalar computeScalarQuantile(const Scalar prob, const Bool tail = false) const override;
+
+  /** Compute the mean of the distribution */
+  void computeMean() const override;
+
+  /** Compute the covariance of the distribution */
+  void computeCovariance() const override;
+
+  /** Size of the permutations */
+  UnsignedInteger k_;
+
+  /** Size of the base set */
+  UnsignedInteger n_;
+
+  /** Log PDF value */
+  Scalar logPDFValue_;
+
+}; /* class CombinationsDistribution */
+
+
+END_NAMESPACE_OPENTURNS
+
+#endif /* OPENTURNS_KPERMUTATIONSDISTRIBUTION_HXX */
diff --git a/lib/src/Uncertainty/Distribution/openturns/OTDistribution.hxx b/lib/src/Uncertainty/Distribution/openturns/OTDistribution.hxx
index eaf18b0ac5..0f2eefc2c9 100644
--- a/lib/src/Uncertainty/Distribution/openturns/OTDistribution.hxx
+++ b/lib/src/Uncertainty/Distribution/openturns/OTDistribution.hxx
@@ -45,6 +45,7 @@
 #include "openturns/ChiFactory.hxx"
 #include "openturns/ClaytonCopula.hxx"
 #include "openturns/ClaytonCopulaFactory.hxx"
+#include "openturns/CombinationsDistribution.hxx"
 #include "openturns/BlockIndependentCopula.hxx"
 #include "openturns/JointDistribution.hxx"
 #include "openturns/CompositeDistribution.hxx"
diff --git a/lib/test/CMakeLists.txt b/lib/test/CMakeLists.txt
index ccf74e675f..f8a5b3e7b8 100644
--- a/lib/test/CMakeLists.txt
+++ b/lib/test/CMakeLists.txt
@@ -356,6 +356,7 @@ ot_check_test (ChiSquareFactory_std)
 ot_check_test (BlockIndependentCopula_std)
 ot_check_test (JointDistribution_std)
 ot_check_test (JointDistribution_large)
+ot_check_test (CombinationsDistribution_std IGNOREOUT)
 ot_check_test (CompositeDistribution_std)
 ot_check_test (ConditionalDistribution_std)
 ot_check_test (CumulativeDistributionNetwork_std)
diff --git a/lib/test/t_CombinationsDistribution_std.cxx b/lib/test/t_CombinationsDistribution_std.cxx
new file mode 100644
index 0000000000..71b6df5246
--- /dev/null
+++ b/lib/test/t_CombinationsDistribution_std.cxx
@@ -0,0 +1,114 @@
+//                                               -*- C++ -*-
+/**
+ *  @brief The test file of class CombinationsDistribution for standard methods
+ *
+ *  Copyright 2005-2024 Airbus-EDF-IMACS-ONERA-Phimeca
+ *
+ *  This library is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU Lesser General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This library is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU Lesser General Public License for more details.
+ *
+ *  You should have received a copy of the GNU Lesser General Public License
+ *  along with this library.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+#include "openturns/OT.hxx"
+#include "openturns/OTtestcode.hxx"
+
+using namespace OT;
+using namespace OT::Test;
+
+class TestObject : public CombinationsDistribution
+{
+public:
+  TestObject() : CombinationsDistribution(5, 12) {}
+  virtual ~TestObject() {}
+};
+
+
+int main(int, char *[])
+{
+  TESTPREAMBLE;
+  OStream fullprint(std::cout);
+  setRandomGenerator();
+  try
+  {
+    // Test basic functionnalities
+    checkClassWithClassName<TestObject>();
+
+    // Instantiate one distribution object
+    CombinationsDistribution distribution(5, 12);
+    fullprint << "Distribution " << distribution << std::endl;
+    std::cout << "Distribution " << distribution << std::endl;
+    UserDefined ref(distribution.getSupport());
+    // Is this distribution elliptical ?
+    fullprint << "Elliptical = " << (distribution.isElliptical() ? "true" : "false") << std::endl;
+
+    // Is this distribution continuous ?
+    fullprint << "Continuous = " << (distribution.isContinuous() ? "true" : "false") << std::endl;
+
+    // Test for realization of distribution
+    Point oneRealization = distribution.getRealization();
+    fullprint << "oneRealization=" << oneRealization << std::endl;
+
+    // Test for sampling
+    UnsignedInteger size = 10000;
+    Sample oneSample = distribution.getSample( size );
+    fullprint << "oneSample first=" << oneSample[0] << " last=" << oneSample[size - 1] << std::endl;
+    fullprint << "mean=" << oneSample.computeMean() << std::endl;
+    fullprint << "covariance=" << oneSample.computeCovariance() << std::endl;
+    // Define a point
+    Point point({1, 3, 6, 8, 10});
+    fullprint << "Point= " << point << std::endl;
+
+    // Show PDF and CDF of point
+    Scalar LPDF = distribution.computeLogPDF( point );
+    Scalar PDF = distribution.computePDF( point );
+    fullprint << "pdf     =" << PDF << std::endl;
+    fullprint << "pdf ref =" << ref.computePDF(point) << std::endl;
+    Scalar CDF = distribution.computeCDF( point );
+    fullprint << "cdf     =" << CDF << std::endl;
+    fullprint << "cdf ref =" << ref.computeCDF(point) << std::endl;
+    Scalar CCDF = distribution.computeComplementaryCDF( point );
+    fullprint << "ccdf     =" << CCDF << std::endl;
+    fullprint << "ccdf ref =" << ref.computeComplementaryCDF(point) << std::endl;
+    //    Scalar Survival = distribution.computeSurvivalFunction( point );
+    //    fullprint << "survival=" << Survival << std::endl;
+    Point quantile = distribution.computeQuantile( 0.95 );
+    fullprint << "quantile     =" << quantile << std::endl;
+    fullprint << "quantile ref =" << ref.computeQuantile(0.95) << std::endl;
+    fullprint << "cdf(quantile)=" << distribution.computeCDF(quantile) << std::endl;
+    fullprint << "entropy     =" << distribution.computeEntropy() << std::endl;
+    fullprint << "entropy ref =" << ref.computeEntropy() << std::endl;
+    fullprint << "entropy (MC)=" << -distribution.computeLogPDF(distribution.getSample(1000000)).computeMean()[0] << std::endl;
+    Point mean = distribution.getMean();
+    fullprint << "mean     =" << mean << std::endl;
+    fullprint << "mean ref =" << ref.getMean() << std::endl;
+    CovarianceMatrix covariance = distribution.getCovariance();
+    fullprint << "covariance     =" << covariance << std::endl;
+    fullprint << "covariance ref =" << ref.getCovariance() << std::endl;
+    CovarianceMatrix correlation = distribution.getCorrelation();
+    fullprint << "correlation     =" << correlation << std::endl;
+    fullprint << "correlation ref =" << ref.getCorrelation() << std::endl;
+    //    CovarianceMatrix spearman = distribution.getSpearmanCorrelation();
+    //    fullprint << "spearman=" << spearman << std::endl;
+    //    CovarianceMatrix kendall = distribution.getKendallTau();
+    //    fullprint << "kendall=" << kendall << std::endl;
+    CombinationsDistribution::PointWithDescriptionCollection parameters = distribution.getParametersCollection();
+    fullprint << "parameters=" << parameters << std::endl;
+  }
+  catch (TestFailed & ex)
+  {
+    std::cerr << ex << std::endl;
+    return ExitCode::Error;
+  }
+
+
+  return ExitCode::Success;
+}
diff --git a/python/doc/pyplots/CombinationsDistribution.py b/python/doc/pyplots/CombinationsDistribution.py
new file mode 100644
index 0000000000..1f731509c0
--- /dev/null
+++ b/python/doc/pyplots/CombinationsDistribution.py
@@ -0,0 +1,24 @@
+import openturns as ot
+import openturns.experimental as otexp
+from openturns.viewer import View
+
+grid = ot.GridLayout(1, 2)
+pdf_2d = ot.Graph("CombinationsDistribution, PDF", "x1", "x2", True)
+cdf_2d = ot.Graph("CombinationsDistribution, CDF", "x1", "x2", True)
+
+distribution_2d = otexp.CombinationsDistribution(2, 10)
+
+xMin = [-0.1] * 2
+xMax = [10.1] * 2
+
+pdf_2d.add(distribution_2d.drawPDF(xMin, xMax))
+pdf_2d.setLegends([""])
+cdf_2d.add(distribution_2d.drawCDF(xMin, xMax, [101] * 2))
+cdf_2d.setLegends([""])
+grid.setGraph(0, 0, pdf_2d)
+grid.setGraph(0, 1, cdf_2d)
+grid.setTitle("CombinationsDistribution")
+grid.setLegendPosition("upper right")
+v = View(grid)
+fig = v.getFigure()
+fig.axes[1].legend(loc="best")
diff --git a/python/doc/user_manual/probabilistic_modelling.rst b/python/doc/user_manual/probabilistic_modelling.rst
index c1e15eeecc..27e7387602 100644
--- a/python/doc/user_manual/probabilistic_modelling.rst
+++ b/python/doc/user_manual/probabilistic_modelling.rst
@@ -111,6 +111,11 @@ Discrete parametric distributions
 
     Bernoulli
     Binomial
+
+    :template: classWithPlot.rst_t
+    experimental.CombinationsDistribution
+
+    :template: Distribution.rst_t
     Dirac
     Geometric
     Hypergeometric
diff --git a/python/src/CMakeLists.txt b/python/src/CMakeLists.txt
index c7a87967f0..e1070831f8 100644
--- a/python/src/CMakeLists.txt
+++ b/python/src/CMakeLists.txt
@@ -632,6 +632,7 @@ ot_add_python_module(dist_bundle1 dist_bundle1_module.i
                       ChiSquareFactory.i ChiSquareFactory_doc.i.in
                       Chi.i Chi_doc.i.in
                       ChiFactory.i ChiFactory_doc.i.in
+                      CombinationsDistribution.i CombinationsDistribution_doc.i.in
                       CompositeDistribution.i CompositeDistribution_doc.i.in
                       Dirac.i Dirac_doc.i.in
                       DiracFactory.i DiracFactory_doc.i.in
diff --git a/python/src/CombinationsDistribution.i b/python/src/CombinationsDistribution.i
new file mode 100644
index 0000000000..0d595ebec4
--- /dev/null
+++ b/python/src/CombinationsDistribution.i
@@ -0,0 +1,11 @@
+// SWIG file CombinationsDistribution.i
+
+%{
+#include "openturns/CombinationsDistribution.hxx"
+%}
+
+%include CombinationsDistribution_doc.i
+
+%copyctor OT::CombinationsDistribution;
+
+%include openturns/CombinationsDistribution.hxx
diff --git a/python/src/CombinationsDistribution_doc.i.in b/python/src/CombinationsDistribution_doc.i.in
new file mode 100644
index 0000000000..84f07dc2f7
--- /dev/null
+++ b/python/src/CombinationsDistribution_doc.i.in
@@ -0,0 +1,85 @@
+%feature("docstring") OT::CombinationsDistribution
+"Combinations distribution.
+
+Parameters
+----------
+k : int, :math:`k > 0`
+n : int, :math:`n > 0`
+
+See also
+--------
+:class:`~openturns.Combinations`
+
+Notes
+-----
+:class:`~openturns.experimental.CombinationsDistribution` is the discrete uniform
+distribution on the set of strictly increasing functions :math:`(i_0, \hdots, i_{k_1})`
+from :math:`\{0, \hdots, k-1\}` into :math:`\{0, \hdots, n-1\}`.
+Its probability distribution function is defined as:
+
+.. math::
+
+    \Prob{\vect{X} = (i_0, \hdots, i_{k-1})} = \frac{1}{\binom{n}{k}}
+
+Its first moments are:
+
+.. math::
+    :nowrap:
+
+    \begin{eqnarray*}
+        \Expect{\vect{X}} & = & \frac{n - 1}{2}\\
+        \Cov{\vect{X}} & = & \left\{
+                             \begin{array}{ll}
+                                 \displaystyle \frac{n^2-1}{12} \\
+                                 \displaystyle -\frac{1 + n}{12}
+                             \end{array}
+                             \right.
+    \end{eqnarray*}
+
+Examples
+--------
+Create a distribution:
+
+>>> import openturns as ot
+>>> import openturns.experimental as otexp
+>>> distribution = otexp.CombinationsDistribution(6, 8)
+
+Draw a sample:
+
+>>> sample = distribution.getSample(10)"
+
+// ---------------------------------------------------------------------
+
+%feature("docstring") OT::CombinationsDistribution::getN
+"Accessor to the parameter :math:`n`.
+
+Returns
+-------
+n : int"
+
+// ---------------------------------------------------------------------
+
+%feature("docstring") OT::CombinationsDistribution::getK
+"Accessor to the parameter :math:`k`.
+
+Returns
+-------
+k : int"
+
+// ---------------------------------------------------------------------
+
+%feature("docstring") OT::CombinationsDistribution::setN
+"Accessor to the parameter :math:`n`.
+
+Parameters
+----------
+n : int, :math:`n > 0`"
+
+// ---------------------------------------------------------------------
+
+%feature("docstring") OT::CombinationsDistribution::setK
+"Accessor to the parameter :math:`k`.
+
+Parameters
+----------
+k : int, :math:`k > 0`"
diff --git a/python/src/dist_bundle1_module.i b/python/src/dist_bundle1_module.i
index 4297efa083..95bd5bae6a 100644
--- a/python/src/dist_bundle1_module.i
+++ b/python/src/dist_bundle1_module.i
@@ -41,6 +41,7 @@
 %include ChiSquare.i
 %include ChiSquareFactory.i
 %include CompositeDistribution.i
+%include CombinationsDistribution.i
 %include Dirac.i
 %include DiracFactory.i
 %include Dirichlet.i