Efficient integration of bell-shaped functions#
- quad_bellshaped(f, *points, **kwargs)#
where
ctxisipm,dec,mpm, orgmp.Computes a single, double or triple integral over a given, 1D interval, 2D rectangle, or 3D cuboid, providing optionalle a guaranteed error bound.
Introduction#
Some text explaining search for optimum and using 4 standard subintervals
Error function#
- ctx.real_quad_erf(a, x)#
where
ctxisipm,dec,mpm, orgmp.Returns the error function of \(\text{erf}(z)\). See also BoostMath [80], Wikipedia [1340] , MathWorld [916], NIST [843], Mpmath [571]. The error function defined by
\[\text{erf}(z) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-z^2} \mathrm{d} t,\]References for integration with erf(z) as argument:
https://www.wolframalpha.com/input/?i=integrate+exp%28x%5E2%29
http://fredrikj.net/blog/2016/03/taking-the-error-out-of-the-error-function/
Bounds on the complex error function:
\[\text{erf}{(x+i y)} = \text{erf}{x} \pm i e^{-x^2} \text{erfi}{y} \pm e^{-x^2} \text{erfi}{y}.\]
Lower non-normalised incomplete gamma function, \(\gamma(a,x)\) (Continued fractions)#
- ctx.real_quad_gamma_lower(a, x)#
where
ctxisipm,dec,mpm, orgmp.Returns the real lower non-normalised incomplete gamma function \(\gamma(a,x)\). See also Wikipedia [1322], MathWorld [921], NIST [482], BoostMath [87], Mpmath [576]. The function is defined as
\[\gamma(a,x)= \int_0^x t^{a-1} e^{-t}\mathrm{d} t\]for \(a \geq 0\) and \(x \geq 0\).
Real upper non-normalised incomplete gamma function, \(\Gamma(a,x)\) (Quadrature)#
- ctx.real_quad_gamma_upper(a, x)#
where
ctxisipm,dec,mpm, orgmp.Returns the real lower non-normalised incomplete gamma function \(\Gamma(a,x)\). See also Wikipedia [1322], MathWorld [921], NIST [482], BoostMath [87], Mpmath [576]. The function is defined as:
\[\Gamma(a,x) = \int_x^{\infty} t^{a-1} e^{-t}\mathrm{d} t\]for \(a \geq 0\) and \(x \geq 0\).
Normalised incomplete beta function, \(I_{x}(a,b)\)#
- ctx.real_quad_ibeta(a, b, x)#
where
ctxisipm,dec,mpm, orgmp.Returns the normalised incomplete beta function \(I_x(a,b)\) for \(a>0\), \(b>0\), and \(0 \leq x \leq 1\):
\[I_x(a,b) = \frac{B_x(a,b)}{B(a,b)}, \quad B_x(a,b) = \int_0^x t^{a-1} (1-t)^{b-1} \mathrm{d} t.\]See also Wikipedia [1321], MathWorld [927], NIST [481], BoostMath [85], BoostMath [93], Mpmath [575].
>>> from mpformula import mp4 >>> mp4.dps = 30; >>> matrix(2) matrix( [['0.0', '0.0'], ['0.0', '0.0']]) >>> matrix(2, 3)
Non-central chi-square cdf and sf (Chow)#
- ctx.chi_squared_nc_quad_cdf(n, x, lambda)#
where
ctxisipm,dec,mpm, orgmp.[Chou1985` gives the following representation for \(n \geq 2\) and \(\lambda \geq 0\):
\[F_{\chi^2}\left(n, x; \lambda\right) = \frac{2^{(1-n)/2}\sqrt{2\pi}}{ \Gamma((n-1)/2))} \int_{0}^{x} y^{(n-3)/2} \phi \left(\sqrt{y}\right) \left[\Phi \left(\sqrt{x-y}-\sqrt{\lambda}\right) - \Phi \left(-\sqrt{x-y}-\sqrt{\lambda}\right) \right]\mathrm{d} y\]where \(\phi(\cdot)\) denotes the pdf of the normal distribution (see section ref{sec:NormalDistribution_pdf}) and \(\Phi(\cdot)\) denotes the cdf of the normal distribution (see section ref{sec:NormalDistribution_CDF}).
Marcum \(Q\) function#
- ctx.marcumq1(a, b)#
where
ctxisfpm,mpm,ipm,dec,gmporapm.The Marcum Q-function is defined as
\[Q(a,b)=\int _{b}^{\infty} x \exp \left(-{\frac {x^{2}+a^{2}}{2}}\right)I_{0}(ax)\,\mathrm{d} x\]where \(b\geq 0\), \(a>0\) and \(I_{0}\) is the modified Bessel function of first kind of order 0. It is a special case of the generalized Marcum \(Q_m(a,b)\) function with \(m=1\). It is made available as a separate function because it can be calculated in a particularly simple way.
Special cases are
\[Q(a,0)=1, \quad Q(0,b)=e^{-b^2/2}, \quad Q(a,a) = \tfrac{1}{2} \left(1 + \frac{I_0(a^2)}{\exp(a^2)} \right)\]The function can also be evaluated as
\[\begin{split}Q(a,b) =\begin{cases} H(a,b), & a<b,\\ \tfrac{1}{2} + H(a,a), & a=b,\\ 1+H(a,b) & a>b, \end{cases}\end{split}\]where
\[H(a,b) = \frac{1}{\pi} \exp \left( -\frac{a^2+b^2}{2} \right) \int _{0}^{\pi} G(t) \exp(a b \cos(t)) \mathrm{d} t, \quad \text{and}\]\[\begin{split}G(t) =\begin{cases} \tfrac{1}{2}, & t=0, a=b\\ \frac{1-z \cos(t)}{1-2 z \cos(t)+z^2} & \text{otherwise, with } z = \frac{a}{b} \end{cases}\end{split}\]Available in Amath.
See also: https://en.wikipedia.org/wiki/Marcum_Q-function
See also: Short, 2012
See also Morales-Jimenez, 2013
See Marcum, 1960
>>> from mpfunlab import dec, mpm, ipm >>> mpm.dps = 40; a = '3.7'; b = '10.3' >>> dx = dec.marcumq(a, b); mx = mpm.marcumq(a, b); ix = ipm.marcumq(a, b) >>> mpm.show([dx, mx, ix]) dec: 3.465377030836928861042229075522437177483E-11 mpm: 3.465377030836928861042229075522437177483e-11 ipm: 3.465377030836928861042229075522437177486e-11 (9.641e-40%) >>> from mpfunlab import mpm, fpm, gmp, apm >>> mpm.dps = 40; a = '3.7'; b = '10.3' >>> fx = fpm.marcumq(a, b); gx = gmp.marcumq(a, b); ax = apm.marcumq(a, b) >>> mpm.show([gx, fx, ax]) gmp: 3.465377030836928861042229075522437177483E-11 fpm: 3.46537703083691E-11 apm: 3.465377030836928861042229075522437177486e-11 (9.641e-40%)
Owen’s \(T\) function#
- ctx.owent(h, a, boost=True)#
where
ctxisfpm,mpm,ipm,dec,gmporapm.Returns Owen’s \(T\) function \(T(h, a)\). See also Wikipedia [1460], MathWorld [1088], Owen [476], and Patefield and Tand [486]. The function is defined by
\[T(h,a) = \frac {1}{2\pi } \int _{0}^{a} f(x) \mathrm{d} x = \frac {a}{4\pi } \int _{-1}^{1} f(ax) \mathrm{d} x, \quad f(x) = {\frac {e^{-{\frac {1}{2}}h^{2}(1+x^{2})}}{1+x^{2}}}, \quad \left(-\infty <h,a<+\infty \right).\]It has the following properties:
\[T(h,0)=0\]\[T(0,a)={\frac {1}{2\pi }}\arctan(a)\]\[T(-h,a)=T(h,a)\]\[T(h,-a)=-T(h,a)\]\[T(h,a)+T(ah,{\frac {1}{a}})={\frac {1}{2}}\left(\Phi (h)+\Phi (ah)\right)-\Phi (h)\Phi (ah)\quad {\text{if}}\quad a\geq 0\]\[T(h,a)+T(ah,\frac {1}{a}) = \frac {1}{2} \left(\Phi (h)+\Phi (ah)\right)-\Phi (h)\Phi (ah)-{\frac {1}{2}} \quad \text{if}\quad a<0\]\[T(h,1) = \frac{1}{2} \Phi(h) \left( 1-\Phi(h) \right).\]>>> from mpfunlab import dec, mpm, ipm >>> mpm.dps = 40; h = 3.7; a = 10.3 >>> dx = dec.owent(h, a); mx = mpm.owent(h, a); ix = ipm.owent(h, a) >>> mpm.show([dx, mx, ix]) dec: 5.389986673869416846873471643522652586897E-5 mpm: 5.389986673869413074066777467941111347102e-5 ipm: 5.389986673869413074066777467941111347102e-5 (6.5e-40%) >>> from mpfunlab import mpm, fpm, gmp, apm >>> mpm.dps = 40; h = 3.7; a = 10.3 >>> fx = fpm.owent(h, a); gx = gmp.owent(h, a); ax = apm.owent(h, a) >>> mpm.show([gx, fx, ax]) gmp: 5.389986673869416846873471643522652586897E-05 fpm: 5.38998667386941E-05 apm: 5.389986673869413074066777467941111347102e-5 (6.5e-40%) >>> fx = fpm.owent(h, a); mpm.show([fx]) # boost fpm: 5.38998667386941E-05 >>> fx = fpm.owent(h, a, False); mpm.show([fx]) # amath fpm: 5.38998667386941E-05