sagemathgh-35939: Support for exponential and gamma distribution

    
I think that the exponential and gamma distributions are essential.
Instead of filing an issue I tried to edit the code. Please have a look
at the changes.


### 📚 Description

Support for exponential and gamma distribution from GSL.


### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. It should be `[x]` not `[x
]`. -->

- [x] The title is concise, informative, and self-explanatory.
- [x] The description explains in detail what this PR is about.
- [ ] I have linked a relevant issue or discussion.
- [ ] I have created tests covering the changes.
- [x] I have updated the documentation accordingly.

### ⌛ Dependencies

<!-- List all open PRs that this PR logically depends on
- sagemath#12345: short description why this is a dependency
- sagemath#34567: ...
-->

<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
    
URL: sagemath#35939
Reported by: Gerald Teschl
Reviewer(s): Matthias Köppe

src/sage/probability/probability_distribution.pyx

@@ -46,8 +46,7 @@ import sage.rings.real_double
 from sage.modules.free_module_element import vector
 
 # TODO: Add more distributions available in gsl
-# available but not currently wrapped are exponential, laplace, cauchy, landau, gamma,
-# gamma, beta logistic.
+# available but not currently wrapped are laplace, cauchy, landau, logistic.
 
 cdef enum:
     uniform
@@ -61,6 +60,8 @@ cdef enum:
     exppow
     weibull
     beta
+    exponential
+    gamma
 
 cdef class ProbabilityDistribution:
     r"""
@@ -501,6 +502,31 @@ cdef class RealDistribution(ProbabilityDistribution):
         sage: T.cum_distribution_function(1)
         1.0
 
+    The exponential distribution has one parameter ``mu``::
+
+        sage: mu = 2
+        sage: T = RealDistribution('exponential', mu)
+        sage: s = T.get_random_element()
+        sage: 0 <= s
+        True
+        sage: s.parent()
+        Real Double Field
+        sage: T.distribution_function(0)
+        0.5
+
+    The gamma distribution has two parameters ``a`` and ``b``::
+
+        sage: a = 2
+        sage: b = 2
+        sage: T = RealDistribution('gamma', [a, b])
+        sage: s = T.get_random_element()
+        sage: 0 <= s
+        True
+        sage: s.parent()
+        Real Double Field
+        sage: T.distribution_function(0)
+        0.0
+
     The weibull distribution has two parameters ``a`` and ``b``::
 
         sage: a = 1
@@ -682,6 +708,10 @@ cdef class RealDistribution(ProbabilityDistribution):
             result = gsl_ran_weibull(self.r, self.parameters[0], self.parameters[1])
         elif self.distribution_type == beta:
             result = gsl_ran_beta(self.r, self.parameters[0], self.parameters[1])
+        elif self.distribution_type == exponential:
+            result = gsl_ran_exponential(self.r, self.parameters[0])
+        elif self.distribution_type == gamma:
+            result = gsl_ran_gamma(self.r, self.parameters[0], self.parameters[1])
         else:
             raise TypeError("Not a supported probability distribution")
 
@@ -732,7 +762,6 @@ cdef class RealDistribution(ProbabilityDistribution):
             self.parameters[1] = float(parameters[1])
             self.distribution_type = pareto
         elif name == 'rayleigh':
-            self.distribution_type = rayleigh
             try:
                 float(parameters)
             except Exception:
@@ -813,6 +842,25 @@ cdef class RealDistribution(ProbabilityDistribution):
             self.parameters[0] = float(parameters[0])
             self.parameters[1] = float(parameters[1])
             self.distribution_type = beta
+        elif name == 'exponential':
+            try:
+                float(parameters)
+            except Exception:
+                raise TypeError("exponential distribution requires parameter mu coercible to float")
+            self.parameters = <double*>sig_malloc(sizeof(double))
+            self.parameters[0] = float(parameters)
+            self.distribution_type = exponential
+        elif name == 'gamma':
+            if len(parameters) != 2:
+                raise TypeError("gamma distribution requires two real parameters")
+            try:
+                map(float, parameters)
+            except Exception:
+                raise TypeError("gamma distribution requires real parameters")
+            self.parameters = <double *>sig_malloc(sizeof(double)*2)
+            self.parameters[0] = float(parameters[0])
+            self.parameters[1] = float(parameters[1])
+            self.distribution_type = gamma
         else:
             raise TypeError("Not a supported probability distribution")
 
@@ -878,6 +926,10 @@ cdef class RealDistribution(ProbabilityDistribution):
             return sage.rings.real_double.RDF(gsl_ran_weibull_pdf(x, self.parameters[0], self.parameters[1]))
         elif self.distribution_type == beta:
             return sage.rings.real_double.RDF(gsl_ran_beta_pdf(x, self.parameters[0], self.parameters[1]))
+        elif self.distribution_type == exponential:
+            return sage.rings.real_double.RDF(gsl_ran_exponential_pdf(x, self.parameters[0]))
+        elif self.distribution_type == gamma:
+            return sage.rings.real_double.RDF(gsl_ran_gamma_pdf(x, self.parameters[0], self.parameters[1]))
         else:
             raise TypeError("Not a supported probability distribution")
 
@@ -914,6 +966,10 @@ cdef class RealDistribution(ProbabilityDistribution):
             return sage.rings.real_double.RDF(gsl_cdf_weibull_P(x, self.parameters[0], self.parameters[1]))
         elif self.distribution_type == beta:
             return sage.rings.real_double.RDF(gsl_cdf_beta_P(x, self.parameters[0], self.parameters[1]))
+        elif self.distribution_type == exponential:
+            return sage.rings.real_double.RDF(gsl_cdf_exponential_P(x, self.parameters[0]))
+        elif self.distribution_type == gamma:
+            return sage.rings.real_double.RDF(gsl_cdf_gamma_P(x, self.parameters[0], self.parameters[1]))
         else:
             raise TypeError("Not a supported probability distribution")
 
@@ -951,6 +1007,10 @@ cdef class RealDistribution(ProbabilityDistribution):
             return sage.rings.real_double.RDF(gsl_cdf_weibull_Pinv(x, self.parameters[0], self.parameters[1]))
         elif self.distribution_type == beta:
             return sage.rings.real_double.RDF(gsl_cdf_beta_Pinv(x, self.parameters[0], self.parameters[1]))
+        elif self.distribution_type == exponential:
+            return sage.rings.real_double.RDF(gsl_cdf_exponential_Pinv(x, self.parameters[0]))
+        elif self.distribution_type == gamma:
+            return sage.rings.real_double.RDF(gsl_cdf_gamma_Pinv(x, self.parameters[0], self.parameters[1]))
         else:
             raise TypeError("Not a supported probability distribution")