SciPy
es paquete que incluye una colección de algoritmos matemáticos y funciones construidas sobre el paquete NumPy
. En definitiva, nos proporciona una serie de comandos de alto nivel para realizar acciones como:
Incluso otras más especializadas como:
Cada una de estas tareas pertenece a un subpaquete distinto de SciPy
. En la clase de hoy nos centraremos en las tres primeras, pero seguro que si quieres hincarle el diente a alguna de las otras no tienes ningún problema.
Subpackage | Description |
---|---|
constants | Physical and mathematical constants |
integrate | Integration and ordinary differential equation solvers |
interpolate | Interpolation and smoothing splines |
linalg | Linear algebra |
signal | Signal processing |
special | Special functions |
stats | Statistical distributions and functions |
Como siempre lo primero es lo primero, importemos lo paquetes que vamos a utilizar:
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
¿Cómo importamos SciPy?
import scipy as ...
Bueno, en realidad SciPy
nos ofrece muchas funciones y rara vez querremos usarlo todo en una misma sesión. Típicamente usaremos:
from scipy import ...
from scipy import special
np.info(special)
======================================== Special functions (:mod:`scipy.special`) ======================================== .. module:: scipy.special Nearly all of the functions below are universal functions and follow broadcasting and automatic array-looping rules. Exceptions are noted. Error handling ============== Errors are handled by returning nans, or other appropriate values. Some of the special function routines will emit warnings when an error occurs. By default this is disabled. To enable such messages use ``errprint(1)``, and to disable such messages use ``errprint(0)``. Example: >>> print scipy.special.bdtr(-1,10,0.3) >>> scipy.special.errprint(1) >>> print scipy.special.bdtr(-1,10,0.3) .. autosummary:: :toctree: generated/ errprint Available functions =================== Airy functions -------------- .. autosummary:: :toctree: generated/ airy -- Airy functions and their derivatives. airye -- Exponentially scaled Airy functions ai_zeros -- [+]Zeros of Airy functions Ai(x) and Ai'(x) bi_zeros -- [+]Zeros of Airy functions Bi(x) and Bi'(x) Elliptic Functions and Integrals -------------------------------- .. autosummary:: :toctree: generated/ ellipj -- Jacobian elliptic functions ellipk -- Complete elliptic integral of the first kind. ellipkm1 -- ellipkm1(x) == ellipk(1 - x) ellipkinc -- Incomplete elliptic integral of the first kind. ellipe -- Complete elliptic integral of the second kind. ellipeinc -- Incomplete elliptic integral of the second kind. Bessel Functions ---------------- .. autosummary:: :toctree: generated/ jn -- Bessel function of integer order and real argument. jv -- Bessel function of real-valued order and complex argument. jve -- Exponentially scaled Bessel function. yn -- Bessel function of second kind (integer order). yv -- Bessel function of the second kind (real-valued order). yve -- Exponentially scaled Bessel function of the second kind. kn -- Modified Bessel function of the second kind (integer order). kv -- Modified Bessel function of the second kind (real order). kve -- Exponentially scaled modified Bessel function of the second kind. iv -- Modified Bessel function. ive -- Exponentially scaled modified Bessel function. hankel1 -- Hankel function of the first kind. hankel1e -- Exponentially scaled Hankel function of the first kind. hankel2 -- Hankel function of the second kind. hankel2e -- Exponentially scaled Hankel function of the second kind. The following is not an universal function: .. autosummary:: :toctree: generated/ lmbda -- [+]Sequence of lambda functions with arbitrary order v. Zeros of Bessel Functions ^^^^^^^^^^^^^^^^^^^^^^^^^ These are not universal functions: .. autosummary:: :toctree: generated/ jnjnp_zeros -- [+]Zeros of integer-order Bessel functions and derivatives sorted in order. jnyn_zeros -- [+]Zeros of integer-order Bessel functions and derivatives as separate arrays. jn_zeros -- [+]Zeros of Jn(x) jnp_zeros -- [+]Zeros of Jn'(x) yn_zeros -- [+]Zeros of Yn(x) ynp_zeros -- [+]Zeros of Yn'(x) y0_zeros -- [+]Complex zeros: Y0(z0)=0 and values of Y0'(z0) y1_zeros -- [+]Complex zeros: Y1(z1)=0 and values of Y1'(z1) y1p_zeros -- [+]Complex zeros of Y1'(z1')=0 and values of Y1(z1') Faster versions of common Bessel Functions ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. autosummary:: :toctree: generated/ j0 -- Bessel function of order 0. j1 -- Bessel function of order 1. y0 -- Bessel function of second kind of order 0. y1 -- Bessel function of second kind of order 1. i0 -- Modified Bessel function of order 0. i0e -- Exponentially scaled modified Bessel function of order 0. i1 -- Modified Bessel function of order 1. i1e -- Exponentially scaled modified Bessel function of order 1. k0 -- Modified Bessel function of the second kind of order 0. k0e -- Exponentially scaled modified Bessel function of the second kind of order 0. k1 -- Modified Bessel function of the second kind of order 1. k1e -- Exponentially scaled modified Bessel function of the second kind of order 1. Integrals of Bessel Functions ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. autosummary:: :toctree: generated/ itj0y0 -- Basic integrals of j0 and y0 from 0 to x. it2j0y0 -- Integrals of (1-j0(t))/t from 0 to x and y0(t)/t from x to inf. iti0k0 -- Basic integrals of i0 and k0 from 0 to x. it2i0k0 -- Integrals of (i0(t)-1)/t from 0 to x and k0(t)/t from x to inf. besselpoly -- Integral of a Bessel function: Jv(2* a* x) * x[+]lambda from x=0 to 1. Derivatives of Bessel Functions ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. autosummary:: :toctree: generated/ jvp -- Nth derivative of Jv(v,z) yvp -- Nth derivative of Yv(v,z) kvp -- Nth derivative of Kv(v,z) ivp -- Nth derivative of Iv(v,z) h1vp -- Nth derivative of H1v(v,z) h2vp -- Nth derivative of H2v(v,z) Spherical Bessel Functions ^^^^^^^^^^^^^^^^^^^^^^^^^^ These are not universal functions: .. autosummary:: :toctree: generated/ sph_jn -- [+]Sequence of spherical Bessel functions, jn(z) sph_yn -- [+]Sequence of spherical Bessel functions, yn(z) sph_jnyn -- [+]Sequence of spherical Bessel functions, jn(z) and yn(z) sph_in -- [+]Sequence of spherical Bessel functions, in(z) sph_kn -- [+]Sequence of spherical Bessel functions, kn(z) sph_inkn -- [+]Sequence of spherical Bessel functions, in(z) and kn(z) Riccati-Bessel Functions ^^^^^^^^^^^^^^^^^^^^^^^^ These are not universal functions: .. autosummary:: :toctree: generated/ riccati_jn -- [+]Sequence of Ricatti-Bessel functions of first kind. riccati_yn -- [+]Sequence of Ricatti-Bessel functions of second kind. Struve Functions ---------------- .. autosummary:: :toctree: generated/ struve -- Struve function --- Hv(x) modstruve -- Modified Struve function --- Lv(x) itstruve0 -- Integral of H0(t) from 0 to x it2struve0 -- Integral of H0(t)/t from x to Inf. itmodstruve0 -- Integral of L0(t) from 0 to x. Raw Statistical Functions ------------------------- .. seealso:: :mod:`scipy.stats`: Friendly versions of these functions. .. autosummary:: :toctree: generated/ bdtr -- Sum of terms 0 through k of the binomial pdf. bdtrc -- Sum of terms k+1 through n of the binomial pdf. bdtri -- Inverse of bdtr btdtr -- Integral from 0 to x of beta pdf. btdtri -- Quantiles of beta distribution fdtr -- Integral from 0 to x of F pdf. fdtrc -- Integral from x to infinity under F pdf. fdtri -- Inverse of fdtrc gdtr -- Integral from 0 to x of gamma pdf. gdtrc -- Integral from x to infinity under gamma pdf. gdtria -- Inverse with respect to `a` of gdtr. gdtrib -- Inverse with respect to `b` of gdtr. gdtrix -- Inverse with respect to `x` of gdtr. nbdtr -- Sum of terms 0 through k of the negative binomial pdf. nbdtrc -- Sum of terms k+1 to infinity under negative binomial pdf. nbdtri -- Inverse of nbdtr pdtr -- Sum of terms 0 through k of the Poisson pdf. pdtrc -- Sum of terms k+1 to infinity of the Poisson pdf. pdtri -- Inverse of pdtr stdtr -- Integral from -infinity to t of the Student-t pdf. stdtridf -- stdtrit -- chdtr -- Integral from 0 to x of the Chi-square pdf. chdtrc -- Integral from x to infnity of Chi-square pdf. chdtri -- Inverse of chdtrc. ndtr -- Integral from -infinity to x of standard normal pdf ndtri -- Inverse of ndtr (quantiles) smirnov -- Kolmogorov-Smirnov complementary CDF for one-sided test statistic (Dn+ or Dn-) smirnovi -- Inverse of smirnov. kolmogorov -- The complementary CDF of the (scaled) two-sided test statistic (Kn*) valid for large n. kolmogi -- Inverse of kolmogorov tklmbda -- Tukey-Lambda CDF logit -- expit -- boxcox -- Compute the Box-Cox transformation. boxcox1p -- Compute the Box-Cox transformation. Gamma and Related Functions --------------------------- .. autosummary:: :toctree: generated/ gamma -- Gamma function. gammaln -- Log of the absolute value of the gamma function. gammasgn -- Sign of the gamma function. gammainc -- Incomplete gamma integral. gammaincinv -- Inverse of gammainc. gammaincc -- Complemented incomplete gamma integral. gammainccinv -- Inverse of gammaincc. beta -- Beta function. betaln -- Log of the absolute value of the beta function. betainc -- Incomplete beta integral. betaincinv -- Inverse of betainc. psi -- Logarithmic derivative of the gamma function. rgamma -- One divided by the gamma function. polygamma -- Nth derivative of psi function. multigammaln Error Function and Fresnel Integrals ------------------------------------ .. autosummary:: :toctree: generated/ erf -- Error function. erfc -- Complemented error function (1- erf(x)) erfcx -- Scaled complemented error function exp(x**2)*erfc(x) erfi -- Imaginary error function, -i erf(i x) erfinv -- Inverse of error function erfcinv -- Inverse of erfc wofz -- Fadeeva function. dawsn -- Dawson's integral. fresnel -- Fresnel sine and cosine integrals. fresnel_zeros -- Complex zeros of both Fresnel integrals modfresnelp -- Modified Fresnel integrals F_+(x) and K_+(x) modfresnelm -- Modified Fresnel integrals F_-(x) and K_-(x) These are not universal functions: .. autosummary:: :toctree: generated/ erf_zeros -- [+]Complex zeros of erf(z) fresnelc_zeros -- [+]Complex zeros of Fresnel cosine integrals fresnels_zeros -- [+]Complex zeros of Fresnel sine integrals Legendre Functions ------------------ .. autosummary:: :toctree: generated/ lpmv -- Associated Legendre Function of arbitrary non-negative degree v. sph_harm -- Spherical Harmonics (complex-valued) Y^m_n(theta,phi) These are not universal functions: .. autosummary:: :toctree: generated/ clpmn -- [+]Associated Legendre Function of the first kind for complex arguments. lpn -- [+]Legendre Functions (polynomials) of the first kind lqn -- [+]Legendre Functions of the second kind. lpmn -- [+]Associated Legendre Function of the first kind for real arguments. lqmn -- [+]Associated Legendre Function of the second kind. Orthogonal polynomials ---------------------- The following functions evaluate values of orthogonal polynomials: .. autosummary:: :toctree: generated/ eval_legendre eval_chebyt eval_chebyu eval_chebyc eval_chebys eval_jacobi eval_laguerre eval_genlaguerre eval_hermite eval_hermitenorm eval_gegenbauer eval_sh_legendre eval_sh_chebyt eval_sh_chebyu eval_sh_jacobi The functions below, in turn, return :ref:`orthopoly1d` objects, which functions similarly as :ref:`numpy.poly1d`. The :ref:`orthopoly1d` class also has an attribute ``weights`` which returns the roots, weights, and total weights for the appropriate form of Gaussian quadrature. These are returned in an ``n x 3`` array with roots in the first column, weights in the second column, and total weights in the final column. .. autosummary:: :toctree: generated/ legendre -- [+]Legendre polynomial P_n(x) (lpn -- for function). chebyt -- [+]Chebyshev polynomial T_n(x) chebyu -- [+]Chebyshev polynomial U_n(x) chebyc -- [+]Chebyshev polynomial C_n(x) chebys -- [+]Chebyshev polynomial S_n(x) jacobi -- [+]Jacobi polynomial P^(alpha,beta)_n(x) laguerre -- [+]Laguerre polynomial, L_n(x) genlaguerre -- [+]Generalized (Associated) Laguerre polynomial, L^alpha_n(x) hermite -- [+]Hermite polynomial H_n(x) hermitenorm -- [+]Normalized Hermite polynomial, He_n(x) gegenbauer -- [+]Gegenbauer (Ultraspherical) polynomials, C^(alpha)_n(x) sh_legendre -- [+]shifted Legendre polynomial, P*_n(x) sh_chebyt -- [+]shifted Chebyshev polynomial, T*_n(x) sh_chebyu -- [+]shifted Chebyshev polynomial, U*_n(x) sh_jacobi -- [+]shifted Jacobi polynomial, J*_n(x) = G^(p,q)_n(x) .. warning:: Large-order polynomials obtained from these functions are numerically unstable. ``orthopoly1d`` objects are converted to ``poly1d``, when doing arithmetic. ``numpy.poly1d`` works in power basis and cannot represent high-order polynomials accurately, which can cause significant inaccuracy. Hypergeometric Functions ------------------------ .. autosummary:: :toctree: generated/ hyp2f1 -- Gauss hypergeometric function (2F1) hyp1f1 -- Confluent hypergeometric function (1F1) hyperu -- Confluent hypergeometric function (U) hyp0f1 -- Confluent hypergeometric limit function (0F1) hyp2f0 -- Hypergeometric function (2F0) hyp1f2 -- Hypergeometric function (1F2) hyp3f0 -- Hypergeometric function (3F0) Parabolic Cylinder Functions ---------------------------- .. autosummary:: :toctree: generated/ pbdv -- Parabolic cylinder function Dv(x) and derivative. pbvv -- Parabolic cylinder function Vv(x) and derivative. pbwa -- Parabolic cylinder function W(a,x) and derivative. These are not universal functions: .. autosummary:: :toctree: generated/ pbdv_seq -- [+]Sequence of parabolic cylinder functions Dv(x) pbvv_seq -- [+]Sequence of parabolic cylinder functions Vv(x) pbdn_seq -- [+]Sequence of parabolic cylinder functions Dn(z), complex z Mathieu and Related Functions ----------------------------- .. autosummary:: :toctree: generated/ mathieu_a -- Characteristic values for even solution (ce_m) mathieu_b -- Characteristic values for odd solution (se_m) These are not universal functions: .. autosummary:: :toctree: generated/ mathieu_even_coef -- [+]sequence of expansion coefficients for even solution mathieu_odd_coef -- [+]sequence of expansion coefficients for odd solution The following return both function and first derivative: .. autosummary:: :toctree: generated/ mathieu_cem -- Even Mathieu function mathieu_sem -- Odd Mathieu function mathieu_modcem1 -- Even modified Mathieu function of the first kind mathieu_modcem2 -- Even modified Mathieu function of the second kind mathieu_modsem1 -- Odd modified Mathieu function of the first kind mathieu_modsem2 -- Odd modified Mathieu function of the second kind Spheroidal Wave Functions ------------------------- .. autosummary:: :toctree: generated/ pro_ang1 -- Prolate spheroidal angular function of the first kind pro_rad1 -- Prolate spheroidal radial function of the first kind pro_rad2 -- Prolate spheroidal radial function of the second kind obl_ang1 -- Oblate spheroidal angular function of the first kind obl_rad1 -- Oblate spheroidal radial function of the first kind obl_rad2 -- Oblate spheroidal radial function of the second kind pro_cv -- Compute characteristic value for prolate functions obl_cv -- Compute characteristic value for oblate functions pro_cv_seq -- Compute sequence of prolate characteristic values obl_cv_seq -- Compute sequence of oblate characteristic values The following functions require pre-computed characteristic value: .. autosummary:: :toctree: generated/ pro_ang1_cv -- Prolate spheroidal angular function of the first kind pro_rad1_cv -- Prolate spheroidal radial function of the first kind pro_rad2_cv -- Prolate spheroidal radial function of the second kind obl_ang1_cv -- Oblate spheroidal angular function of the first kind obl_rad1_cv -- Oblate spheroidal radial function of the first kind obl_rad2_cv -- Oblate spheroidal radial function of the second kind Kelvin Functions ---------------- .. autosummary:: :toctree: generated/ kelvin -- All Kelvin functions (order 0) and derivatives. kelvin_zeros -- [+]Zeros of All Kelvin functions (order 0) and derivatives ber -- Kelvin function ber x bei -- Kelvin function bei x berp -- Derivative of Kelvin function ber x beip -- Derivative of Kelvin function bei x ker -- Kelvin function ker x kei -- Kelvin function kei x kerp -- Derivative of Kelvin function ker x keip -- Derivative of Kelvin function kei x These are not universal functions: .. autosummary:: :toctree: generated/ ber_zeros -- [+]Zeros of Kelvin function bei x bei_zeros -- [+]Zeros of Kelvin function ber x berp_zeros -- [+]Zeros of derivative of Kelvin function ber x beip_zeros -- [+]Zeros of derivative of Kelvin function bei x ker_zeros -- [+]Zeros of Kelvin function kei x kei_zeros -- [+]Zeros of Kelvin function ker x kerp_zeros -- [+]Zeros of derivative of Kelvin function ker x keip_zeros -- [+]Zeros of derivative of Kelvin function kei x Combinatorics ------------- .. autosummary:: :toctree: generated/ comb -- [+]Combinations of N things taken k at a time, "N choose k" perm -- [+]Permutations of N things taken k at a time, "k-permutations of N" Other Special Functions ----------------------- .. autosummary:: :toctree: generated/ binom -- Binomial coefficient. expn -- Exponential integral. exp1 -- Exponential integral of order 1 (for complex argument) expi -- Another exponential integral -- Ei(x) factorial -- The factorial function, n! = special.gamma(n+1) factorial2 -- Double factorial, (n!)! factorialk -- [+](...((n!)!)!...)! where there are k '!' shichi -- Hyperbolic sine and cosine integrals. sici -- Integral of the sinc and "cosinc" functions. spence -- Dilogarithm integral. lambertw -- Lambert W function zeta -- Riemann zeta function of two arguments. zetac -- Standard Riemann zeta function minus 1. Convenience Functions --------------------- .. autosummary:: :toctree: generated/ cbrt -- Cube root. exp10 -- 10 raised to the x power. exp2 -- 2 raised to the x power. radian -- radian angle given degrees, minutes, and seconds. cosdg -- cosine of the angle given in degrees. sindg -- sine of the angle given in degrees. tandg -- tangent of the angle given in degrees. cotdg -- cotangent of the angle given in degrees. log1p -- log(1+x) expm1 -- exp(x)-1 cosm1 -- cos(x)-1 round -- round the argument to the nearest integer. If argument ends in 0.5 exactly, pick the nearest even integer. xlogy -- x*log(y) xlog1py -- x*log1p(y) .. [+] in the description indicates a function which is not a universal .. function and does not follow broadcasting and automatic .. array-looping rules.
x = np.linspace(0,20,100)
for ii in range(5):
plt.plot(x, special.jn(ii, x))
plt.grid(True)
Este subpaquete de SciPy
proporciona algunas técnicas de integración tanto de funciones como de ecuaciones diferenciales. En primer lugar importémoslo y ejecutemos la ayuda para ver cuáles son estas funciones:
from scipy import integrate
from IPython.display import HTML
HTML('<iframe src="http://docs.scipy.org/doc/scipy/reference/integrate.html#module-scipy.integrate" width="800" height="600"></iframe>')
help(integrate)
Help on package scipy.integrate in scipy: NAME scipy.integrate DESCRIPTION ============================================= Integration and ODEs (:mod:`scipy.integrate`) ============================================= .. currentmodule:: scipy.integrate Integrating functions, given function object ============================================ .. autosummary:: :toctree: generated/ quad -- General purpose integration dblquad -- General purpose double integration tplquad -- General purpose triple integration nquad -- General purpose n-dimensional integration fixed_quad -- Integrate func(x) using Gaussian quadrature of order n quadrature -- Integrate with given tolerance using Gaussian quadrature romberg -- Integrate func using Romberg integration Integrating functions, given fixed samples ========================================== .. autosummary:: :toctree: generated/ cumtrapz -- Use trapezoidal rule to cumulatively compute integral. simps -- Use Simpson's rule to compute integral from samples. romb -- Use Romberg Integration to compute integral from -- (2**k + 1) evenly-spaced samples. .. seealso:: :mod:`scipy.special` for orthogonal polynomials (special) for Gaussian quadrature roots and weights for other weighting factors and regions. Integrators of ODE systems ========================== .. autosummary:: :toctree: generated/ odeint -- General integration of ordinary differential equations. ode -- Integrate ODE using VODE and ZVODE routines. complex_ode -- Convert a complex-valued ODE to real-valued and integrate. PACKAGE CONTENTS _dop _ode _odepack _quadpack lsoda odepack quadpack quadrature setup vode CLASSES builtins.UserWarning(builtins.Warning) scipy.integrate.quadpack.IntegrationWarning builtins.object scipy.integrate._ode.ode scipy.integrate._ode.complex_ode class IntegrationWarning(builtins.UserWarning) | Method resolution order: | IntegrationWarning | builtins.UserWarning | builtins.Warning | builtins.Exception | builtins.BaseException | builtins.object | | Data descriptors defined here: | | __weakref__ | list of weak references to the object (if defined) | | ---------------------------------------------------------------------- | Methods inherited from builtins.UserWarning: | | __init__(self, /, *args, **kwargs) | Initialize self. See help(type(self)) for accurate signature. | | __new__(*args, **kwargs) from builtins.type | Create and return a new object. See help(type) for accurate signature. | | ---------------------------------------------------------------------- | Methods inherited from builtins.BaseException: | | __delattr__(self, name, /) | Implement delattr(self, name). | | __getattribute__(self, name, /) | Return getattr(self, name). | | __reduce__(...) | | __repr__(self, /) | Return repr(self). | | __setattr__(self, name, value, /) | Implement setattr(self, name, value). | | __setstate__(...) | | __str__(self, /) | Return str(self). | | with_traceback(...) | Exception.with_traceback(tb) -- | set self.__traceback__ to tb and return self. | | ---------------------------------------------------------------------- | Data descriptors inherited from builtins.BaseException: | | __cause__ | exception cause | | __context__ | exception context | | __dict__ | | __suppress_context__ | | __traceback__ | | args class complex_ode(ode) | A wrapper of ode for complex systems. | | This functions similarly as `ode`, but re-maps a complex-valued | equation system to a real-valued one before using the integrators. | | Parameters | ---------- | f : callable ``f(t, y, *f_args)`` | Rhs of the equation. t is a scalar, ``y.shape == (n,)``. | ``f_args`` is set by calling ``set_f_params(*args)``. | jac : callable ``jac(t, y, *jac_args)`` | Jacobian of the rhs, ``jac[i,j] = d f[i] / d y[j]``. | ``jac_args`` is set by calling ``set_f_params(*args)``. | | Attributes | ---------- | t : float | Current time. | y : ndarray | Current variable values. | | Examples | -------- | For usage examples, see `ode`. | | Method resolution order: | complex_ode | ode | builtins.object | | Methods defined here: | | __init__(self, f, jac=None) | | integrate(self, t, step=0, relax=0) | Find y=y(t), set y as an initial condition, and return y. | | set_initial_value(self, y, t=0.0) | Set initial conditions y(t) = y. | | set_integrator(self, name, **integrator_params) | Set integrator by name. | | Parameters | ---------- | name : str | Name of the integrator | integrator_params : | Additional parameters for the integrator. | | set_solout(self, solout) | Set callable to be called at every successful integration step. | | Parameters | ---------- | solout : callable | ``solout(t, y)`` is called at each internal integrator step, | t is a scalar providing the current independent position | y is the current soloution ``y.shape == (n,)`` | solout should return -1 to stop integration | otherwise it should return None or 0 | | ---------------------------------------------------------------------- | Data descriptors defined here: | | y | | ---------------------------------------------------------------------- | Methods inherited from ode: | | set_f_params(self, *args) | Set extra parameters for user-supplied function f. | | set_jac_params(self, *args) | Set extra parameters for user-supplied function jac. | | successful(self) | Check if integration was successful. | | ---------------------------------------------------------------------- | Data descriptors inherited from ode: | | __dict__ | dictionary for instance variables (if defined) | | __weakref__ | list of weak references to the object (if defined) class ode(builtins.object) | A generic interface class to numeric integrators. | | Solve an equation system :math:`y'(t) = f(t,y)` with (optional) ``jac = df/dy``. | | Parameters | ---------- | f : callable ``f(t, y, *f_args)`` | Rhs of the equation. t is a scalar, ``y.shape == (n,)``. | ``f_args`` is set by calling ``set_f_params(*args)``. | `f` should return a scalar, array or list (not a tuple). | jac : callable ``jac(t, y, *jac_args)`` | Jacobian of the rhs, ``jac[i,j] = d f[i] / d y[j]``. | ``jac_args`` is set by calling ``set_f_params(*args)``. | | Attributes | ---------- | t : float | Current time. | y : ndarray | Current variable values. | | See also | -------- | odeint : an integrator with a simpler interface based on lsoda from ODEPACK | quad : for finding the area under a curve | | Notes | ----- | Available integrators are listed below. They can be selected using | the `set_integrator` method. | | "vode" | | Real-valued Variable-coefficient Ordinary Differential Equation | solver, with fixed-leading-coefficient implementation. It provides | implicit Adams method (for non-stiff problems) and a method based on | backward differentiation formulas (BDF) (for stiff problems). | | Source: http://www.netlib.org/ode/vode.f | | .. warning:: | | This integrator is not re-entrant. You cannot have two `ode` | instances using the "vode" integrator at the same time. | | This integrator accepts the following parameters in `set_integrator` | method of the `ode` class: | | - atol : float or sequence | absolute tolerance for solution | - rtol : float or sequence | relative tolerance for solution | - lband : None or int | - rband : None or int | Jacobian band width, jac[i,j] != 0 for i-lband <= j <= i+rband. | Setting these requires your jac routine to return the jacobian | in packed format, jac_packed[i-j+lband, j] = jac[i,j]. | - method: 'adams' or 'bdf' | Which solver to use, Adams (non-stiff) or BDF (stiff) | - with_jacobian : bool | Whether to use the jacobian | - nsteps : int | Maximum number of (internally defined) steps allowed during one | call to the solver. | - first_step : float | - min_step : float | - max_step : float | Limits for the step sizes used by the integrator. | - order : int | Maximum order used by the integrator, | order <= 12 for Adams, <= 5 for BDF. | | "zvode" | | Complex-valued Variable-coefficient Ordinary Differential Equation | solver, with fixed-leading-coefficient implementation. It provides | implicit Adams method (for non-stiff problems) and a method based on | backward differentiation formulas (BDF) (for stiff problems). | | Source: http://www.netlib.org/ode/zvode.f | | .. warning:: | | This integrator is not re-entrant. You cannot have two `ode` | instances using the "zvode" integrator at the same time. | | This integrator accepts the same parameters in `set_integrator` | as the "vode" solver. | | .. note:: | | When using ZVODE for a stiff system, it should only be used for | the case in which the function f is analytic, that is, when each f(i) | is an analytic function of each y(j). Analyticity means that the | partial derivative df(i)/dy(j) is a unique complex number, and this | fact is critical in the way ZVODE solves the dense or banded linear | systems that arise in the stiff case. For a complex stiff ODE system | in which f is not analytic, ZVODE is likely to have convergence | failures, and for this problem one should instead use DVODE on the | equivalent real system (in the real and imaginary parts of y). | | "lsoda" | | Real-valued Variable-coefficient Ordinary Differential Equation | solver, with fixed-leading-coefficient implementation. It provides | automatic method switching between implicit Adams method (for non-stiff | problems) and a method based on backward differentiation formulas (BDF) | (for stiff problems). | | Source: http://www.netlib.org/odepack | | .. warning:: | | This integrator is not re-entrant. You cannot have two `ode` | instances using the "lsoda" integrator at the same time. | | This integrator accepts the following parameters in `set_integrator` | method of the `ode` class: | | - atol : float or sequence | absolute tolerance for solution | - rtol : float or sequence | relative tolerance for solution | - lband : None or int | - rband : None or int | Jacobian band width, jac[i,j] != 0 for i-lband <= j <= i+rband. | Setting these requires your jac routine to return the jacobian | in packed format, jac_packed[i-j+lband, j] = jac[i,j]. | - with_jacobian : bool | Whether to use the jacobian | - nsteps : int | Maximum number of (internally defined) steps allowed during one | call to the solver. | - first_step : float | - min_step : float | - max_step : float | Limits for the step sizes used by the integrator. | - max_order_ns : int | Maximum order used in the nonstiff case (default 12). | - max_order_s : int | Maximum order used in the stiff case (default 5). | - max_hnil : int | Maximum number of messages reporting too small step size (t + h = t) | (default 0) | - ixpr : int | Whether to generate extra printing at method switches (default False). | | "dopri5" | | This is an explicit runge-kutta method of order (4)5 due to Dormand & | Prince (with stepsize control and dense output). | | Authors: | | E. Hairer and G. Wanner | Universite de Geneve, Dept. de Mathematiques | CH-1211 Geneve 24, Switzerland | e-mail: ernst.hairer@math.unige.ch, gerhard.wanner@math.unige.ch | | This code is described in [HNW93]_. | | This integrator accepts the following parameters in set_integrator() | method of the ode class: | | - atol : float or sequence | absolute tolerance for solution | - rtol : float or sequence | relative tolerance for solution | - nsteps : int | Maximum number of (internally defined) steps allowed during one | call to the solver. | - first_step : float | - max_step : float | - safety : float | Safety factor on new step selection (default 0.9) | - ifactor : float | - dfactor : float | Maximum factor to increase/decrease step size by in one step | - beta : float | Beta parameter for stabilised step size control. | - verbosity : int | Switch for printing messages (< 0 for no messages). | | "dop853" | | This is an explicit runge-kutta method of order 8(5,3) due to Dormand | & Prince (with stepsize control and dense output). | | Options and references the same as "dopri5". | | Examples | -------- | | A problem to integrate and the corresponding jacobian: | | >>> from scipy.integrate import ode | >>> | >>> y0, t0 = [1.0j, 2.0], 0 | >>> | >>> def f(t, y, arg1): | >>> return [1j*arg1*y[0] + y[1], -arg1*y[1]**2] | >>> def jac(t, y, arg1): | >>> return [[1j*arg1, 1], [0, -arg1*2*y[1]]] | | The integration: | | >>> r = ode(f, jac).set_integrator('zvode', method='bdf', with_jacobian=True) | >>> r.set_initial_value(y0, t0).set_f_params(2.0).set_jac_params(2.0) | >>> t1 = 10 | >>> dt = 1 | >>> while r.successful() and r.t < t1: | >>> r.integrate(r.t+dt) | >>> print("%g %g" % (r.t, r.y)) | | References | ---------- | .. [HNW93] E. Hairer, S.P. Norsett and G. Wanner, Solving Ordinary | Differential Equations i. Nonstiff Problems. 2nd edition. | Springer Series in Computational Mathematics, | Springer-Verlag (1993) | | Methods defined here: | | __init__(self, f, jac=None) | | integrate(self, t, step=0, relax=0) | Find y=y(t), set y as an initial condition, and return y. | | set_f_params(self, *args) | Set extra parameters for user-supplied function f. | | set_initial_value(self, y, t=0.0) | Set initial conditions y(t) = y. | | set_integrator(self, name, **integrator_params) | Set integrator by name. | | Parameters | ---------- | name : str | Name of the integrator. | integrator_params : | Additional parameters for the integrator. | | set_jac_params(self, *args) | Set extra parameters for user-supplied function jac. | | set_solout(self, solout) | Set callable to be called at every successful integration step. | | Parameters | ---------- | solout : callable | ``solout(t, y)`` is called at each internal integrator step, | t is a scalar providing the current independent position | y is the current soloution ``y.shape == (n,)`` | solout should return -1 to stop integration | otherwise it should return None or 0 | | successful(self) | Check if integration was successful. | | ---------------------------------------------------------------------- | Data descriptors defined here: | | __dict__ | dictionary for instance variables (if defined) | | __weakref__ | list of weak references to the object (if defined) | | y FUNCTIONS cumtrapz(y, x=None, dx=1.0, axis=-1, initial=None) Cumulatively integrate y(x) using the composite trapezoidal rule. Parameters ---------- y : array_like Values to integrate. x : array_like, optional The coordinate to integrate along. If None (default), use spacing `dx` between consecutive elements in `y`. dx : int, optional Spacing between elements of `y`. Only used if `x` is None. axis : int, optional Specifies the axis to cumulate. Default is -1 (last axis). initial : scalar, optional If given, uses this value as the first value in the returned result. Typically this value should be 0. Default is None, which means no value at ``x[0]`` is returned and `res` has one element less than `y` along the axis of integration. Returns ------- res : ndarray The result of cumulative integration of `y` along `axis`. If `initial` is None, the shape is such that the axis of integration has one less value than `y`. If `initial` is given, the shape is equal to that of `y`. See Also -------- numpy.cumsum, numpy.cumprod quad: adaptive quadrature using QUADPACK romberg: adaptive Romberg quadrature quadrature: adaptive Gaussian quadrature fixed_quad: fixed-order Gaussian quadrature dblquad: double integrals tplquad: triple integrals romb: integrators for sampled data ode: ODE integrators odeint: ODE integrators Examples -------- >>> from scipy import integrate >>> import matplotlib.pyplot as plt >>> x = np.linspace(-2, 2, num=20) >>> y = x >>> y_int = integrate.cumtrapz(y, x, initial=0) >>> plt.plot(x, y_int, 'ro', x, y[0] + 0.5 * x**2, 'b-') >>> plt.show() dblquad(func, a, b, gfun, hfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08) Compute a double integral. Return the double (definite) integral of ``func(y, x)`` from ``x = a..b`` and ``y = gfun(x)..hfun(x)``. Parameters ---------- func : callable A Python function or method of at least two variables: y must be the first argument and x the second argument. (a,b) : tuple The limits of integration in x: `a` < `b` gfun : callable The lower boundary curve in y which is a function taking a single floating point argument (x) and returning a floating point result: a lambda function can be useful here. hfun : callable The upper boundary curve in y (same requirements as `gfun`). args : sequence, optional Extra arguments to pass to `func`. epsabs : float, optional Absolute tolerance passed directly to the inner 1-D quadrature integration. Default is 1.49e-8. epsrel : float Relative tolerance of the inner 1-D integrals. Default is 1.49e-8. Returns ------- y : float The resultant integral. abserr : float An estimate of the error. See also -------- quad : single integral tplquad : triple integral nquad : N-dimensional integrals fixed_quad : fixed-order Gaussian quadrature quadrature : adaptive Gaussian quadrature odeint : ODE integrator ode : ODE integrator simps : integrator for sampled data romb : integrator for sampled data scipy.special : for coefficients and roots of orthogonal polynomials fixed_quad(func, a, b, args=(), n=5) Compute a definite integral using fixed-order Gaussian quadrature. Integrate `func` from `a` to `b` using Gaussian quadrature of order `n`. Parameters ---------- func : callable A Python function or method to integrate (must accept vector inputs). a : float Lower limit of integration. b : float Upper limit of integration. args : tuple, optional Extra arguments to pass to function, if any. n : int, optional Order of quadrature integration. Default is 5. Returns ------- val : float Gaussian quadrature approximation to the integral See Also -------- quad : adaptive quadrature using QUADPACK dblquad : double integrals tplquad : triple integrals romberg : adaptive Romberg quadrature quadrature : adaptive Gaussian quadrature romb : integrators for sampled data simps : integrators for sampled data cumtrapz : cumulative integration for sampled data ode : ODE integrator odeint : ODE integrator newton_cotes(rn, equal=0) Return weights and error coefficient for Newton-Cotes integration. Suppose we have (N+1) samples of f at the positions x_0, x_1, ..., x_N. Then an N-point Newton-Cotes formula for the integral between x_0 and x_N is: :math:`\int_{x_0}^{x_N} f(x)dx = \Delta x \sum_{i=0}^{N} a_i f(x_i) + B_N (\Delta x)^{N+2} f^{N+1} (\xi)` where :math:`\xi \in [x_0,x_N]` and :math:`\Delta x = \frac{x_N-x_0}{N}` is the averages samples spacing. If the samples are equally-spaced and N is even, then the error term is :math:`B_N (\Delta x)^{N+3} f^{N+2}(\xi)`. Parameters ---------- rn : int The integer order for equally-spaced data or the relative positions of the samples with the first sample at 0 and the last at N, where N+1 is the length of `rn`. N is the order of the Newton-Cotes integration. equal : int, optional Set to 1 to enforce equally spaced data. Returns ------- an : ndarray 1-D array of weights to apply to the function at the provided sample positions. B : float Error coefficient. Notes ----- Normally, the Newton-Cotes rules are used on smaller integration regions and a composite rule is used to return the total integral. nquad(func, ranges, args=None, opts=None) Integration over multiple variables. Wraps `quad` to enable integration over multiple variables. Various options allow improved integration of discontinuous functions, as well as the use of weighted integration, and generally finer control of the integration process. Parameters ---------- func : callable The function to be integrated. Has arguments of ``x0, ... xn``, ``t0, tm``, where integration is carried out over ``x0, ... xn``, which must be floats. Function signature should be ``func(x0, x1, ..., xn, t0, t1, ..., tm)``. Integration is carried out in order. That is, integration over ``x0`` is the innermost integral, and ``xn`` is the outermost. ranges : iterable object Each element of ranges may be either a sequence of 2 numbers, or else a callable that returns such a sequence. ``ranges[0]`` corresponds to integration over x0, and so on. If an element of ranges is a callable, then it will be called with all of the integration arguments available. e.g. if ``func = f(x0, x1, x2)``, then ``ranges[0]`` may be defined as either ``(a, b)`` or else as ``(a, b) = range0(x1, x2)``. args : iterable object, optional Additional arguments ``t0, ..., tn``, required by `func`. opts : iterable object or dict, optional Options to be passed to `quad`. May be empty, a dict, or a sequence of dicts or functions that return a dict. If empty, the default options from scipy.integrate.quadare used. If a dict, the same options are used for all levels of integraion. If a sequence, then each element of the sequence corresponds to a particular integration. e.g. opts[0] corresponds to integration over x0, and so on. The available options together with their default values are: - epsabs = 1.49e-08 - epsrel = 1.49e-08 - limit = 50 - points = None - weight = None - wvar = None - wopts = None The ``full_output`` option from `quad` is unavailable, due to the complexity of handling the large amount of data such an option would return for this kind of nested integration. For more information on these options, see `quad` and `quad_explain`. Returns ------- result : float The result of the integration. abserr : float The maximum of the estimates of the absolute error in the various integration results. See Also -------- quad : 1-dimensional numerical integration dblquad, tplquad : double and triple integrals fixed_quad : fixed-order Gaussian quadrature quadrature : adaptive Gaussian quadrature Examples -------- >>> from scipy import integrate >>> func = lambda x0,x1,x2,x3 : x0**2 + x1*x2 - x3**3 + np.sin(x0) + ( ... 1 if (x0-.2*x3-.5-.25*x1>0) else 0) >>> points = [[lambda (x1,x2,x3) : 0.2*x3 + 0.5 + 0.25*x1], [], [], []] >>> def opts0(*args, **kwargs): ... return {'points':[0.2*args[2] + 0.5 + 0.25*args[0]]} >>> integrate.nquad(func, [[0,1], [-1,1], [.13,.8], [-.15,1]], ... opts=[opts0,{},{},{}]) (1.5267454070738633, 2.9437360001402324e-14) >>> scale = .1 >>> def func2(x0, x1, x2, x3, t0, t1): ... return x0*x1*x3**2 + np.sin(x2) + 1 + (1 if x0+t1*x1-t0>0 else 0) >>> def lim0(x1, x2, x3, t0, t1): ... return [scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) - 1, ... scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) + 1] >>> def lim1(x2, x3, t0, t1): ... return [scale * (t0*x2 + t1*x3) - 1, ... scale * (t0*x2 + t1*x3) + 1] >>> def lim2(x3, t0, t1): ... return [scale * (x3 + t0**2*t1**3) - 1, ... scale * (x3 + t0**2*t1**3) + 1] >>> def lim3(t0, t1): ... return [scale * (t0+t1) - 1, scale * (t0+t1) + 1] >>> def opts0(x1, x2, x3, t0, t1): ... return {'points' : [t0 - t1*x1]} >>> def opts1(x2, x3, t0, t1): ... return {} >>> def opts2(x3, t0, t1): ... return {} >>> def opts3(t0, t1): ... return {} >>> integrate.nquad(func2, [lim0, lim1, lim2, lim3], args=(0,0), opts=[opts0, opts1, opts2, opts3]) (25.066666666666666, 2.7829590483937256e-13) odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0) Integrate a system of ordinary differential equations. Solve a system of ordinary differential equations using lsoda from the FORTRAN library odepack. Solves the initial value problem for stiff or non-stiff systems of first order ode-s:: dy/dt = func(y,t0,...) where y can be a vector. Parameters ---------- func : callable(y, t0, ...) Computes the derivative of y at t0. y0 : array Initial condition on y (can be a vector). t : array A sequence of time points for which to solve for y. The initial value point should be the first element of this sequence. args : tuple, optional Extra arguments to pass to function. Dfun : callable(y, t0, ...) Gradient (Jacobian) of `func`. col_deriv : bool, optional True if `Dfun` defines derivatives down columns (faster), otherwise `Dfun` should define derivatives across rows. full_output : bool, optional True if to return a dictionary of optional outputs as the second output printmessg : bool, optional Whether to print the convergence message Returns ------- y : array, shape (len(t), len(y0)) Array containing the value of y for each desired time in t, with the initial value `y0` in the first row. infodict : dict, only returned if full_output == True Dictionary containing additional output information ======= ============================================================ key meaning ======= ============================================================ 'hu' vector of step sizes successfully used for each time step. 'tcur' vector with the value of t reached for each time step. (will always be at least as large as the input times). 'tolsf' vector of tolerance scale factors, greater than 1.0, computed when a request for too much accuracy was detected. 'tsw' value of t at the time of the last method switch (given for each time step) 'nst' cumulative number of time steps 'nfe' cumulative number of function evaluations for each time step 'nje' cumulative number of jacobian evaluations for each time step 'nqu' a vector of method orders for each successful step. 'imxer' index of the component of largest magnitude in the weighted local error vector (e / ewt) on an error return, -1 otherwise. 'lenrw' the length of the double work array required. 'leniw' the length of integer work array required. 'mused' a vector of method indicators for each successful time step: 1: adams (nonstiff), 2: bdf (stiff) ======= ============================================================ Other Parameters ---------------- ml, mu : int, optional If either of these are not None or non-negative, then the Jacobian is assumed to be banded. These give the number of lower and upper non-zero diagonals in this banded matrix. For the banded case, `Dfun` should return a matrix whose rows contain the non-zero bands (starting with the lowest diagonal). Thus, the return matrix `jac` from `Dfun` should have shape ``(ml + mu + 1, len(y0))`` when ``ml >=0`` or ``mu >=0``. The data in `jac` must be stored such that ``jac[i - j + mu, j]`` holds the derivative of the `i`th equation with respect to the `j`th state variable. If `col_deriv` is True, the transpose of this `jac` must be returned. rtol, atol : float, optional The input parameters `rtol` and `atol` determine the error control performed by the solver. The solver will control the vector, e, of estimated local errors in y, according to an inequality of the form ``max-norm of (e / ewt) <= 1``, where ewt is a vector of positive error weights computed as ``ewt = rtol * abs(y) + atol``. rtol and atol can be either vectors the same length as y or scalars. Defaults to 1.49012e-8. tcrit : ndarray, optional Vector of critical points (e.g. singularities) where integration care should be taken. h0 : float, (0: solver-determined), optional The step size to be attempted on the first step. hmax : float, (0: solver-determined), optional The maximum absolute step size allowed. hmin : float, (0: solver-determined), optional The minimum absolute step size allowed. ixpr : bool, optional Whether to generate extra printing at method switches. mxstep : int, (0: solver-determined), optional Maximum number of (internally defined) steps allowed for each integration point in t. mxhnil : int, (0: solver-determined), optional Maximum number of messages printed. mxordn : int, (0: solver-determined), optional Maximum order to be allowed for the non-stiff (Adams) method. mxords : int, (0: solver-determined), optional Maximum order to be allowed for the stiff (BDF) method. See Also -------- ode : a more object-oriented integrator based on VODE. quad : for finding the area under a curve. quad(func, a, b, args=(), full_output=0, epsabs=1.49e-08, epsrel=1.49e-08, limit=50, points=None, weight=None, wvar=None, wopts=None, maxp1=50, limlst=50) Compute a definite integral. Integrate func from `a` to `b` (possibly infinite interval) using a technique from the Fortran library QUADPACK. Parameters ---------- func : function A Python function or method to integrate. If `func` takes many arguments, it is integrated along the axis corresponding to the first argument. a : float Lower limit of integration (use -numpy.inf for -infinity). b : float Upper limit of integration (use numpy.inf for +infinity). args : tuple, optional Extra arguments to pass to `func`. full_output : int, optional Non-zero to return a dictionary of integration information. If non-zero, warning messages are also suppressed and the message is appended to the output tuple. Returns ------- y : float The integral of func from `a` to `b`. abserr : float An estimate of the absolute error in the result. infodict : dict A dictionary containing additional information. Run scipy.integrate.quad_explain() for more information. message : A convergence message. explain : Appended only with 'cos' or 'sin' weighting and infinite integration limits, it contains an explanation of the codes in infodict['ierlst'] Other Parameters ---------------- epsabs : float or int, optional Absolute error tolerance. epsrel : float or int, optional Relative error tolerance. limit : float or int, optional An upper bound on the number of subintervals used in the adaptive algorithm. points : (sequence of floats,ints), optional A sequence of break points in the bounded integration interval where local difficulties of the integrand may occur (e.g., singularities, discontinuities). The sequence does not have to be sorted. weight : float or int, optional String indicating weighting function. Full explanation for this and the remaining arguments can be found below. wvar : optional Variables for use with weighting functions. wopts : optional Optional input for reusing Chebyshev moments. maxp1 : float or int, optional An upper bound on the number of Chebyshev moments. limlst : int, optional Upper bound on the number of cycles (>=3) for use with a sinusoidal weighting and an infinite end-point. See Also -------- dblquad : double integral tplquad : triple integral nquad : n-dimensional integrals (uses `quad` recursively) fixed_quad : fixed-order Gaussian quadrature quadrature : adaptive Gaussian quadrature odeint : ODE integrator ode : ODE integrator simps : integrator for sampled data romb : integrator for sampled data scipy.special : for coefficients and roots of orthogonal polynomials Notes ----- **Extra information for quad() inputs and outputs** If full_output is non-zero, then the third output argument (infodict) is a dictionary with entries as tabulated below. For infinite limits, the range is transformed to (0,1) and the optional outputs are given with respect to this transformed range. Let M be the input argument limit and let K be infodict['last']. The entries are: 'neval' The number of function evaluations. 'last' The number, K, of subintervals produced in the subdivision process. 'alist' A rank-1 array of length M, the first K elements of which are the left end points of the subintervals in the partition of the integration range. 'blist' A rank-1 array of length M, the first K elements of which are the right end points of the subintervals. 'rlist' A rank-1 array of length M, the first K elements of which are the integral approximations on the subintervals. 'elist' A rank-1 array of length M, the first K elements of which are the moduli of the absolute error estimates on the subintervals. 'iord' A rank-1 integer array of length M, the first L elements of which are pointers to the error estimates over the subintervals with L=K if K<=M/2+2 or L=M+1-K otherwise. Let I be the sequence infodict['iord'] and let E be the sequence infodict['elist']. Then E[I[1]], ..., E[I[L]] forms a decreasing sequence. If the input argument points is provided (i.e. it is not None), the following additional outputs are placed in the output dictionary. Assume the points sequence is of length P. 'pts' A rank-1 array of length P+2 containing the integration limits and the break points of the intervals in ascending order. This is an array giving the subintervals over which integration will occur. 'level' A rank-1 integer array of length M (=limit), containing the subdivision levels of the subintervals, i.e., if (aa,bb) is a subinterval of (pts[1], pts[2]) where pts[0] and pts[2] are adjacent elements of infodict['pts'], then (aa,bb) has level l if |bb-aa|=|pts[2]-pts[1]| * 2**(-l). 'ndin' A rank-1 integer array of length P+2. After the first integration over the intervals (pts[1], pts[2]), the error estimates over some of the intervals may have been increased artificially in order to put their subdivision forward. This array has ones in slots corresponding to the subintervals for which this happens. **Weighting the integrand** The input variables, *weight* and *wvar*, are used to weight the integrand by a select list of functions. Different integration methods are used to compute the integral with these weighting functions. The possible values of weight and the corresponding weighting functions are. ========== =================================== ===================== ``weight`` Weight function used ``wvar`` ========== =================================== ===================== 'cos' cos(w*x) wvar = w 'sin' sin(w*x) wvar = w 'alg' g(x) = ((x-a)**alpha)*((b-x)**beta) wvar = (alpha, beta) 'alg-loga' g(x)*log(x-a) wvar = (alpha, beta) 'alg-logb' g(x)*log(b-x) wvar = (alpha, beta) 'alg-log' g(x)*log(x-a)*log(b-x) wvar = (alpha, beta) 'cauchy' 1/(x-c) wvar = c ========== =================================== ===================== wvar holds the parameter w, (alpha, beta), or c depending on the weight selected. In these expressions, a and b are the integration limits. For the 'cos' and 'sin' weighting, additional inputs and outputs are available. For finite integration limits, the integration is performed using a Clenshaw-Curtis method which uses Chebyshev moments. For repeated calculations, these moments are saved in the output dictionary: 'momcom' The maximum level of Chebyshev moments that have been computed, i.e., if M_c is infodict['momcom'] then the moments have been computed for intervals of length |b-a|* 2**(-l), l=0,1,...,M_c. 'nnlog' A rank-1 integer array of length M(=limit), containing the subdivision levels of the subintervals, i.e., an element of this array is equal to l if the corresponding subinterval is |b-a|* 2**(-l). 'chebmo' A rank-2 array of shape (25, maxp1) containing the computed Chebyshev moments. These can be passed on to an integration over the same interval by passing this array as the second element of the sequence wopts and passing infodict['momcom'] as the first element. If one of the integration limits is infinite, then a Fourier integral is computed (assuming w neq 0). If full_output is 1 and a numerical error is encountered, besides the error message attached to the output tuple, a dictionary is also appended to the output tuple which translates the error codes in the array info['ierlst'] to English messages. The output information dictionary contains the following entries instead of 'last', 'alist', 'blist', 'rlist', and 'elist': 'lst' The number of subintervals needed for the integration (call it K_f). 'rslst' A rank-1 array of length M_f=limlst, whose first K_f elements contain the integral contribution over the interval (a+(k-1)c, a+kc) where c = (2*floor(|w|) + 1) * pi / |w| and k=1,2,...,K_f. 'erlst' A rank-1 array of length M_f containing the error estimate corresponding to the interval in the same position in infodict['rslist']. 'ierlst' A rank-1 integer array of length M_f containing an error flag corresponding to the interval in the same position in infodict['rslist']. See the explanation dictionary (last entry in the output tuple) for the meaning of the codes. Examples -------- Calculate :math:`\int^4_0 x^2 dx` and compare with an analytic result >>> from scipy import integrate >>> x2 = lambda x: x**2 >>> integrate.quad(x2, 0, 4) (21.333333333333332, 2.3684757858670003e-13) >>> print(4**3 / 3.) # analytical result 21.3333333333 Calculate :math:`\int^\infty_0 e^{-x} dx` >>> invexp = lambda x: np.exp(-x) >>> integrate.quad(invexp, 0, np.inf) (1.0, 5.842605999138044e-11) >>> f = lambda x,a : a*x >>> y, err = integrate.quad(f, 0, 1, args=(1,)) >>> y 0.5 >>> y, err = integrate.quad(f, 0, 1, args=(3,)) >>> y 1.5 quad_explain(output=<IPython.kernel.zmq.iostream.OutStream object at 0x7f3f287e60f0>) Print extra information about integrate.quad() parameters and returns. Parameters ---------- output : instance with "write" method Information about `quad` is passed to ``output.write()``. Default is ``sys.stdout``. Returns ------- None quadrature(func, a, b, args=(), tol=1.49e-08, rtol=1.49e-08, maxiter=50, vec_func=True, miniter=1) Compute a definite integral using fixed-tolerance Gaussian quadrature. Integrate `func` from `a` to `b` using Gaussian quadrature with absolute tolerance `tol`. Parameters ---------- func : function A Python function or method to integrate. a : float Lower limit of integration. b : float Upper limit of integration. args : tuple, optional Extra arguments to pass to function. tol, rol : float, optional Iteration stops when error between last two iterates is less than `tol` OR the relative change is less than `rtol`. maxiter : int, optional Maximum order of Gaussian quadrature. vec_func : bool, optional True or False if func handles arrays as arguments (is a "vector" function). Default is True. miniter : int, optional Minimum order of Gaussian quadrature. Returns ------- val : float Gaussian quadrature approximation (within tolerance) to integral. err : float Difference between last two estimates of the integral. See also -------- romberg: adaptive Romberg quadrature fixed_quad: fixed-order Gaussian quadrature quad: adaptive quadrature using QUADPACK dblquad: double integrals tplquad: triple integrals romb: integrator for sampled data simps: integrator for sampled data cumtrapz: cumulative integration for sampled data ode: ODE integrator odeint: ODE integrator romb(y, dx=1.0, axis=-1, show=False) Romberg integration using samples of a function. Parameters ---------- y : array_like A vector of ``2**k + 1`` equally-spaced samples of a function. dx : array_like, optional The sample spacing. Default is 1. axis : int, optional The axis along which to integrate. Default is -1 (last axis). show : bool, optional When `y` is a single 1-D array, then if this argument is True print the table showing Richardson extrapolation from the samples. Default is False. Returns ------- romb : ndarray The integrated result for `axis`. See also -------- quad : adaptive quadrature using QUADPACK romberg : adaptive Romberg quadrature quadrature : adaptive Gaussian quadrature fixed_quad : fixed-order Gaussian quadrature dblquad : double integrals tplquad : triple integrals simps : integrators for sampled data cumtrapz : cumulative integration for sampled data ode : ODE integrators odeint : ODE integrators romberg(function, a, b, args=(), tol=1.48e-08, rtol=1.48e-08, show=False, divmax=10, vec_func=False) Romberg integration of a callable function or method. Returns the integral of `function` (a function of one variable) over the interval (`a`, `b`). If `show` is 1, the triangular array of the intermediate results will be printed. If `vec_func` is True (default is False), then `function` is assumed to support vector arguments. Parameters ---------- function : callable Function to be integrated. a : float Lower limit of integration. b : float Upper limit of integration. Returns ------- results : float Result of the integration. Other Parameters ---------------- args : tuple, optional Extra arguments to pass to function. Each element of `args` will be passed as a single argument to `func`. Default is to pass no extra arguments. tol, rtol : float, optional The desired absolute and relative tolerances. Defaults are 1.48e-8. show : bool, optional Whether to print the results. Default is False. divmax : int, optional Maximum order of extrapolation. Default is 10. vec_func : bool, optional Whether `func` handles arrays as arguments (i.e whether it is a "vector" function). Default is False. See Also -------- fixed_quad : Fixed-order Gaussian quadrature. quad : Adaptive quadrature using QUADPACK. dblquad : Double integrals. tplquad : Triple integrals. romb : Integrators for sampled data. simps : Integrators for sampled data. cumtrapz : Cumulative integration for sampled data. ode : ODE integrator. odeint : ODE integrator. References ---------- .. [1] 'Romberg's method' http://en.wikipedia.org/wiki/Romberg%27s_method Examples -------- Integrate a gaussian from 0 to 1 and compare to the error function. >>> from scipy import integrate >>> from scipy.special import erf >>> gaussian = lambda x: 1/np.sqrt(np.pi) * np.exp(-x**2) >>> result = integrate.romberg(gaussian, 0, 1, show=True) Romberg integration of <function vfunc at ...> from [0, 1] :: Steps StepSize Results 1 1.000000 0.385872 2 0.500000 0.412631 0.421551 4 0.250000 0.419184 0.421368 0.421356 8 0.125000 0.420810 0.421352 0.421350 0.421350 16 0.062500 0.421215 0.421350 0.421350 0.421350 0.421350 32 0.031250 0.421317 0.421350 0.421350 0.421350 0.421350 0.421350 The final result is 0.421350396475 after 33 function evaluations. >>> print("%g %g" % (2*result, erf(1))) 0.842701 0.842701 simps(y, x=None, dx=1, axis=-1, even='avg') Integrate y(x) using samples along the given axis and the composite Simpson's rule. If x is None, spacing of dx is assumed. If there are an even number of samples, N, then there are an odd number of intervals (N-1), but Simpson's rule requires an even number of intervals. The parameter 'even' controls how this is handled. Parameters ---------- y : array_like Array to be integrated. x : array_like, optional If given, the points at which `y` is sampled. dx : int, optional Spacing of integration points along axis of `y`. Only used when `x` is None. Default is 1. axis : int, optional Axis along which to integrate. Default is the last axis. even : {'avg', 'first', 'str'}, optional 'avg' : Average two results:1) use the first N-2 intervals with a trapezoidal rule on the last interval and 2) use the last N-2 intervals with a trapezoidal rule on the first interval. 'first' : Use Simpson's rule for the first N-2 intervals with a trapezoidal rule on the last interval. 'last' : Use Simpson's rule for the last N-2 intervals with a trapezoidal rule on the first interval. See Also -------- quad: adaptive quadrature using QUADPACK romberg: adaptive Romberg quadrature quadrature: adaptive Gaussian quadrature fixed_quad: fixed-order Gaussian quadrature dblquad: double integrals tplquad: triple integrals romb: integrators for sampled data cumtrapz: cumulative integration for sampled data ode: ODE integrators odeint: ODE integrators Notes ----- For an odd number of samples that are equally spaced the result is exact if the function is a polynomial of order 3 or less. If the samples are not equally spaced, then the result is exact only if the function is a polynomial of order 2 or less. tplquad(func, a, b, gfun, hfun, qfun, rfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08) Compute a triple (definite) integral. Return the triple integral of ``func(z, y, x)`` from ``x = a..b``, ``y = gfun(x)..hfun(x)``, and ``z = qfun(x,y)..rfun(x,y)``. Parameters ---------- func : function A Python function or method of at least three variables in the order (z, y, x). (a,b) : tuple The limits of integration in x: `a` < `b` gfun : function The lower boundary curve in y which is a function taking a single floating point argument (x) and returning a floating point result: a lambda function can be useful here. hfun : function The upper boundary curve in y (same requirements as `gfun`). qfun : function The lower boundary surface in z. It must be a function that takes two floats in the order (x, y) and returns a float. rfun : function The upper boundary surface in z. (Same requirements as `qfun`.) args : Arguments Extra arguments to pass to `func`. epsabs : float, optional Absolute tolerance passed directly to the innermost 1-D quadrature integration. Default is 1.49e-8. epsrel : float, optional Relative tolerance of the innermost 1-D integrals. Default is 1.49e-8. Returns ------- y : float The resultant integral. abserr : float An estimate of the error. See Also -------- quad: Adaptive quadrature using QUADPACK quadrature: Adaptive Gaussian quadrature fixed_quad: Fixed-order Gaussian quadrature dblquad: Double integrals nquad : N-dimensional integrals romb: Integrators for sampled data simps: Integrators for sampled data ode: ODE integrators odeint: ODE integrators scipy.special: For coefficients and roots of orthogonal polynomials trapz(y, x=None, dx=1.0, axis=-1) Integrate along the given axis using the composite trapezoidal rule. Integrate `y` (`x`) along given axis. Parameters ---------- y : array_like Input array to integrate. x : array_like, optional If `x` is None, then spacing between all `y` elements is `dx`. dx : scalar, optional If `x` is None, spacing given by `dx` is assumed. Default is 1. axis : int, optional Specify the axis. Returns ------- trapz : float Definite integral as approximated by trapezoidal rule. See Also -------- sum, cumsum Notes ----- Image [2]_ illustrates trapezoidal rule -- y-axis locations of points will be taken from `y` array, by default x-axis distances between points will be 1.0, alternatively they can be provided with `x` array or with `dx` scalar. Return value will be equal to combined area under the red lines. References ---------- .. [1] Wikipedia page: http://en.wikipedia.org/wiki/Trapezoidal_rule .. [2] Illustration image: http://en.wikipedia.org/wiki/File:Composite_trapezoidal_rule_illustration.png Examples -------- >>> np.trapz([1,2,3]) 4.0 >>> np.trapz([1,2,3], x=[4,6,8]) 8.0 >>> np.trapz([1,2,3], dx=2) 8.0 >>> a = np.arange(6).reshape(2, 3) >>> a array([[0, 1, 2], [3, 4, 5]]) >>> np.trapz(a, axis=0) array([ 1.5, 2.5, 3.5]) >>> np.trapz(a, axis=1) array([ 2., 8.]) DATA __all__ = ['IntegrationWarning', 'absolute_import', 'complex_ode', 'cu... absolute_import = _Feature((2, 5, 0, 'alpha', 1), (3, 0, 0, 'alpha', 0... division = _Feature((2, 2, 0, 'alpha', 2), (3, 0, 0, 'alpha', 0), 8192... print_function = _Feature((2, 6, 0, 'alpha', 2), (3, 0, 0, 'alpha', 0)... FILE /home/alex/anaconda3/lib/python3.4/site-packages/scipy/integrate/__init__.py
Como se puede ver en la ayuda, si queremos realizar una integración numérica de una función de una variable, debemos utilizar quad
(aunque también podemos usar trapz
, simps
... La forma de acceder a ella tal y como hemos importado el paquete sería ejecutando integrate.quad
. Sin emabrgo, sería más normal importar del siguiete modo:
from scipy.integrate import quad
De este modo se puede usar la función quad, simplemente como quad
. Pero todavía no sabemos cómo funciona, ¿te atreves a investigarlo tú?
help(quad)
Help on function quad in module scipy.integrate.quadpack: quad(func, a, b, args=(), full_output=0, epsabs=1.49e-08, epsrel=1.49e-08, limit=50, points=None, weight=None, wvar=None, wopts=None, maxp1=50, limlst=50) Compute a definite integral. Integrate func from `a` to `b` (possibly infinite interval) using a technique from the Fortran library QUADPACK. Parameters ---------- func : function A Python function or method to integrate. If `func` takes many arguments, it is integrated along the axis corresponding to the first argument. a : float Lower limit of integration (use -numpy.inf for -infinity). b : float Upper limit of integration (use numpy.inf for +infinity). args : tuple, optional Extra arguments to pass to `func`. full_output : int, optional Non-zero to return a dictionary of integration information. If non-zero, warning messages are also suppressed and the message is appended to the output tuple. Returns ------- y : float The integral of func from `a` to `b`. abserr : float An estimate of the absolute error in the result. infodict : dict A dictionary containing additional information. Run scipy.integrate.quad_explain() for more information. message : A convergence message. explain : Appended only with 'cos' or 'sin' weighting and infinite integration limits, it contains an explanation of the codes in infodict['ierlst'] Other Parameters ---------------- epsabs : float or int, optional Absolute error tolerance. epsrel : float or int, optional Relative error tolerance. limit : float or int, optional An upper bound on the number of subintervals used in the adaptive algorithm. points : (sequence of floats,ints), optional A sequence of break points in the bounded integration interval where local difficulties of the integrand may occur (e.g., singularities, discontinuities). The sequence does not have to be sorted. weight : float or int, optional String indicating weighting function. Full explanation for this and the remaining arguments can be found below. wvar : optional Variables for use with weighting functions. wopts : optional Optional input for reusing Chebyshev moments. maxp1 : float or int, optional An upper bound on the number of Chebyshev moments. limlst : int, optional Upper bound on the number of cycles (>=3) for use with a sinusoidal weighting and an infinite end-point. See Also -------- dblquad : double integral tplquad : triple integral nquad : n-dimensional integrals (uses `quad` recursively) fixed_quad : fixed-order Gaussian quadrature quadrature : adaptive Gaussian quadrature odeint : ODE integrator ode : ODE integrator simps : integrator for sampled data romb : integrator for sampled data scipy.special : for coefficients and roots of orthogonal polynomials Notes ----- **Extra information for quad() inputs and outputs** If full_output is non-zero, then the third output argument (infodict) is a dictionary with entries as tabulated below. For infinite limits, the range is transformed to (0,1) and the optional outputs are given with respect to this transformed range. Let M be the input argument limit and let K be infodict['last']. The entries are: 'neval' The number of function evaluations. 'last' The number, K, of subintervals produced in the subdivision process. 'alist' A rank-1 array of length M, the first K elements of which are the left end points of the subintervals in the partition of the integration range. 'blist' A rank-1 array of length M, the first K elements of which are the right end points of the subintervals. 'rlist' A rank-1 array of length M, the first K elements of which are the integral approximations on the subintervals. 'elist' A rank-1 array of length M, the first K elements of which are the moduli of the absolute error estimates on the subintervals. 'iord' A rank-1 integer array of length M, the first L elements of which are pointers to the error estimates over the subintervals with L=K if K<=M/2+2 or L=M+1-K otherwise. Let I be the sequence infodict['iord'] and let E be the sequence infodict['elist']. Then E[I[1]], ..., E[I[L]] forms a decreasing sequence. If the input argument points is provided (i.e. it is not None), the following additional outputs are placed in the output dictionary. Assume the points sequence is of length P. 'pts' A rank-1 array of length P+2 containing the integration limits and the break points of the intervals in ascending order. This is an array giving the subintervals over which integration will occur. 'level' A rank-1 integer array of length M (=limit), containing the subdivision levels of the subintervals, i.e., if (aa,bb) is a subinterval of (pts[1], pts[2]) where pts[0] and pts[2] are adjacent elements of infodict['pts'], then (aa,bb) has level l if |bb-aa|=|pts[2]-pts[1]| * 2**(-l). 'ndin' A rank-1 integer array of length P+2. After the first integration over the intervals (pts[1], pts[2]), the error estimates over some of the intervals may have been increased artificially in order to put their subdivision forward. This array has ones in slots corresponding to the subintervals for which this happens. **Weighting the integrand** The input variables, *weight* and *wvar*, are used to weight the integrand by a select list of functions. Different integration methods are used to compute the integral with these weighting functions. The possible values of weight and the corresponding weighting functions are. ========== =================================== ===================== ``weight`` Weight function used ``wvar`` ========== =================================== ===================== 'cos' cos(w*x) wvar = w 'sin' sin(w*x) wvar = w 'alg' g(x) = ((x-a)**alpha)*((b-x)**beta) wvar = (alpha, beta) 'alg-loga' g(x)*log(x-a) wvar = (alpha, beta) 'alg-logb' g(x)*log(b-x) wvar = (alpha, beta) 'alg-log' g(x)*log(x-a)*log(b-x) wvar = (alpha, beta) 'cauchy' 1/(x-c) wvar = c ========== =================================== ===================== wvar holds the parameter w, (alpha, beta), or c depending on the weight selected. In these expressions, a and b are the integration limits. For the 'cos' and 'sin' weighting, additional inputs and outputs are available. For finite integration limits, the integration is performed using a Clenshaw-Curtis method which uses Chebyshev moments. For repeated calculations, these moments are saved in the output dictionary: 'momcom' The maximum level of Chebyshev moments that have been computed, i.e., if M_c is infodict['momcom'] then the moments have been computed for intervals of length |b-a|* 2**(-l), l=0,1,...,M_c. 'nnlog' A rank-1 integer array of length M(=limit), containing the subdivision levels of the subintervals, i.e., an element of this array is equal to l if the corresponding subinterval is |b-a|* 2**(-l). 'chebmo' A rank-2 array of shape (25, maxp1) containing the computed Chebyshev moments. These can be passed on to an integration over the same interval by passing this array as the second element of the sequence wopts and passing infodict['momcom'] as the first element. If one of the integration limits is infinite, then a Fourier integral is computed (assuming w neq 0). If full_output is 1 and a numerical error is encountered, besides the error message attached to the output tuple, a dictionary is also appended to the output tuple which translates the error codes in the array info['ierlst'] to English messages. The output information dictionary contains the following entries instead of 'last', 'alist', 'blist', 'rlist', and 'elist': 'lst' The number of subintervals needed for the integration (call it K_f). 'rslst' A rank-1 array of length M_f=limlst, whose first K_f elements contain the integral contribution over the interval (a+(k-1)c, a+kc) where c = (2*floor(|w|) + 1) * pi / |w| and k=1,2,...,K_f. 'erlst' A rank-1 array of length M_f containing the error estimate corresponding to the interval in the same position in infodict['rslist']. 'ierlst' A rank-1 integer array of length M_f containing an error flag corresponding to the interval in the same position in infodict['rslist']. See the explanation dictionary (last entry in the output tuple) for the meaning of the codes. Examples -------- Calculate :math:`\int^4_0 x^2 dx` and compare with an analytic result >>> from scipy import integrate >>> x2 = lambda x: x**2 >>> integrate.quad(x2, 0, 4) (21.333333333333332, 2.3684757858670003e-13) >>> print(4**3 / 3.) # analytical result 21.3333333333 Calculate :math:`\int^\infty_0 e^{-x} dx` >>> invexp = lambda x: np.exp(-x) >>> integrate.quad(invexp, 0, np.inf) (1.0, 5.842605999138044e-11) >>> f = lambda x,a : a*x >>> y, err = integrate.quad(f, 0, 1, args=(1,)) >>> y 0.5 >>> y, err = integrate.quad(f, 0, 1, args=(3,)) >>> y 1.5
Quizá esta ayuda te resulte más atractiva.
¿Qué es lo primero que necesitamos hacer para integrar una función? Pues sí, la función... definamos una:
$$f(x) = x \cdot sin(x)$$def fun(x):
return x * np.sin(x)
Antes de integrarla genera esta gráfica:
x = np.linspace(0,10,100)
y = fun(x)
plt.title('$y = x sin(x)$', fontsize = 25)
plt.plot(x,y, linewidth = 2)
x_fill = np.linspace(2,9,100)
y_fill = fun(x_fill)
plt.fill_between(x_fill, y_fill, color='gray', alpha=0.5)
plt.grid()
quad
¶Integremos la función en el intervalo $[2, 9]$. Recuerda que esto te calcula la integral, no el área:
value, err = quad(fun, 2, 9)
print("El resultado es: ", value, "con un error de: ", err)
El resultado es: 6.870699742283883 con un error de: 2.864870105641461e-13
Según figura en la documentación a estos métodos hay que pasarles las coordenadas de los puntos (no la función). Esto puede ser útil si no disponemos de una función, sino de una serie da valores, que por ejemplo, provienen de un experimento.
x = np.linspace(2,9,100)
value = integrate.trapz(fun(x), x)
print("El resultado es: ", value)
El resultado es: 6.86742266171
x = np.linspace(2,9,100)
value = integrate.simps(fun(x), x)
print("El resultado es: ", value)
El resultado es: 6.8705759095
Este módulo contiene un gran número de distribuciones de probabilidad, tanto continuas como discretas, así como un creciente número de funciones estadísticas.
# Importar el módulo entero
import scipy.stats as st
# Información
np.info(st)
========================================== Statistical functions (:mod:`scipy.stats`) ========================================== .. module:: scipy.stats This module contains a large number of probability distributions as well as a growing library of statistical functions. Each included distribution is an instance of the class rv_continous: For each given name the following methods are available: .. autosummary:: :toctree: generated/ rv_continuous rv_continuous.pdf rv_continuous.logpdf rv_continuous.cdf rv_continuous.logcdf rv_continuous.sf rv_continuous.logsf rv_continuous.ppf rv_continuous.isf rv_continuous.moment rv_continuous.stats rv_continuous.entropy rv_continuous.fit rv_continuous.expect Calling the instance as a function returns a frozen pdf whose shape, location, and scale parameters are fixed. Similarly, each discrete distribution is an instance of the class rv_discrete: .. autosummary:: :toctree: generated/ rv_discrete rv_discrete.rvs rv_discrete.pmf rv_discrete.logpmf rv_discrete.cdf rv_discrete.logcdf rv_discrete.sf rv_discrete.logsf rv_discrete.ppf rv_discrete.isf rv_discrete.stats rv_discrete.moment rv_discrete.entropy rv_discrete.expect Continuous distributions ======================== .. autosummary:: :toctree: generated/ alpha -- Alpha anglit -- Anglit arcsine -- Arcsine beta -- Beta betaprime -- Beta Prime bradford -- Bradford burr -- Burr cauchy -- Cauchy chi -- Chi chi2 -- Chi-squared cosine -- Cosine dgamma -- Double Gamma dweibull -- Double Weibull erlang -- Erlang expon -- Exponential exponweib -- Exponentiated Weibull exponpow -- Exponential Power f -- F (Snecdor F) fatiguelife -- Fatigue Life (Birnbaum-Sanders) fisk -- Fisk foldcauchy -- Folded Cauchy foldnorm -- Folded Normal frechet_r -- Frechet Right Sided, Extreme Value Type II (Extreme LB) or weibull_min frechet_l -- Frechet Left Sided, Weibull_max genlogistic -- Generalized Logistic genpareto -- Generalized Pareto genexpon -- Generalized Exponential genextreme -- Generalized Extreme Value gausshyper -- Gauss Hypergeometric gamma -- Gamma gengamma -- Generalized gamma genhalflogistic -- Generalized Half Logistic gilbrat -- Gilbrat gompertz -- Gompertz (Truncated Gumbel) gumbel_r -- Right Sided Gumbel, Log-Weibull, Fisher-Tippett, Extreme Value Type I gumbel_l -- Left Sided Gumbel, etc. halfcauchy -- Half Cauchy halflogistic -- Half Logistic halfnorm -- Half Normal hypsecant -- Hyperbolic Secant invgamma -- Inverse Gamma invgauss -- Inverse Gaussian invweibull -- Inverse Weibull johnsonsb -- Johnson SB johnsonsu -- Johnson SU ksone -- Kolmogorov-Smirnov one-sided (no stats) kstwobign -- Kolmogorov-Smirnov two-sided test for Large N (no stats) laplace -- Laplace logistic -- Logistic loggamma -- Log-Gamma loglaplace -- Log-Laplace (Log Double Exponential) lognorm -- Log-Normal lomax -- Lomax (Pareto of the second kind) maxwell -- Maxwell mielke -- Mielke's Beta-Kappa nakagami -- Nakagami ncx2 -- Non-central chi-squared ncf -- Non-central F nct -- Non-central Student's T norm -- Normal (Gaussian) pareto -- Pareto pearson3 -- Pearson type III powerlaw -- Power-function powerlognorm -- Power log normal powernorm -- Power normal rdist -- R-distribution reciprocal -- Reciprocal rayleigh -- Rayleigh rice -- Rice recipinvgauss -- Reciprocal Inverse Gaussian semicircular -- Semicircular t -- Student's T triang -- Triangular truncexpon -- Truncated Exponential truncnorm -- Truncated Normal tukeylambda -- Tukey-Lambda uniform -- Uniform vonmises -- Von-Mises (Circular) wald -- Wald weibull_min -- Minimum Weibull (see Frechet) weibull_max -- Maximum Weibull (see Frechet) wrapcauchy -- Wrapped Cauchy Multivariate distributions ========================== .. autosummary:: :toctree: generated/ multivariate_normal -- Multivariate normal distribution Discrete distributions ====================== .. autosummary:: :toctree: generated/ bernoulli -- Bernoulli binom -- Binomial boltzmann -- Boltzmann (Truncated Discrete Exponential) dlaplace -- Discrete Laplacian geom -- Geometric hypergeom -- Hypergeometric logser -- Logarithmic (Log-Series, Series) nbinom -- Negative Binomial planck -- Planck (Discrete Exponential) poisson -- Poisson randint -- Discrete Uniform skellam -- Skellam zipf -- Zipf Statistical functions ===================== Several of these functions have a similar version in scipy.stats.mstats which work for masked arrays. .. autosummary:: :toctree: generated/ describe -- Descriptive statistics gmean -- Geometric mean hmean -- Harmonic mean kurtosis -- Fisher or Pearson kurtosis kurtosistest -- mode -- Modal value moment -- Central moment normaltest -- skew -- Skewness skewtest -- tmean -- Truncated arithmetic mean tvar -- Truncated variance tmin -- tmax -- tstd -- tsem -- nanmean -- Mean, ignoring NaN values nanstd -- Standard deviation, ignoring NaN values nanmedian -- Median, ignoring NaN values variation -- Coefficient of variation .. autosummary:: :toctree: generated/ cumfreq _ histogram2 _ histogram _ itemfreq _ percentileofscore _ scoreatpercentile _ relfreq _ .. autosummary:: :toctree: generated/ binned_statistic -- Compute a binned statistic for a set of data. binned_statistic_2d -- Compute a 2-D binned statistic for a set of data. binned_statistic_dd -- Compute a d-D binned statistic for a set of data. .. autosummary:: :toctree: generated/ obrientransform signaltonoise bayes_mvs sem zmap zscore .. autosummary:: :toctree: generated/ sigmaclip threshold trimboth trim1 .. autosummary:: :toctree: generated/ f_oneway pearsonr spearmanr pointbiserialr kendalltau linregress .. autosummary:: :toctree: generated/ ttest_1samp ttest_ind ttest_rel kstest chisquare power_divergence ks_2samp mannwhitneyu tiecorrect rankdata ranksums wilcoxon kruskal friedmanchisquare .. autosummary:: :toctree: generated/ ansari bartlett levene shapiro anderson anderson_ksamp binom_test fligner mood .. autosummary:: :toctree: generated/ boxcox boxcox_normmax boxcox_llf entropy Contingency table functions =========================== .. autosummary:: :toctree: generated/ chi2_contingency contingency.expected_freq contingency.margins fisher_exact Plot-tests ========== .. autosummary:: :toctree: generated/ ppcc_max ppcc_plot probplot boxcox_normplot Masked statistics functions =========================== .. toctree:: stats.mstats Univariate and multivariate kernel density estimation (:mod:`scipy.stats.kde`) ============================================================================== .. autosummary:: :toctree: generated/ gaussian_kde For many more stat related functions install the software R and the interface package rpy.
Carguemos unos datos, por ejemplo unas notas de la carrera, y veamos cómo podemos aprovechar las funciones de scipy.stats
.
# esta línea no funciona en Windows
!head ../static/notas.csv
# Leemos el archivo
datos = np.loadtxt("../static/notas.csv", skiprows=1)
datos
array([ 2.9, 4.3, 3.9, 0. , 4.1, 7.3, 2.3, 5.6, 2.9, 3.9, 4.6, 6.3, 2.1, 2.1, 6.5, 1.9, 0. , 6.5, 2.5, 5.1, 5.3, 6.3, 5.4, 5.3, 5.3, 2. , 3.5, 4.4, 5.5, 3.6, 3.9, 2.5, 4.1, 3. , 4.6, 4. , 6.3, 0.6, 2.4, 6.5, 2.3, 4.6, 6.9, 5.1, 5.4, 5.3, 4.5, 6.5, 2.1, 5.5, 3.4, 8.1, 4. , 1.9, 1.6, 4.3, 4.6, 5.4, 1. , 6.5, 5.5, 4.9, 4. , 5.3, 3.5, 4.4, 2.8, 5.4, 3.5, 2.3, 4.8, 2.1, 6.6, 0.5, 2.1, 3.1, 3.4, 5.9, 3.4, 4.3, 1.5, 5.5, 4.4, 1.9, 4.4, 2.9, 3.9, 5.8, 2.8, 3. , 1.5, 2.6, 2.9, 3.4, 5.4, 3.6, 4.6, 5. , 1.4, 4.3, 4.6, 3.1, 2. , 3.6, 4. , 2.5, 3. , 5.1, 6.4, 3.5, 5.8, 4.1, 5.9, 4. , 6.4, 2.3, 7. , 1.4, 3.5, 4.4, 2.9, 5.1, 3.4, 4.8, 4.6, 4.3, 6.9, 5.4, 4. , 3.3, 1.4, 1.9, 3.8, 3.4, 3.6, 3.8, 6.3, 4.8, 4. , 6.8, 4. , 3.6, 4.4, 4.1, 6. , 4.1, 5.6, 3.9, 4.6, 5. , 6.5, 3.5, 5.5, 4.6, 4.8, 4.6, 6.5, 4.1, 4.4, 5.3, 3.6, 7.1, 4.6, 2.1, 3.3, 3.9, 4. , 4.4, 0.9, 4.3, 2.4, 2.9, 3.6, 1.4, 2.8, 2.5, 6.6, 0. , 5.1, 0. , 5.1, 4. , 2.6, 4.1, 4.6, 3.1, 4.4, 2.8, 2.8, 5.6, 3.9, 4.4, 4.1, 0.3, 2.4, 3.3, 2.5, 4.3, 2.5, 4.5, 4.8, 4.3, 3.3, 6. , 2.3, 4.5, 3.4, 4.5, 2.3, 4.5, 2.5, 6.4, 7. , 5.8, 3.4, 5.1, 4.9, 8.5, 3. , 3.3, 3.3, 2.6, 1.8, 2.9, 4.3, 2.1, 4.8, 5.3, 3.1, 5. , 4.6, 3.1, 5.6, 5.5, 4.3, 4.5, 5.4, 4.3, 1.5, 3.9, 5.4, 1.8, 3.5, 5.3, 2.6, 3.5, 3. , 4.9, 3.5, 3.4, 2. , 4.5, 5.3, 3. , 1.6, 4.8, 3.8, 0.6, 4.1, 5.4, 7.5, 0.6, 2.6, 2.1, 3.8, 3.3, 5.4, 3.6, 2.4, 4.6, 0.9, 5.8, 4.4, 1. , 3.6, 6.3, 4.4, 7.5, 5.9, 0.5, 4.3, 2.4, 6. , 4.6, 5.1, 6.1, 4. , 4.5, 7.9, 3.5, 3.1, 5. , 6.3, 4.9, 5.9, 4.6, 5.8, 5.4, 0.1, 1.8, 5.1, 4. , 2.4, 4.6, 4.9, 3.1, 1.4, 4. , 3.6, 4.3, 3.8, 4.4, 4.8, 5.1, 4. , 0. , 2.1, 5.9, 6.3, 3.1, 6. , 3.4, 1.9, 5.6, 5.3, 4.8, 2.6, 5.6, 4.8, 5.4, 3.4, 5.3, 4.1, 3.8, 3.6, 3.9, 2. , 3.5, 4. , 3.6, 0.6, 2.4, 3.9, 4.1, 2.8, 3. , 5.1, 5.4, 3.9, 3.3, 3.8, 2.5, 0.6, 2.8, 2.9, 7. , 6. , 2.8, 4. , 4.9, 4.8, 2.8, 2. , 4. , 2.6, 3.1, 2.9, 6.5, 4.3, 2.1, 3.9, 4.3, 0. , 7.4, 3.9])
# Descripción rápida de los datos
st.describe(datos)
(375, (0.0, 8.5), 3.9706666666666668, 2.5927736185383243, -0.13203546994646295, -0.06611485627230884)
# Histograma con st
st.histogram(datos, numbins=10, defaultlimits=(0,10))
(array([ 17., 20., 61., 78., 98., 61., 29., 9., 2., 0.]), 0, 1.0, 0)
# Pintemos un histograma con plt
plt.hist(datos, range(0,11,))
plt.xticks(range(0,11))
plt.grid(True)
plt.vlines(5, 0, 100, lw=5, colors='red', alpha=0.8)
plt.fill_between([0, 5], [100, 100], color='red', alpha=0.5)
<matplotlib.collections.PolyCollection at 0x7f3f07891358>
# Pintemos un histograma acumulado con plt
plt.hist(datos, range(0,11), cumulative=True)
plt.xticks(range(0,11))
plt.vlines(5, 0, 400, lw=5, colors='red', alpha=0.8)
plt.fill_between([0, 5], [400, 400], color='red', alpha=0.5)
plt.grid(True)
# Percentil
st.percentileofscore(datos, 5)
73.733333333333334
# Nota de un percentil
st.scoreatpercentile(datos, 50)
4.0
¿Te parecen normales estas notas? No, no me refiero a si te gustan o no... Me refiero a que si crees que estas notas se distribuyen de manera gaussiana.
# Parámetros
med = st.nanmean(datos)
des_tip = st.nanstd(datos)
# Distribución normal
dist_normal = st.norm(loc=med, scale=des_tip)
Ahora podemos ver:
pdf
cdf
De esta manera, nos ahorramos definir funciones como:
$$N(\mu,\sigma)=\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2}$$$$\phi(x)=\intop_{-\infty}^{x}\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}dx$$Para calcular probabilidades $P[a\leq X\leq b]=\intop_{a}^{b}f(x)dx$
# Calculamos la pdf
x = np.linspace(0, 10, 100)
y1 = dist_normal.pdf(x)
y1
array([ 1.18466112e-02, 1.38013272e-02, 1.60154276e-02, 1.85117381e-02, 2.13131100e-02, 2.44420387e-02, 2.79202305e-02, 3.17681221e-02, 3.60043569e-02, 4.06452262e-02, 4.57040834e-02, 5.11907430e-02, 5.71108776e-02, 6.34654270e-02, 7.02500358e-02, 7.74545375e-02, 8.50625018e-02, 9.30508640e-02, 1.01389654e-01, 1.10041839e-01, 1.18963303e-01, 1.28102958e-01, 1.37403018e-01, 1.46799426e-01, 1.56222440e-01, 1.65597374e-01, 1.74845495e-01, 1.83885049e-01, 1.92632412e-01, 2.01003342e-01, 2.08914298e-01, 2.16283816e-01, 2.23033887e-01, 2.29091330e-01, 2.34389108e-01, 2.38867558e-01, 2.42475512e-01, 2.45171267e-01, 2.46923389e-01, 2.47711323e-01, 2.47525793e-01, 2.46368984e-01, 2.44254502e-01, 2.41207102e-01, 2.37262214e-01, 2.32465247e-01, 2.26870729e-01, 2.20541268e-01, 2.13546393e-01, 2.05961280e-01, 1.97865419e-01, 1.89341226e-01, 1.80472670e-01, 1.71343910e-01, 1.62038000e-01, 1.52635673e-01, 1.43214236e-01, 1.33846588e-01, 1.24600387e-01, 1.15537363e-01, 1.06712791e-01, 9.81751290e-02, 8.99658014e-02, 8.21191408e-02, 7.46624629e-02, 6.76162683e-02, 6.09945536e-02, 5.48052151e-02, 4.90505281e-02, 4.37276824e-02, 3.88293574e-02, 3.43443198e-02, 3.02580270e-02, 2.65532243e-02, 2.32105212e-02, 2.02089381e-02, 1.75264140e-02, 1.51402697e-02, 1.30276213e-02, 1.11657426e-02, 9.53237395e-03, 8.10597949e-03, 6.86595419e-03, 5.79278319e-03, 4.86815747e-03, 4.07504992e-03, 3.39775625e-03, 2.82190558e-03, 2.33444533e-03, 1.92360505e-03, 1.57884349e-03, 1.29078296e-03, 1.05113462e-03, 8.52617863e-04, 6.88876683e-04, 5.54395268e-04, 4.44414802e-04, 3.54853006e-04, 2.82227567e-04, 2.23584322e-04])
# La representamos
plt.plot(x, y1)
plt.grid(True)
# Calculamos la cdf
y2 = dist_normal.cdf(x)
y2
array([ 0.00683286, 0.00812612, 0.00962974, 0.01137106, 0.01337976, 0.01568777, 0.0183293 , 0.02134067, 0.02476018, 0.02862791, 0.03298542, 0.03787547, 0.04334158, 0.04942765, 0.05617737, 0.06363372, 0.07183834, 0.08083088, 0.0906483 , 0.1013242 , 0.11288805, 0.12536454, 0.13877285, 0.15312598, 0.1684302 , 0.18468445, 0.20187994, 0.21999975, 0.2390186 , 0.25890272, 0.27960985, 0.30108939, 0.32328269, 0.34612342, 0.36953818, 0.39344715, 0.41776493, 0.44240138, 0.46726269, 0.49225242, 0.5172726 , 0.54222491, 0.56701185, 0.59153783, 0.61571033, 0.63944089, 0.66264614, 0.68524861, 0.70717752, 0.72836942, 0.74876866, 0.76832781, 0.78700782, 0.80477817, 0.82161679, 0.83750988, 0.85245164, 0.86644386, 0.87949545, 0.89162183, 0.90284437, 0.91318964, 0.92268879, 0.93137677, 0.93929165, 0.94647392, 0.95296581, 0.95881062, 0.96405221, 0.96873436, 0.97290038, 0.9765926 , 0.97985205, 0.98271818, 0.98522854, 0.98741866, 0.98932189, 0.99096933, 0.99238974, 0.99360961, 0.99465314, 0.99554231, 0.99629699, 0.99693499, 0.99747225, 0.99792288, 0.99829939, 0.99861271, 0.99887244, 0.9990869 , 0.99926327, 0.99940776, 0.99952567, 0.9996215 , 0.99969908, 0.99976165, 0.99981191, 0.99985212, 0.99988416, 0.9999096 ])
# La representamos
plt.plot(x, y2)
plt.grid(True)
Del mismo modo se pueden usar otras distribuciones continuas o discretas e incluso, definir distribuciones propias. Pero sigamos con las notas...
Ahora que ya hemos visualizado la distribución de las notas y que sabemos generar distribuciones normales. ¿Por qué no hacemos un test de Kolmogórov-Smirnov?
Se trata de ver lo bien o lo mal que se ajusta la distribución a una normal con $\mu=3.97$ y $\sigma²=2.57$
bars = st.histogram(datos, numbins=10, defaultlimits=(0,10))[0]
bars /= 375
plt.bar(np.arange(0,10), bars, alpha=0.5, width=1)
plt.plot(x, y1, c='black', lw=2)
plt.grid(True)
plt.hist(datos, np.linspace(0,10,51), cumulative=True, alpha=0.5)
plt.plot(x, y2 * 375, lw=2, c='black')
plt.xticks(range(0,11))
plt.grid(True)
datos2 = dist_normal.cdf
st.kstest(datos, dist_normal.cdf)
(0.04783071674813294, 0.34838712365988389)
Se rechaza la hipótesis nula si el valor p asociado al resultado observado es igual o menor que el nivel de significación establecido, convencionalmente 0,05 ó 0,01. Es decir, el valor p nos muestra la probabilidad de haber obtenido el resultado que hemos obtenido si suponemos que la hipótesis nula es cierta. (Wikipedia)
Si probamos con st.normaltest
que también comprueba la bondad del ajuste obtenemos un valor-p más alto:
st.normaltest(datos)
(1.1315306710610515, 0.56792532695583819)
En definitiva, parece que las notas esta vez siguieron una normal con $\mu=3.97$ y $\sigma²=2.57$
Las siguientes celdas contienen configuración del Notebook
Para visualizar y utlizar los enlaces a Twitter el notebook debe ejecutarse como seguro
File > Trusted Notebook
%%html
<a href="https://twitter.com/Pybonacci" class="twitter-follow-button" data-show-count="false">Follow @Pybonacci</a>
<script>!function(d,s,id){var js,fjs=d.getElementsByTagName(s)[0],p=/^http:/.test(d.location)?'http':'https';if(!d.getElementById(id)){js=d.createElement(s);js.id=id;js.src=p+'://platform.twitter.com/widgets.js';fjs.parentNode.insertBefore(js,fjs);}}(document, 'script', 'twitter-wjs');</script>
# Esta celda da el estilo al notebook
from IPython.core.display import HTML
css_file = '../static/styles/style.css'
HTML(open(css_file, "r").read())