B dGl@s|dZddlZddlmZmZddddd gZdd d Zdd dddZdd dddZ dd dddZ dddZ ddd Z dS)za This module contains functions for computing robust statistics using Tukey's biweight function. N) _expand_dimsmedian_absolute_deviationbiweight_locationbiweight_scalebiweight_midvariancebiweight_midcovariancebiweight_midcorrelationFcCsFt|tjjr tjj}tjj}n|r2tj}tj}n tj}tj}||fS)N) isinstancenpma MaskedArraymediansumZ nanmedianZnansum)data ignore_nan median_funcsum_funcrc/work/yifan.wang/ringdown/master-ringdown-env/lib/python3.7/site-packages/astropy/stats/biweight.py_stat_functionss r@)rc Cst||d\}}t|tjjr:|r:tjjt||dd}t|tj }|dkr`|||d}|dk rtt ||d}||}t |||d}|dkr|dkst|r|S|dk rt ||d}tj ddd |||} WdQRXtj dd t | d k} WdQRXd | d d } d | | <tj ddd l|||| |d|| |d} t| rd| Stj} t|tjjrtjj} | |d k|| SQRXdS)a# Compute the biweight location. The biweight location is a robust statistic for determining the central location of a distribution. It is given by: .. math:: \zeta_{biloc}= M + \frac{\sum_{|u_i|<1} \ (x_i - M) (1 - u_i^2)^2} {\sum_{|u_i|<1} \ (1 - u_i^2)^2} where :math:`x` is the input data, :math:`M` is the sample median (or the input initial location guess) and :math:`u_i` is given by: .. math:: u_{i} = \frac{(x_i - M)}{c * MAD} where :math:`c` is the tuning constant and :math:`MAD` is the `median absolute deviation `_. The biweight location tuning constant ``c`` is typically 6.0 (the default). If :math:`MAD` is zero, then the median will be returned. Parameters ---------- data : array-like Input array or object that can be converted to an array. ``data`` can be a `~numpy.ma.MaskedArray`. c : float, optional Tuning constant for the biweight estimator (default = 6.0). M : float or array-like, optional Initial guess for the location. If ``M`` is a scalar value, then its value will be used for the entire array (or along each ``axis``, if specified). If ``M`` is an array, then its must be an array containing the initial location estimate along each ``axis`` of the input array. If `None` (default), then the median of the input array will be used (or along each ``axis``, if specified). axis : None, int, or tuple of int, optional The axis or axes along which the biweight locations are computed. If `None` (default), then the biweight location of the flattened input array will be computed. ignore_nan : bool, optional Whether to ignore NaN values in the input ``data``. Returns ------- biweight_location : float or `~numpy.ndarray` The biweight location of the input data. If ``axis`` is `None` then a scalar will be returned, otherwise a `~numpy.ndarray` will be returned. See Also -------- biweight_scale, biweight_midvariance, biweight_midcovariance References ---------- .. [1] Beers, Flynn, and Gebhardt (1990; AJ 100, 32) (https://ui.adsabs.harvard.edu/abs/1990AJ....100...32B) .. [2] https://www.itl.nist.gov/div898/software/dataplot/refman2/auxillar/biwloc.htm Examples -------- Generate random variates from a Gaussian distribution and return the biweight location of the distribution: >>> import numpy as np >>> from astropy.stats import biweight_location >>> rand = np.random.default_rng(12345) >>> biloc = biweight_location(rand.standard_normal(1000)) >>> print(biloc) # doctest: +FLOAT_CMP 0.01535330525461019 )rT)copyN)axis)rrgignore)divideinvalid)rrr)rr r r r masked_whereisnan asanyarrayastypefloat64rrerrstateabssqueezeisscalarwhere) rcMrrrrdmadumaskvalue where_funcrrrrs:O    "@c Cstt||||||dS)a Compute the biweight scale. The biweight scale is a robust statistic for determining the standard deviation of a distribution. It is the square root of the `biweight midvariance `_. It is given by: .. math:: \zeta_{biscl} = \sqrt{n} \ \frac{\sqrt{\sum_{|u_i| < 1} \ (x_i - M)^2 (1 - u_i^2)^4}} {|(\sum_{|u_i| < 1} \ (1 - u_i^2) (1 - 5u_i^2))|} where :math:`x` is the input data, :math:`M` is the sample median (or the input location) and :math:`u_i` is given by: .. math:: u_{i} = \frac{(x_i - M)}{c * MAD} where :math:`c` is the tuning constant and :math:`MAD` is the `median absolute deviation `_. The biweight midvariance tuning constant ``c`` is typically 9.0 (the default). If :math:`MAD` is zero, then zero will be returned. For the standard definition of biweight scale, :math:`n` is the total number of points in the array (or along the input ``axis``, if specified). That definition is used if ``modify_sample_size`` is `False`, which is the default. However, if ``modify_sample_size = True``, then :math:`n` is the number of points for which :math:`|u_i| < 1` (i.e. the total number of non-rejected values), i.e. .. math:: n = \sum_{|u_i| < 1} \ 1 which results in a value closer to the true standard deviation for small sample sizes or for a large number of rejected values. Parameters ---------- data : array-like Input array or object that can be converted to an array. ``data`` can be a `~numpy.ma.MaskedArray`. c : float, optional Tuning constant for the biweight estimator (default = 9.0). M : float or array-like, optional The location estimate. If ``M`` is a scalar value, then its value will be used for the entire array (or along each ``axis``, if specified). If ``M`` is an array, then its must be an array containing the location estimate along each ``axis`` of the input array. If `None` (default), then the median of the input array will be used (or along each ``axis``, if specified). axis : None, int, or tuple of int, optional The axis or axes along which the biweight scales are computed. If `None` (default), then the biweight scale of the flattened input array will be computed. modify_sample_size : bool, optional If `False` (default), then the sample size used is the total number of elements in the array (or along the input ``axis``, if specified), which follows the standard definition of biweight scale. If `True`, then the sample size is reduced to correct for any rejected values (i.e. the sample size used includes only the non-rejected values), which results in a value closer to the true standard deviation for small sample sizes or for a large number of rejected values. ignore_nan : bool, optional Whether to ignore NaN values in the input ``data``. Returns ------- biweight_scale : float or `~numpy.ndarray` The biweight scale of the input data. If ``axis`` is `None` then a scalar will be returned, otherwise a `~numpy.ndarray` will be returned. See Also -------- biweight_midvariance, biweight_midcovariance, biweight_location, astropy.stats.mad_std, astropy.stats.median_absolute_deviation References ---------- .. [1] Beers, Flynn, and Gebhardt (1990; AJ 100, 32) (https://ui.adsabs.harvard.edu/abs/1990AJ....100...32B) .. [2] https://www.itl.nist.gov/div898/software/dataplot/refman2/auxillar/biwscale.htm Examples -------- Generate random variates from a Gaussian distribution and return the biweight scale of the distribution: >>> import numpy as np >>> from astropy.stats import biweight_scale >>> rand = np.random.default_rng(12345) >>> biscl = biweight_scale(rand.standard_normal(1000)) >>> print(biscl) # doctest: +FLOAT_CMP 1.0239311812635818 )r(r)rmodify_sample_sizer)r sqrtr)rr(r)rr1rrrrrsl c CsPt||d\}}t|tjjr:|r:tjjt||dd}t|tj }|dkr`|||d}|dk rtt ||d}||}t |||d} |dkr| dkrdSt| rtj Sn t | |d} tj ddd ||| } WdQRXtj dd t| d k} WdQRXt| tjjr$| jd d } | d} |r@|| |d} nHt|j} t|tjjrfd| |j<|rzd| t|<tj| |d} ||d| d}d|| <|||d}d| dd| }d|| <ttj||dd}tj ddd P| ||}t|r|Stj}t|tjjr2tjj}|| dkd|SQRXdS)a Compute the biweight midvariance. The biweight midvariance is a robust statistic for determining the variance of a distribution. Its square root is a robust estimator of scale (i.e. standard deviation). It is given by: .. math:: \zeta_{bivar} = n \ \frac{\sum_{|u_i| < 1} \ (x_i - M)^2 (1 - u_i^2)^4} {(\sum_{|u_i| < 1} \ (1 - u_i^2) (1 - 5u_i^2))^2} where :math:`x` is the input data, :math:`M` is the sample median (or the input location) and :math:`u_i` is given by: .. math:: u_{i} = \frac{(x_i - M)}{c * MAD} where :math:`c` is the tuning constant and :math:`MAD` is the `median absolute deviation `_. The biweight midvariance tuning constant ``c`` is typically 9.0 (the default). If :math:`MAD` is zero, then zero will be returned. For the standard definition of `biweight midvariance `_, :math:`n` is the total number of points in the array (or along the input ``axis``, if specified). That definition is used if ``modify_sample_size`` is `False`, which is the default. However, if ``modify_sample_size = True``, then :math:`n` is the number of points for which :math:`|u_i| < 1` (i.e. the total number of non-rejected values), i.e. .. math:: n = \sum_{|u_i| < 1} \ 1 which results in a value closer to the true variance for small sample sizes or for a large number of rejected values. Parameters ---------- data : array-like Input array or object that can be converted to an array. ``data`` can be a `~numpy.ma.MaskedArray`. c : float, optional Tuning constant for the biweight estimator (default = 9.0). M : float or array-like, optional The location estimate. If ``M`` is a scalar value, then its value will be used for the entire array (or along each ``axis``, if specified). If ``M`` is an array, then its must be an array containing the location estimate along each ``axis`` of the input array. If `None` (default), then the median of the input array will be used (or along each ``axis``, if specified). axis : None, int, or tuple of int, optional The axis or axes along which the biweight midvariances are computed. If `None` (default), then the biweight midvariance of the flattened input array will be computed. modify_sample_size : bool, optional If `False` (default), then the sample size used is the total number of elements in the array (or along the input ``axis``, if specified), which follows the standard definition of biweight midvariance. If `True`, then the sample size is reduced to correct for any rejected values (i.e. the sample size used includes only the non-rejected values), which results in a value closer to the true variance for small sample sizes or for a large number of rejected values. ignore_nan : bool, optional Whether to ignore NaN values in the input ``data``. Returns ------- biweight_midvariance : float or `~numpy.ndarray` The biweight midvariance of the input data. If ``axis`` is `None` then a scalar will be returned, otherwise a `~numpy.ndarray` will be returned. See Also -------- biweight_midcovariance, biweight_midcorrelation, astropy.stats.mad_std, astropy.stats.median_absolute_deviation References ---------- .. [1] https://en.wikipedia.org/wiki/Robust_measures_of_scale#The_biweight_midvariance .. [2] Beers, Flynn, and Gebhardt (1990; AJ 100, 32) (https://ui.adsabs.harvard.edu/abs/1990AJ....100...32B) Examples -------- Generate random variates from a Gaussian distribution and return the biweight midvariance of the distribution: >>> import numpy as np >>> from astropy.stats import biweight_midvariance >>> rand = np.random.default_rng(12345) >>> bivar = biweight_midvariance(rand.standard_normal(1000)) >>> print(bivar) # doctest: +FLOAT_CMP 1.0484350639638342 )rT)rN)r)rrgr)rr)rrF)Z fill_valuerrg?g@)rr r r r rrr r!r"rrnanr#r$ZfilledZonesshaper-rr&r'r%)rr(r)rr1rrrr*r+r,r-nZ include_maskf1f2r.r/rrrrsZk            c Cst|tj}|jdkr.|tjddf}|jdkr@td|dkrVtj|dd}t|}|jdkrrtd|j|j}t |dd}tj ddd|j||j}WdQRXtj dd t |dk}WdQRX|d}|r|t }t ||} n |d j} d |} d d |} d | |<tj ddd|| d} | | jddddtjf} t| | j}t| | j}| ||}t|d kd }d ||ddf<d |dd|f<|SQRXdS)a5 Compute the biweight midcovariance between pairs of multiple variables. The biweight midcovariance is a robust and resistant estimator of the covariance between two variables. This function computes the biweight midcovariance between all pairs of the input variables (rows) in the input data. The output array will have a shape of (N_variables, N_variables). The diagonal elements will be the biweight midvariances of each input variable (see :func:`biweight_midvariance`). The off-diagonal elements will be the biweight midcovariances between each pair of input variables. For example, if the input array ``data`` contains three variables (rows) ``x``, ``y``, and ``z``, the output `~numpy.ndarray` midcovariance matrix will be: .. math:: \begin{pmatrix} \zeta_{xx} & \zeta_{xy} & \zeta_{xz} \\ \zeta_{yx} & \zeta_{yy} & \zeta_{yz} \\ \zeta_{zx} & \zeta_{zy} & \zeta_{zz} \end{pmatrix} where :math:`\zeta_{xx}`, :math:`\zeta_{yy}`, and :math:`\zeta_{zz}` are the biweight midvariances of each variable. The biweight midcovariance between :math:`x` and :math:`y` is :math:`\zeta_{xy}` (:math:`= \zeta_{yx}`). The biweight midcovariance between :math:`x` and :math:`z` is :math:`\zeta_{xz}` (:math:`= \zeta_{zx}`). The biweight midcovariance between :math:`y` and :math:`z` is :math:`\zeta_{yz}` (:math:`= \zeta_{zy}`). The biweight midcovariance between two variables :math:`x` and :math:`y` is given by: .. math:: \zeta_{xy} = n_{xy} \ \frac{\sum_{|u_i| < 1, \ |v_i| < 1} \ (x_i - M_x) (1 - u_i^2)^2 (y_i - M_y) (1 - v_i^2)^2} {(\sum_{|u_i| < 1} \ (1 - u_i^2) (1 - 5u_i^2)) (\sum_{|v_i| < 1} \ (1 - v_i^2) (1 - 5v_i^2))} where :math:`M_x` and :math:`M_y` are the medians (or the input locations) of the two variables and :math:`u_i` and :math:`v_i` are given by: .. math:: u_{i} = \frac{(x_i - M_x)}{c * MAD_x} v_{i} = \frac{(y_i - M_y)}{c * MAD_y} where :math:`c` is the biweight tuning constant and :math:`MAD_x` and :math:`MAD_y` are the `median absolute deviation `_ of the :math:`x` and :math:`y` variables. The biweight midvariance tuning constant ``c`` is typically 9.0 (the default). If :math:`MAD_x` or :math:`MAD_y` are zero, then zero will be returned for that element. For the standard definition of biweight midcovariance, :math:`n_{xy}` is the total number of observations of each variable. That definition is used if ``modify_sample_size`` is `False`, which is the default. However, if ``modify_sample_size = True``, then :math:`n_{xy}` is the number of observations for which :math:`|u_i| < 1` and/or :math:`|v_i| < 1`, i.e. .. math:: n_{xx} = \sum_{|u_i| < 1} \ 1 .. math:: n_{xy} = n_{yx} = \sum_{|u_i| < 1, \ |v_i| < 1} \ 1 .. math:: n_{yy} = \sum_{|v_i| < 1} \ 1 which results in a value closer to the true variance for small sample sizes or for a large number of rejected values. Parameters ---------- data : 2D or 1D array-like Input data either as a 2D or 1D array. For a 2D array, it should have a shape (N_variables, N_observations). A 1D array may be input for observations of a single variable, in which case the biweight midvariance will be calculated (no covariance). Each row of ``data`` represents a variable, and each column a single observation of all those variables (same as the `numpy.cov` convention). c : float, optional Tuning constant for the biweight estimator (default = 9.0). M : float or 1D array-like, optional The location estimate of each variable, either as a scalar or array. If ``M`` is an array, then its must be a 1D array containing the location estimate of each row (i.e. ``a.ndim`` elements). If ``M`` is a scalar value, then its value will be used for each variable (row). If `None` (default), then the median of each variable (row) will be used. modify_sample_size : bool, optional If `False` (default), then the sample size used is the total number of observations of each variable, which follows the standard definition of biweight midcovariance. If `True`, then the sample size is reduced to correct for any rejected values (see formula above), which results in a value closer to the true covariance for small sample sizes or for a large number of rejected values. Returns ------- biweight_midcovariance : ndarray A 2D array representing the biweight midcovariances between each pair of the variables (rows) in the input array. The output array will have a shape of (N_variables, N_variables). The diagonal elements will be the biweight midvariances of each input variable. The off-diagonal elements will be the biweight midcovariances between each pair of input variables. See Also -------- biweight_midvariance, biweight_midcorrelation, biweight_scale, biweight_location References ---------- .. [1] https://www.itl.nist.gov/div898/software/dataplot/refman2/auxillar/biwmidc.htm Examples -------- Compute the biweight midcovariance between two random variables: >>> import numpy as np >>> from astropy.stats import biweight_midcovariance >>> # Generate two random variables x and y >>> rng = np.random.default_rng(1) >>> x = rng.normal(0, 1, 200) >>> y = rng.normal(0, 3, 200) >>> # Introduce an obvious outlier >>> x[0] = 30.0 >>> # Calculate the biweight midcovariances between x and y >>> bicov = biweight_midcovariance([x, y]) >>> print(bicov) # doctest: +FLOAT_CMP [[0.83435568 0.02379316] [0.02379316 7.15665769]] >>> # Print standard deviation estimates >>> print(np.sqrt(bicov.diagonal())) # doctest: +FLOAT_CMP [0.91343072 2.67519302] rNrz!The input array must be 2D or 1D.)rzM must be a scalar or 1D array.r)rr)rrg?g@g)r r r!r"ndimZnewaxis ValueErrorrTrr#r$floatinnersizerdotr')rr(r)r1r*r+r,r-Zmaskfr6Zusub1Zusub5 numerator denominatorZnumerator_matrixZdenominator_matrixr.idxrrrrsF             cCs~t|}t|}|jdkr&td|jdkr8td|j|jkrLtdt||g|||d}|dt|d|dS) a Compute the biweight midcorrelation between two variables. The `biweight midcorrelation `_ is a measure of similarity between samples. It is given by: .. math:: r_{bicorr} = \frac{\zeta_{xy}}{\sqrt{\zeta_{xx} \ \zeta_{yy}}} where :math:`\zeta_{xx}` is the biweight midvariance of :math:`x`, :math:`\zeta_{yy}` is the biweight midvariance of :math:`y`, and :math:`\zeta_{xy}` is the biweight midcovariance of :math:`x` and :math:`y`. Parameters ---------- x, y : 1D array-like Input arrays for the two variables. ``x`` and ``y`` must be 1D arrays and have the same number of elements. c : float, optional Tuning constant for the biweight estimator (default = 9.0). See `biweight_midcovariance` for more details. M : float or array-like, optional The location estimate. If ``M`` is a scalar value, then its value will be used for the entire array (or along each ``axis``, if specified). If ``M`` is an array, then its must be an array containing the location estimate along each ``axis`` of the input array. If `None` (default), then the median of the input array will be used (or along each ``axis``, if specified). See `biweight_midcovariance` for more details. modify_sample_size : bool, optional If `False` (default), then the sample size used is the total number of elements in the array (or along the input ``axis``, if specified), which follows the standard definition of biweight midcovariance. If `True`, then the sample size is reduced to correct for any rejected values (i.e. the sample size used includes only the non-rejected values), which results in a value closer to the true midcovariance for small sample sizes or for a large number of rejected values. See `biweight_midcovariance` for more details. Returns ------- biweight_midcorrelation : float The biweight midcorrelation between ``x`` and ``y``. See Also -------- biweight_scale, biweight_midvariance, biweight_midcovariance, biweight_location References ---------- .. [1] https://en.wikipedia.org/wiki/Biweight_midcorrelation Examples -------- Calculate the biweight midcorrelation between two variables: >>> import numpy as np >>> from astropy.stats import biweight_midcorrelation >>> rng = np.random.default_rng(12345) >>> x = rng.normal(0, 1, 200) >>> y = rng.normal(0, 3, 200) >>> # Introduce an obvious outlier >>> x[0] = 30.0 >>> bicorr = biweight_midcorrelation(x, y) >>> print(bicorr) # doctest: +FLOAT_CMP -0.09203238319481295 rzx must be a 1D array.zy must be a 1D array.z!x and y must have the same shape.)r(r)r1)rr)rr)rr)r r r9r:r5rr2)xyr(r)r1Zbicorrrrrr sI      )F)rNN)r0NNF)r0NNF)r0NF)r0NF) __doc__numpyr funcsrr__all__rrrrrr rrrrs q1 U