Python scipy.stats.circmean() Examples

def test_circmean_axis(self):
        x = np.array([[355, 5, 2, 359, 10, 350],
                      [351, 7, 4, 352, 9, 349],
                      [357, 9, 8, 358, 4, 356]])
        M1 = stats.circmean(x, high=360)
        M2 = stats.circmean(x.ravel(), high=360)
        assert_allclose(M1, M2, rtol=1e-14)

        M1 = stats.circmean(x, high=360, axis=1)
        M2 = [stats.circmean(x[i], high=360) for i in range(x.shape[0])]
        assert_allclose(M1, M2, rtol=1e-14)

        M1 = stats.circmean(x, high=360, axis=0)
        M2 = [stats.circmean(x[:, i], high=360) for i in range(x.shape[1])]
        assert_allclose(M1, M2, rtol=1e-14) 

def circ_corr(alpha1, alpha2):
    '''Helper function to compute the circular correlation.'''
    alpha1_bar = stats.circmean(alpha1)
    alpha2_bar = stats.circmean(alpha2)
    num = np.sum(np.sin(alpha1 - alpha1_bar) * np.sin(alpha2 - alpha2_bar))
    den = np.sqrt(np.sum(np.sin(alpha1 - alpha1_bar) ** 2) *
                  np.sum(np.sin(alpha2 - alpha2_bar) ** 2))
    rho = num / den
    return rho 

def _get_dxy(self, p, meta, m):
        """ Returns avg distance between cell edge and dot """
        #return 0,0
        #return di/2, dj/2
        #cropx, cropy, gridx, gridy, cellx, celly = meta
        x = sum(meta[slice(0,len(meta),2)])
        y = sum(meta[slice(1,len(meta),2)])
        dots = self.yd.dots.T/self.yd.imgDpi
        coords = dots[:,(dots[0] >= x) * (dots[0] < x+p.d_i*m.shape[0]) *
            (dots[1] >= y) * (dots[1] < y+p.d_j*m.shape[1])]
        dx = np.average((coords[0]-x)%p.d_i)
        dy = np.average((coords[1]-y)%p.d_j)
        #dx = circmean((coords[0]-x)%p.d_i, high=p.d_i)
        #dy = circmean((coords[1]-y)%p.d_j, high=p.d_j)
        return dx, dy 

def readTDM(self):
        """ 
        Extract TDM from scan 
        @returns: TDM, xOffset, yOffset
        """
        self._print("Extracting tracking dots... ")
        #cv2.imwrite("perspective.png",self.im)
        # get tracking dot matrices
        pp = PrintParser(self.im,ydxArgs=dict(
            inputDpi=self.dpi,interpolation=cv2.INTER_NEAREST,
            rotation=False))
        tdms = list(pp.getAllValidTdms())
        self._print("\t%d valid matrices"%len(tdms))
        
        # select tdm
        # tdms := [closest TDM at edge \forall edge \in Edges]
        # |tdms| = |Edges|
        #pp.tdm = random.choice(tdms)
        width = pp.yd.im.shape[1]*pp.yd.imgDpi
        height = pp.yd.im.shape[0]*pp.yd.imgDpi
        edges = [(0,0),(width,0),(0,height),(width,height)]
        tdms = [
            tdms[np.argmin(
                [abs(e1-tdm.atX)+abs(e2-tdm.atY) for tdm in tdms]
            )] for e1, e2 in edges]
        
        self._print("\tTracking dots pattern found:")
        self._print("\tx=%f, y=%f, trans=%s"%(pp.tdm.atX,pp.tdm.atY,pp.tdm.trans))

        xOffsets = [tdm.xoffset for tdm in tdms]
        yOffsets = [tdm.yoffset for tdm in tdms]
        xOffset = circmean(xOffsets, high=pp.tdm.hps_prototype)
        yOffset = circmean(yOffsets, high=pp.tdm.vps_prototype)
        self._print("TDM offset: %f, %f"%(xOffset,yOffset))
        
        return pp.tdm, xOffset, yOffset 

def radial_average(L):
    # r = round(math.degrees(cmath.phase(sum(cmath.rect(1, math.radians(float(d))) for d in L)/len(L))),2)
    scipy = circmean(L, -180, 180)
    return scipy 

def _hdi(ary, hdi_prob, circular, skipna):
    """Compute hpi over the flattened array."""
    ary = ary.flatten()
    if skipna:
        nans = np.isnan(ary)
        if not nans.all():
            ary = ary[~nans]
    n = len(ary)

    if circular:
        mean = st.circmean(ary, high=np.pi, low=-np.pi)
        ary = ary - mean
        ary = np.arctan2(np.sin(ary), np.cos(ary))

    ary = np.sort(ary)
    interval_idx_inc = int(np.floor(hdi_prob * n))
    n_intervals = n - interval_idx_inc
    interval_width = ary[interval_idx_inc:] - ary[:n_intervals]

    if len(interval_width) == 0:
        raise ValueError("Too few elements for interval calculation. ")

    min_idx = np.argmin(interval_width)
    hdi_min = ary[min_idx]
    hdi_max = ary[min_idx + interval_idx_inc]

    if circular:
        hdi_min = hdi_min + mean
        hdi_max = hdi_max + mean
        hdi_min = np.arctan2(np.sin(hdi_min), np.cos(hdi_min))
        hdi_max = np.arctan2(np.sin(hdi_max), np.cos(hdi_max))

    hdi_interval = np.array([hdi_min, hdi_max])

    return hdi_interval 

def test_can_predict_from_data():
    Invt_model = InvertedEncoding()
    Invt_model.fit(X, y)
    m_reconstruct = []
    for j in np.arange(dim):
        preds = Invt_model.predict(X2[n_*j:n_*(j+1), :])
        tmp = circmean(np.deg2rad(preds))
        m_reconstruct.append(np.rad2deg(tmp))
    logger.info('Reconstructed angles: ' + str(m_reconstruct))


# Show that prediction is invalid when input data is wrong size 

def test_circfuncs_unit8(self):
        # regression test for gh-7255: overflow when working with
        # numpy uint8 data type
        x = np.array([150, 10], dtype='uint8')
        assert_equal(stats.circmean(x, high=180), 170.0)
        assert_allclose(stats.circvar(x, high=180), 437.45871686, rtol=1e-7)
        assert_allclose(stats.circstd(x, high=180), 20.91551378, rtol=1e-7) 

def test_circmean_range(self):
        # regression test for gh-6420: circmean(..., high, low) must be
        # between `high` and `low`
        m = stats.circmean(np.arange(0, 2, 0.1), np.pi, -np.pi)
        assert_(m < np.pi)
        assert_(m > -np.pi) 

def test_empty(self):
        assert_(np.isnan(stats.circmean([])))
        assert_(np.isnan(stats.circstd([])))
        assert_(np.isnan(stats.circvar([]))) 

def test_circfuncs_array_like(self):
        x = [355, 5, 2, 359, 10, 350]
        assert_allclose(stats.circmean(x, high=360), 0.167690146, rtol=1e-7)
        assert_allclose(stats.circvar(x, high=360), 42.51955609, rtol=1e-7)
        assert_allclose(stats.circstd(x, high=360), 6.520702116, rtol=1e-7) 

def test_circmean(self):
        """ Test custom circular mean."""

        ref_mean = scistats.circmean(self.test_angles, **self.circ_kwargs)
        test_mean = pystats.nan_circmean(self.test_angles, **self.circ_kwargs)

        assert ref_mean == test_mean 

def test_circfuncs_small(self):
        x = np.array([20, 21, 22, 18, 19, 20.5, 19.2])
        M1 = x.mean()
        M2 = stats.circmean(x, high=360)
        assert_allclose(M2, M1, rtol=1e-5)

        V1 = x.var()
        V2 = stats.circvar(x, high=360)
        assert_allclose(V2, V1, rtol=1e-4)

        S1 = x.std()
        S2 = stats.circstd(x, high=360)
        assert_allclose(S2, S1, rtol=1e-4) 

def test_circfuncs(self):
        x = np.array([355, 5, 2, 359, 10, 350])
        M = stats.circmean(x, high=360)
        Mval = 0.167690146
        assert_allclose(M, Mval, rtol=1e-7)

        V = stats.circvar(x, high=360)
        Vval = 42.51955609
        assert_allclose(V, Vval, rtol=1e-7)

        S = stats.circstd(x, high=360)
        Sval = 6.520702116
        assert_allclose(S, Sval, rtol=1e-7) 

def test_empty(self):
        assert_(np.isnan(stats.circmean([])))
        assert_(np.isnan(stats.circstd([])))
        assert_(np.isnan(stats.circvar([]))) 

def test_circfuncs_array_like(self):
        x = [355,5,2,359,10,350]
        assert_allclose(stats.circmean(x, high=360), 0.167690146, rtol=1e-7)
        assert_allclose(stats.circvar(x, high=360), 42.51955609, rtol=1e-7)
        assert_allclose(stats.circstd(x, high=360), 6.520702116, rtol=1e-7) 

def test_circmean_axis(self):
        x = np.array([[355,5,2,359,10,350],
                      [351,7,4,352,9,349],
                      [357,9,8,358,4,356]])
        M1 = stats.circmean(x, high=360)
        M2 = stats.circmean(x.ravel(), high=360)
        assert_allclose(M1, M2, rtol=1e-14)

        M1 = stats.circmean(x, high=360, axis=1)
        M2 = [stats.circmean(x[i], high=360) for i in range(x.shape[0])]
        assert_allclose(M1, M2, rtol=1e-14)

        M1 = stats.circmean(x, high=360, axis=0)
        M2 = [stats.circmean(x[:,i], high=360) for i in range(x.shape[1])]
        assert_allclose(M1, M2, rtol=1e-14) 

def test_circfuncs_small(self):
        x = np.array([20,21,22,18,19,20.5,19.2])
        M1 = x.mean()
        M2 = stats.circmean(x, high=360)
        assert_allclose(M2, M1, rtol=1e-5)

        V1 = x.var()
        V2 = stats.circvar(x, high=360)
        assert_allclose(V2, V1, rtol=1e-4)

        S1 = x.std()
        S2 = stats.circstd(x, high=360)
        assert_allclose(S2, S1, rtol=1e-4) 

def test_circfuncs(self):
        x = np.array([355,5,2,359,10,350])
        M = stats.circmean(x, high=360)
        Mval = 0.167690146
        assert_allclose(M, Mval, rtol=1e-7)

        V = stats.circvar(x, high=360)
        Vval = 42.51955609
        assert_allclose(V, Vval, rtol=1e-7)

        S = stats.circstd(x, high=360)
        Sval = 6.520702116
        assert_allclose(S, Sval, rtol=1e-7) 

def test_deprecation_warning_circmean(self):
        """Test if circmean in stats is deprecated"""

        with warnings.catch_warnings(record=True) as war:
            try:
                pystats.nan_circmean(None)
            except TypeError:
                # Setting input to None should produce a TypeError after
                # warning is generated
                pass

        assert len(war) >= 1
        assert war[0].category == DeprecationWarning 

def test_circmean_nan(self):
        """ Test custom circular mean with NaN."""

        ref_mean = scistats.circmean(self.test_angles, **self.circ_kwargs)
        ref_nan = scistats.circmean(self.test_nan, **self.circ_kwargs)
        test_nan = pystats.nan_circmean(self.test_nan, **self.circ_kwargs)

        assert np.isnan(ref_nan)
        assert ref_mean == test_nan 

def _mc_error(ary, batches=5, circular=False):
    """Calculate the simulation standard error, accounting for non-independent samples.

    The trace is divided into batches, and the standard deviation of the batch
    means is calculated.

    Parameters
    ----------
    ary : Numpy array
        An array containing MCMC samples
    batches : integer
        Number of batches
    circular : bool
        Whether to compute the error taking into account `ary` is a circular variable
        (in the range [-np.pi, np.pi]) or not. Defaults to False (i.e non-circular variables).

    Returns
    -------
    mc_error : float
        Simulation standard error
    """
    _numba_flag = Numba.numba_flag
    if ary.ndim > 1:

        dims = np.shape(ary)
        trace = np.transpose([t.ravel() for t in ary])

        return np.reshape([_mc_error(t, batches) for t in trace], dims[1:])

    else:
        if _not_valid(ary, check_shape=False):
            return np.nan
        if batches == 1:
            if circular:
                if _numba_flag:
                    std = _circular_standard_deviation(ary, high=np.pi, low=-np.pi)
                else:
                    std = stats.circstd(ary, high=np.pi, low=-np.pi)
            else:
                if _numba_flag:
                    std = np.float(_sqrt(svar(ary), np.zeros(1)))
                else:
                    std = np.std(ary)
            return std / np.sqrt(len(ary))

        batched_traces = np.resize(ary, (batches, int(len(ary) / batches)))

        if circular:
            means = stats.circmean(batched_traces, high=np.pi, low=-np.pi, axis=1)
            if _numba_flag:
                std = _circular_standard_deviation(means, high=np.pi, low=-np.pi)
            else:
                std = stats.circstd(means, high=np.pi, low=-np.pi)
        else:
            means = np.mean(batched_traces, 1)
            if _numba_flag:
                std = _sqrt(svar(means), np.zeros(1))
            else:
                std = np.std(means)

        return std / np.sqrt(batches) 

def test_helper_angles(self):
        """Test helper circular functions."""
        # Check angles
        _checkangles(a1)
        _checkangles(a2, axis=None)
        with pytest.raises(ValueError):
            _checkangles(a3)
        with pytest.raises(ValueError):
            _checkangles(a3, axis=None)
        # Convert angles
        np.testing.assert_array_almost_equal(a1, convert_angles(a1, low=-np.pi,
                                                                high=np.pi))
        np.testing.assert_array_almost_equal(a2, convert_angles(a2, low=0,
                                                                high=2 * np.pi,
                                                                positive=True))
        _checkangles(convert_angles(a2, low=0, high=2 * np.pi))
        _checkangles(convert_angles(a3, low=0, high=360))
        _checkangles(convert_angles(a3_nan, low=0, high=360))
        _checkangles(convert_angles(a4, low=0, high=24))
        _checkangles(convert_angles(a5, low=0, high=1440))
        _checkangles(convert_angles(a5_nan, low=0, high=1440))
        _checkangles(convert_angles(a5, low=0, high=1440, positive=True))

        convert_angles(a1, low=-np.pi, high=np.pi)
        convert_angles(a2, low=0, high=2 * np.pi)
        convert_angles(a3)
        convert_angles(a4, low=0, high=24)
        convert_angles(a5, low=0, high=1440)
        assert convert_angles(a5, low=0, high=1440, positive=True).min() >= 0
        assert convert_angles(a5, low=0, high=1440).min() <= 0

        # Compare with scipy.stats.circmean
        def assert_circmean(x, low, high, axis=-1):
            m1 = convert_angles(circmean(x, low=low, high=high, axis=axis,
                                         nan_policy='omit'), low, high)
            m2 = circ_mean(convert_angles(x, low, high), axis=axis)
            assert (np.round(m1, 4) == np.round(m2, 4)).all()

        assert_circmean(a1, low=-np.pi, high=np.pi)
        assert_circmean(a2, low=0, high=2 * np.pi)
        assert_circmean(a3, low=0, high=360)
        assert_circmean(a3_nan, low=0, high=360)
        assert_circmean(a4, low=0, high=24)
        assert_circmean(a5, low=0, high=1440, axis=1)
        assert_circmean(a5, low=0, high=1440, axis=0)
        assert_circmean(a5, low=0, high=1440, axis=None)
        assert_circmean(a5_nan, low=0, high=1440, axis=-1)
        assert_circmean(a5_nan, low=0, high=1440, axis=0)
        assert_circmean(a5_nan, low=0, high=1440, axis=None)