Coverage for /usr/lib/python3/dist-packages/sympy/functions/special/beta_functions.py: 33%
113 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from sympy.core import S
2from sympy.core.function import Function, ArgumentIndexError
3from sympy.core.symbol import Dummy
4from sympy.functions.special.gamma_functions import gamma, digamma
5from sympy.functions.combinatorial.numbers import catalan
6from sympy.functions.elementary.complexes import conjugate
8# See mpmath #569 and SymPy #20569
9def betainc_mpmath_fix(a, b, x1, x2, reg=0):
10 from mpmath import betainc, mpf
11 if x1 == x2:
12 return mpf(0)
13 else:
14 return betainc(a, b, x1, x2, reg)
16###############################################################################
17############################ COMPLETE BETA FUNCTION ##########################
18###############################################################################
20class beta(Function):
21 r"""
22 The beta integral is called the Eulerian integral of the first kind by
23 Legendre:
25 .. math::
26 \mathrm{B}(x,y) \int^{1}_{0} t^{x-1} (1-t)^{y-1} \mathrm{d}t.
28 Explanation
29 ===========
31 The Beta function or Euler's first integral is closely associated
32 with the gamma function. The Beta function is often used in probability
33 theory and mathematical statistics. It satisfies properties like:
35 .. math::
36 \mathrm{B}(a,1) = \frac{1}{a} \\
37 \mathrm{B}(a,b) = \mathrm{B}(b,a) \\
38 \mathrm{B}(a,b) = \frac{\Gamma(a) \Gamma(b)}{\Gamma(a+b)}
40 Therefore for integral values of $a$ and $b$:
42 .. math::
43 \mathrm{B} = \frac{(a-1)! (b-1)!}{(a+b-1)!}
45 A special case of the Beta function when `x = y` is the
46 Central Beta function. It satisfies properties like:
48 .. math::
49 \mathrm{B}(x) = 2^{1 - 2x}\mathrm{B}(x, \frac{1}{2})
50 \mathrm{B}(x) = 2^{1 - 2x} cos(\pi x) \mathrm{B}(\frac{1}{2} - x, x)
51 \mathrm{B}(x) = \int_{0}^{1} \frac{t^x}{(1 + t)^{2x}} dt
52 \mathrm{B}(x) = \frac{2}{x} \prod_{n = 1}^{\infty} \frac{n(n + 2x)}{(n + x)^2}
54 Examples
55 ========
57 >>> from sympy import I, pi
58 >>> from sympy.abc import x, y
60 The Beta function obeys the mirror symmetry:
62 >>> from sympy import beta, conjugate
63 >>> conjugate(beta(x, y))
64 beta(conjugate(x), conjugate(y))
66 Differentiation with respect to both $x$ and $y$ is supported:
68 >>> from sympy import beta, diff
69 >>> diff(beta(x, y), x)
70 (polygamma(0, x) - polygamma(0, x + y))*beta(x, y)
72 >>> diff(beta(x, y), y)
73 (polygamma(0, y) - polygamma(0, x + y))*beta(x, y)
75 >>> diff(beta(x), x)
76 2*(polygamma(0, x) - polygamma(0, 2*x))*beta(x, x)
78 We can numerically evaluate the Beta function to
79 arbitrary precision for any complex numbers x and y:
81 >>> from sympy import beta
82 >>> beta(pi).evalf(40)
83 0.02671848900111377452242355235388489324562
85 >>> beta(1 + I).evalf(20)
86 -0.2112723729365330143 - 0.7655283165378005676*I
88 See Also
89 ========
91 gamma: Gamma function.
92 uppergamma: Upper incomplete gamma function.
93 lowergamma: Lower incomplete gamma function.
94 polygamma: Polygamma function.
95 loggamma: Log Gamma function.
96 digamma: Digamma function.
97 trigamma: Trigamma function.
99 References
100 ==========
102 .. [1] https://en.wikipedia.org/wiki/Beta_function
103 .. [2] https://mathworld.wolfram.com/BetaFunction.html
104 .. [3] https://dlmf.nist.gov/5.12
106 """
107 unbranched = True
109 def fdiff(self, argindex):
110 x, y = self.args
111 if argindex == 1:
112 # Diff wrt x
113 return beta(x, y)*(digamma(x) - digamma(x + y))
114 elif argindex == 2:
115 # Diff wrt y
116 return beta(x, y)*(digamma(y) - digamma(x + y))
117 else:
118 raise ArgumentIndexError(self, argindex)
120 @classmethod
121 def eval(cls, x, y=None):
122 if y is None:
123 return beta(x, x)
124 if x.is_Number and y.is_Number:
125 return beta(x, y, evaluate=False).doit()
127 def doit(self, **hints):
128 x = xold = self.args[0]
129 # Deal with unevaluated single argument beta
130 single_argument = len(self.args) == 1
131 y = yold = self.args[0] if single_argument else self.args[1]
132 if hints.get('deep', True):
133 x = x.doit(**hints)
134 y = y.doit(**hints)
135 if y.is_zero or x.is_zero:
136 return S.ComplexInfinity
137 if y is S.One:
138 return 1/x
139 if x is S.One:
140 return 1/y
141 if y == x + 1:
142 return 1/(x*y*catalan(x))
143 s = x + y
144 if (s.is_integer and s.is_negative and x.is_integer is False and
145 y.is_integer is False):
146 return S.Zero
147 if x == xold and y == yold and not single_argument:
148 return self
149 return beta(x, y)
151 def _eval_expand_func(self, **hints):
152 x, y = self.args
153 return gamma(x)*gamma(y) / gamma(x + y)
155 def _eval_is_real(self):
156 return self.args[0].is_real and self.args[1].is_real
158 def _eval_conjugate(self):
159 return self.func(self.args[0].conjugate(), self.args[1].conjugate())
161 def _eval_rewrite_as_gamma(self, x, y, piecewise=True, **kwargs):
162 return self._eval_expand_func(**kwargs)
164 def _eval_rewrite_as_Integral(self, x, y, **kwargs):
165 from sympy.integrals.integrals import Integral
166 t = Dummy('t')
167 return Integral(t**(x - 1)*(1 - t)**(y - 1), (t, 0, 1))
169###############################################################################
170########################## INCOMPLETE BETA FUNCTION ###########################
171###############################################################################
173class betainc(Function):
174 r"""
175 The Generalized Incomplete Beta function is defined as
177 .. math::
178 \mathrm{B}_{(x_1, x_2)}(a, b) = \int_{x_1}^{x_2} t^{a - 1} (1 - t)^{b - 1} dt
180 The Incomplete Beta function is a special case
181 of the Generalized Incomplete Beta function :
183 .. math:: \mathrm{B}_z (a, b) = \mathrm{B}_{(0, z)}(a, b)
185 The Incomplete Beta function satisfies :
187 .. math:: \mathrm{B}_z (a, b) = (-1)^a \mathrm{B}_{\frac{z}{z - 1}} (a, 1 - a - b)
189 The Beta function is a special case of the Incomplete Beta function :
191 .. math:: \mathrm{B}(a, b) = \mathrm{B}_{1}(a, b)
193 Examples
194 ========
196 >>> from sympy import betainc, symbols, conjugate
197 >>> a, b, x, x1, x2 = symbols('a b x x1 x2')
199 The Generalized Incomplete Beta function is given by:
201 >>> betainc(a, b, x1, x2)
202 betainc(a, b, x1, x2)
204 The Incomplete Beta function can be obtained as follows:
206 >>> betainc(a, b, 0, x)
207 betainc(a, b, 0, x)
209 The Incomplete Beta function obeys the mirror symmetry:
211 >>> conjugate(betainc(a, b, x1, x2))
212 betainc(conjugate(a), conjugate(b), conjugate(x1), conjugate(x2))
214 We can numerically evaluate the Incomplete Beta function to
215 arbitrary precision for any complex numbers a, b, x1 and x2:
217 >>> from sympy import betainc, I
218 >>> betainc(2, 3, 4, 5).evalf(10)
219 56.08333333
220 >>> betainc(0.75, 1 - 4*I, 0, 2 + 3*I).evalf(25)
221 0.2241657956955709603655887 + 0.3619619242700451992411724*I
223 The Generalized Incomplete Beta function can be expressed
224 in terms of the Generalized Hypergeometric function.
226 >>> from sympy import hyper
227 >>> betainc(a, b, x1, x2).rewrite(hyper)
228 (-x1**a*hyper((a, 1 - b), (a + 1,), x1) + x2**a*hyper((a, 1 - b), (a + 1,), x2))/a
230 See Also
231 ========
233 beta: Beta function
234 hyper: Generalized Hypergeometric function
236 References
237 ==========
239 .. [1] https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function
240 .. [2] https://dlmf.nist.gov/8.17
241 .. [3] https://functions.wolfram.com/GammaBetaErf/Beta4/
242 .. [4] https://functions.wolfram.com/GammaBetaErf/BetaRegularized4/02/
244 """
245 nargs = 4
246 unbranched = True
248 def fdiff(self, argindex):
249 a, b, x1, x2 = self.args
250 if argindex == 3:
251 # Diff wrt x1
252 return -(1 - x1)**(b - 1)*x1**(a - 1)
253 elif argindex == 4:
254 # Diff wrt x2
255 return (1 - x2)**(b - 1)*x2**(a - 1)
256 else:
257 raise ArgumentIndexError(self, argindex)
259 def _eval_mpmath(self):
260 return betainc_mpmath_fix, self.args
262 def _eval_is_real(self):
263 if all(arg.is_real for arg in self.args):
264 return True
266 def _eval_conjugate(self):
267 return self.func(*map(conjugate, self.args))
269 def _eval_rewrite_as_Integral(self, a, b, x1, x2, **kwargs):
270 from sympy.integrals.integrals import Integral
271 t = Dummy('t')
272 return Integral(t**(a - 1)*(1 - t)**(b - 1), (t, x1, x2))
274 def _eval_rewrite_as_hyper(self, a, b, x1, x2, **kwargs):
275 from sympy.functions.special.hyper import hyper
276 return (x2**a * hyper((a, 1 - b), (a + 1,), x2) - x1**a * hyper((a, 1 - b), (a + 1,), x1)) / a
278###############################################################################
279#################### REGULARIZED INCOMPLETE BETA FUNCTION #####################
280###############################################################################
282class betainc_regularized(Function):
283 r"""
284 The Generalized Regularized Incomplete Beta function is given by
286 .. math::
287 \mathrm{I}_{(x_1, x_2)}(a, b) = \frac{\mathrm{B}_{(x_1, x_2)}(a, b)}{\mathrm{B}(a, b)}
289 The Regularized Incomplete Beta function is a special case
290 of the Generalized Regularized Incomplete Beta function :
292 .. math:: \mathrm{I}_z (a, b) = \mathrm{I}_{(0, z)}(a, b)
294 The Regularized Incomplete Beta function is the cumulative distribution
295 function of the beta distribution.
297 Examples
298 ========
300 >>> from sympy import betainc_regularized, symbols, conjugate
301 >>> a, b, x, x1, x2 = symbols('a b x x1 x2')
303 The Generalized Regularized Incomplete Beta
304 function is given by:
306 >>> betainc_regularized(a, b, x1, x2)
307 betainc_regularized(a, b, x1, x2)
309 The Regularized Incomplete Beta function
310 can be obtained as follows:
312 >>> betainc_regularized(a, b, 0, x)
313 betainc_regularized(a, b, 0, x)
315 The Regularized Incomplete Beta function
316 obeys the mirror symmetry:
318 >>> conjugate(betainc_regularized(a, b, x1, x2))
319 betainc_regularized(conjugate(a), conjugate(b), conjugate(x1), conjugate(x2))
321 We can numerically evaluate the Regularized Incomplete Beta function
322 to arbitrary precision for any complex numbers a, b, x1 and x2:
324 >>> from sympy import betainc_regularized, pi, E
325 >>> betainc_regularized(1, 2, 0, 0.25).evalf(10)
326 0.4375000000
327 >>> betainc_regularized(pi, E, 0, 1).evalf(5)
328 1.00000
330 The Generalized Regularized Incomplete Beta function can be
331 expressed in terms of the Generalized Hypergeometric function.
333 >>> from sympy import hyper
334 >>> betainc_regularized(a, b, x1, x2).rewrite(hyper)
335 (-x1**a*hyper((a, 1 - b), (a + 1,), x1) + x2**a*hyper((a, 1 - b), (a + 1,), x2))/(a*beta(a, b))
337 See Also
338 ========
340 beta: Beta function
341 hyper: Generalized Hypergeometric function
343 References
344 ==========
346 .. [1] https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function
347 .. [2] https://dlmf.nist.gov/8.17
348 .. [3] https://functions.wolfram.com/GammaBetaErf/Beta4/
349 .. [4] https://functions.wolfram.com/GammaBetaErf/BetaRegularized4/02/
351 """
352 nargs = 4
353 unbranched = True
355 def __new__(cls, a, b, x1, x2):
356 return Function.__new__(cls, a, b, x1, x2)
358 def _eval_mpmath(self):
359 return betainc_mpmath_fix, (*self.args, S(1))
361 def fdiff(self, argindex):
362 a, b, x1, x2 = self.args
363 if argindex == 3:
364 # Diff wrt x1
365 return -(1 - x1)**(b - 1)*x1**(a - 1) / beta(a, b)
366 elif argindex == 4:
367 # Diff wrt x2
368 return (1 - x2)**(b - 1)*x2**(a - 1) / beta(a, b)
369 else:
370 raise ArgumentIndexError(self, argindex)
372 def _eval_is_real(self):
373 if all(arg.is_real for arg in self.args):
374 return True
376 def _eval_conjugate(self):
377 return self.func(*map(conjugate, self.args))
379 def _eval_rewrite_as_Integral(self, a, b, x1, x2, **kwargs):
380 from sympy.integrals.integrals import Integral
381 t = Dummy('t')
382 integrand = t**(a - 1)*(1 - t)**(b - 1)
383 expr = Integral(integrand, (t, x1, x2))
384 return expr / Integral(integrand, (t, 0, 1))
386 def _eval_rewrite_as_hyper(self, a, b, x1, x2, **kwargs):
387 from sympy.functions.special.hyper import hyper
388 expr = (x2**a * hyper((a, 1 - b), (a + 1,), x2) - x1**a * hyper((a, 1 - b), (a + 1,), x1)) / a
389 return expr / beta(a, b)