diff --git a/dpnp/backend_statistics.pyx b/dpnp/backend_statistics.pyx
index 264bad74d94e..ce5526f16cf3 100644
--- a/dpnp/backend_statistics.pyx
+++ b/dpnp/backend_statistics.pyx
@@ -38,10 +38,44 @@ from dpnp.dpnp_utils cimport checker_throw_type_error, normalize_axis
 
 
 __all__ += [
+    "dpnp_cov",
     "dpnp_mean"
 ]
 
 
+cpdef dparray dpnp_cov(dparray array1):
+    # behaviour of original numpy
+    if array1.ndim > 2:
+        raise ValueError("array has more than 2 dimensions")
+
+    if array1.ndim < 2:
+        raise NotImplementedError
+
+    # numpy provide result as float64 for any input type
+    cdef dparray mean = dparray(array1.shape[0], dtype=numpy.float64)
+    cdef dparray X = dparray(array1.shape, dtype=numpy.float64)
+
+    # mean(array1, axis=1) #################################
+    # dpmp.mean throws: 'dpnp.dparray.dparray' object is not callable
+    for i in range(array1.shape[0]):
+        sum = 0.0
+        for j in range(array1.shape[1]):
+            sum += array1[i, j]
+        mean[i] = sum / array1.shape[1]
+    ########################################################
+    #X = array1 - mean[:, None]
+    #X = array1 - mean[:, numpy.newaxis]
+    #X = array1 - mean.reshape((array1.shape[0], 1))
+    for i in range(array1.shape[0]):
+        for j in range(array1.shape[1]):
+            X[i, j] = array1[i, j] - mean[i]
+    ########################################################
+    Y = X.transpose()
+    res = dpnp_matmul(X, Y) / (array1.shape[1] - 1)
+
+    return res
+
+
 cpdef dparray dpnp_mean(dparray a, axis):
     cdef dparray_shape_type shape_a = a.shape
     cdef long size_a = a.size
diff --git a/dpnp/dpnp_iface_statistics.py b/dpnp/dpnp_iface_statistics.py
index 7887f3644752..01beafcc852e 100644
--- a/dpnp/dpnp_iface_statistics.py
+++ b/dpnp/dpnp_iface_statistics.py
@@ -47,10 +47,37 @@
 from dpnp.dpnp_utils import checker_throw_value_error, use_origin_backend
 
 __all__ = [
+    'cov',
     'mean'
 ]
 
 
+def cov(in_array1, y=None, rowvar=True, bias=False, ddof=None, fweights=None, aweights=None):
+    """
+    Estimate a covariance matrix, given data and weights.
+    """
+
+    is_dparray1 = isinstance(in_array1, dparray)
+
+    if (not use_origin_backend(in_array1) and is_dparray1):
+        if y is not None:
+            checker_throw_value_error("cov", "y", type(y), None)
+        if rowvar is not True:
+            checker_throw_value_error("cov", "rowvar", rowvar, True)
+        if bias is not False:
+            checker_throw_value_error("cov", "bias", bias, False)
+        if ddof is not None:
+            checker_throw_value_error("cov", "ddof", type(ddof), None)
+        if fweights is not None:
+            checker_throw_value_error("cov", "fweights", type(fweights), None)
+        if aweights is not None:
+            checker_throw_value_error("cov", "aweights", type(aweights), None)
+
+        return dpnp_cov(in_array1)
+
+    return numpy.cov(in_array1)
+
+
 def mean(a, axis=None):
     """
     Compute the arithmetic mean along the specified axis.
diff --git a/tests/third_party/cupy/statistics_tests/__init__.py b/tests/third_party/cupy/statics_tests/__init__.py
similarity index 100%
rename from tests/third_party/cupy/statistics_tests/__init__.py
rename to tests/third_party/cupy/statics_tests/__init__.py
diff --git a/tests/third_party/cupy/statics_tests/test_correlation.py b/tests/third_party/cupy/statics_tests/test_correlation.py
new file mode 100644
index 000000000000..f23d49114277
--- /dev/null
+++ b/tests/third_party/cupy/statics_tests/test_correlation.py
@@ -0,0 +1,166 @@
+import unittest
+
+import numpy
+import pytest
+
+import dpnp as cupy
+from tests.third_party.cupy import testing
+
+
+@testing.gpu
+class TestCorrcoef(unittest.TestCase):
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_allclose()
+    def test_corrcoef(self, xp, dtype):
+        a = testing.shaped_arange((2, 3), xp, dtype)
+        return xp.corrcoef(a)
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_allclose()
+    def test_corrcoef_diag_exception(self, xp, dtype):
+        a = testing.shaped_arange((1, 3), xp, dtype)
+        return xp.corrcoef(a)
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_allclose()
+    def test_corrcoef_y(self, xp, dtype):
+        a = testing.shaped_arange((2, 3), xp, dtype)
+        y = testing.shaped_arange((2, 3), xp, dtype)
+        return xp.corrcoef(a, y=y)
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_allclose()
+    def test_corrcoef_rowvar(self, xp, dtype):
+        a = testing.shaped_arange((2, 3), xp, dtype)
+        y = testing.shaped_arange((2, 3), xp, dtype)
+        return xp.corrcoef(a, y=y, rowvar=False)
+
+
+@testing.gpu
+class TestCov(unittest.TestCase):
+
+    def generate_input(self, a_shape, y_shape, xp, dtype):
+        a = testing.shaped_arange(a_shape, xp, dtype)
+        y = None
+        if y_shape is not None:
+            y = testing.shaped_arange(y_shape, xp, dtype)
+        return a, y
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_allclose()
+    def check(self, a_shape, y_shape=None, rowvar=True, bias=False,
+              ddof=None, xp=None, dtype=None):
+        a, y = self.generate_input(a_shape, y_shape, xp, dtype)
+        return xp.cov(a, y, rowvar, bias, ddof)
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_allclose()
+    def check_warns(self, a_shape, y_shape=None, rowvar=True, bias=False,
+                    ddof=None, xp=None, dtype=None):
+        with testing.assert_warns(RuntimeWarning):
+            a, y = self.generate_input(a_shape, y_shape, xp, dtype)
+            return xp.cov(a, y, rowvar, bias, ddof)
+
+    @testing.for_all_dtypes()
+    def check_raises(self, a_shape, y_shape=None, rowvar=True, bias=False,
+                     ddof=None, dtype=None):
+        for xp in (numpy, cupy):
+            a, y = self.generate_input(a_shape, y_shape, xp, dtype)
+            with pytest.raises(ValueError):
+                xp.cov(a, y, rowvar, bias, ddof)
+
+    def test_cov(self):
+        self.check((2, 3))
+        self.check((2,), (2,))
+        self.check((1, 3), (1, 3), rowvar=False)
+        self.check((2, 3), (2, 3), rowvar=False)
+        self.check((2, 3), bias=True)
+        self.check((2, 3), ddof=2)
+
+    def test_cov_warns(self):
+        self.check_warns((2, 3), ddof=3)
+        self.check_warns((2, 3), ddof=4)
+
+    def test_cov_raises(self):
+        self.check_raises((2, 3), ddof=1.2)
+        self.check_raises((3, 4, 2))
+        self.check_raises((2, 3), (3, 4, 2))
+
+    def test_cov_empty(self):
+        self.check((0, 1))
+
+
+@testing.gpu
+@testing.parameterize(*testing.product({
+    'mode': ['valid', 'same', 'full'],
+    'shape1': [(5,), (6,), (20,), (21,)],
+    'shape2': [(5,), (6,), (20,), (21,)],
+}))
+class TestCorrelateShapeCombination(unittest.TestCase):
+
+    @testing.for_all_dtypes(no_float16=True)
+    @testing.numpy_cupy_allclose(rtol=1e-4)
+    def test_correlate(self, xp, dtype):
+        a = testing.shaped_arange(self.shape1, xp, dtype)
+        b = testing.shaped_arange(self.shape2, xp, dtype)
+        return xp.correlate(a, b, mode=self.mode)
+
+
+@testing.gpu
+@testing.parameterize(*testing.product({
+    'mode': ['valid', 'full', 'same']
+}))
+class TestCorrelate(unittest.TestCase):
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_allclose(rtol=1e-5)
+    def test_correlate_non_contiguous(self, xp, dtype):
+        a = testing.shaped_arange((300,), xp, dtype)
+        b = testing.shaped_arange((100,), xp, dtype)
+        return xp.correlate(a[::200], b[10::70], mode=self.mode)
+
+    @testing.for_all_dtypes(no_float16=True)
+    @testing.numpy_cupy_allclose(rtol=1e-4)
+    def test_correlate_large_non_contiguous(self, xp, dtype):
+        a = testing.shaped_arange((10000,), xp, dtype)
+        b = testing.shaped_arange((1000,), xp, dtype)
+        return xp.correlate(a[200::], b[10::700], mode=self.mode)
+
+    @testing.for_all_dtypes_combination(names=['dtype1', 'dtype2'])
+    @testing.numpy_cupy_allclose(rtol=1e-2)
+    def test_correlate_diff_types(self, xp, dtype1, dtype2):
+        a = testing.shaped_random((200,), xp, dtype1)
+        b = testing.shaped_random((100,), xp, dtype2)
+        return xp.correlate(a, b, mode=self.mode)
+
+
+@testing.gpu
+@testing.parameterize(*testing.product({
+    'mode': ['valid', 'same', 'full']
+}))
+class TestCorrelateInvalid(unittest.TestCase):
+
+    @testing.with_requires('numpy>=1.18')
+    @testing.for_all_dtypes()
+    def test_correlate_empty(self, dtype):
+        for xp in (numpy, cupy):
+            a = xp.zeros((0,), dtype)
+            with pytest.raises(ValueError):
+                xp.correlate(a, a, mode=self.mode)
+
+    @testing.for_all_dtypes()
+    def test_correlate_ndim(self, dtype):
+        for xp in (numpy, cupy):
+            a = testing.shaped_arange((5, 10, 2), xp, dtype)
+            b = testing.shaped_arange((3, 4, 4), xp, dtype)
+            with pytest.raises(ValueError):
+                xp.correlate(a, b, mode=self.mode)
+
+    @testing.for_all_dtypes()
+    def test_correlate_zero_dim(self, dtype):
+        for xp in (numpy, cupy):
+            a = testing.shaped_arange((), xp, dtype)
+            b = testing.shaped_arange((1,), xp, dtype)
+            with pytest.raises(ValueError):
+                xp.correlate(a, b, mode=self.mode)
diff --git a/tests/third_party/cupy/statics_tests/test_histogram.py b/tests/third_party/cupy/statics_tests/test_histogram.py
new file mode 100644
index 000000000000..0560dc00256e
--- /dev/null
+++ b/tests/third_party/cupy/statics_tests/test_histogram.py
@@ -0,0 +1,485 @@
+import sys
+import unittest
+
+import numpy
+import pytest
+
+import dpnp as cupy
+from tests.third_party.cupy import testing
+# from cupy.core import _accelerator
+
+
+# Note that numpy.bincount does not support uint64 on 64-bit environment
+# as it casts an input array to intp.
+# And it does not support uint32, int64 and uint64 on 32-bit environment.
+_all_types = (
+    numpy.float16, numpy.float32, numpy.float64,
+    numpy.int8, numpy.int16, numpy.int32,
+    numpy.uint8, numpy.uint16,
+    numpy.bool_)
+_signed_types = (
+    numpy.int8, numpy.int16, numpy.int32,
+    numpy.bool_)
+
+if sys.maxsize > 2 ** 32:
+    _all_types = _all_types + (numpy.int64, numpy.uint32)
+    _signed_types = _signed_types + (numpy.int64,)
+
+
+def for_all_dtypes_bincount(name='dtype'):
+    return testing.for_dtypes(_all_types, name=name)
+
+
+def for_signed_dtypes_bincount(name='dtype'):
+    return testing.for_dtypes(_signed_types, name=name)
+
+
+def for_all_dtypes_combination_bincount(names):
+    return testing.helper.for_dtypes_combination(_all_types, names=names)
+
+
+@testing.gpu
+class TestHistogram(unittest.TestCase):
+
+    @testing.for_all_dtypes(no_bool=True, no_complex=True)
+    @testing.numpy_cupy_array_equal()
+    def test_histogram(self, xp, dtype):
+        x = testing.shaped_arange((10,), xp, dtype)
+        y, bin_edges = xp.histogram(x)
+        return y, bin_edges
+
+    @testing.for_all_dtypes(no_bool=True, no_complex=True)
+    @testing.numpy_cupy_array_equal()
+    def test_histogram_same_value(self, xp, dtype):
+        x = xp.zeros(10, dtype)
+        y, bin_edges = xp.histogram(x, 3)
+        return y, bin_edges
+
+    @testing.for_all_dtypes(no_bool=True, no_complex=True)
+    @testing.numpy_cupy_array_equal()
+    def test_histogram_density(self, xp, dtype):
+        x = testing.shaped_arange((10,), xp, dtype)
+        y, bin_edges = xp.histogram(x, density=True)
+        # check normalization
+        area = xp.sum(y * xp.diff(bin_edges))
+        testing.assert_allclose(area, 1)
+        return y, bin_edges
+
+    @testing.for_float_dtypes()
+    @testing.numpy_cupy_array_equal()
+    def test_histogram_range_lower_outliers(self, xp, dtype):
+        # Check that lower outliers are not tallied
+        a = xp.arange(10, dtype=dtype) + .5
+        h, b = xp.histogram(a, range=[0, 9])
+        assert int(h.sum()) == 9
+        return h, b
+
+    @testing.for_float_dtypes()
+    @testing.numpy_cupy_array_equal()
+    def test_histogram_range_upper_outliers(self, xp, dtype):
+        # Check that upper outliers are not tallied
+        a = xp.arange(10, dtype=dtype) + .5
+        h, b = xp.histogram(a, range=[1, 10])
+        assert int(h.sum()) == 9
+        return h, b
+
+    @testing.for_float_dtypes()
+    @testing.numpy_cupy_allclose()
+    def test_histogram_range_with_density(self, xp, dtype):
+        a = xp.arange(10, dtype=dtype) + .5
+        h, b = xp.histogram(a, range=[1, 9], density=True)
+        # check normalization
+        testing.assert_allclose(float((h * xp.diff(b)).sum()), 1)
+        return h
+
+    @testing.for_float_dtypes()
+    @testing.numpy_cupy_allclose()
+    def test_histogram_range_with_weights_and_density(self, xp, dtype):
+        a = xp.arange(10, dtype=dtype) + .5
+        w = xp.arange(10, dtype=dtype) + .5
+        h, b = xp.histogram(a, range=[1, 9], weights=w, density=True)
+        testing.assert_allclose(float((h * xp.diff(b)).sum()), 1)
+        return h
+
+    def test_histogram_invalid_range(self):
+        for xp in (numpy, cupy):
+            with pytest.raises(ValueError):
+                # range must be None or have two elements
+                xp.histogram(xp.arange(10), range=[1, 9, 15])
+
+    def test_histogram_invalid_range2(self):
+        for xp in (numpy, cupy):
+            with pytest.raises(TypeError):
+                xp.histogram(xp.arange(10), range=10)
+
+    @testing.for_all_dtypes(no_bool=True, no_complex=True)
+    def test_histogram_weights_mismatch(self, dtype):
+        for xp in (numpy, cupy):
+            a = xp.arange(10, dtype=dtype) + .5
+            w = xp.arange(11, dtype=dtype) + .5
+            with pytest.raises(ValueError):
+                xp.histogram(a, range=[1, 9], weights=w, density=True)
+
+    @testing.for_all_dtypes(no_bool=True, no_complex=True)
+    @testing.numpy_cupy_allclose()
+    def test_histogram_int_weights_dtype(self, xp, dtype):
+        # Check the type of the returned histogram
+        a = xp.arange(10, dtype=dtype)
+        h, b = xp.histogram(a, weights=xp.ones(10, int))
+        assert xp.issubdtype(h.dtype, xp.integer)
+        return h
+
+    @testing.for_all_dtypes(no_bool=True, no_complex=True)
+    @testing.numpy_cupy_allclose()
+    def test_histogram_float_weights_dtype(self, xp, dtype):
+        # Check the type of the returned histogram
+        a = xp.arange(10, dtype=dtype)
+        h, b = xp.histogram(a, weights=xp.ones(10, float))
+        assert xp.issubdtype(h.dtype, xp.floating)
+        return h
+
+    def test_histogram_weights_basic(self):
+        v = cupy.random.rand(100)
+        w = cupy.ones(100) * 5
+        a, b = cupy.histogram(v)
+        na, nb = cupy.histogram(v, density=True)
+        wa, wb = cupy.histogram(v, weights=w)
+        nwa, nwb = cupy.histogram(v, weights=w, density=True)
+        testing.assert_array_almost_equal(a * 5, wa)
+        testing.assert_array_almost_equal(na, nwa)
+
+    @testing.for_float_dtypes()
+    @testing.numpy_cupy_allclose()
+    def test_histogram_float_weights(self, xp, dtype):
+        # Check weights are properly applied.
+        v = xp.linspace(0, 10, 10, dtype=dtype)
+        w = xp.concatenate((xp.zeros(5, dtype=dtype), xp.ones(5, dtype=dtype)))
+        wa, wb = xp.histogram(v, bins=xp.arange(11), weights=w)
+        testing.assert_array_almost_equal(wa, w)
+        return wb
+
+    @testing.for_int_dtypes(no_bool=True)
+    @testing.numpy_cupy_array_equal(type_check=False)
+    def test_histogram_int_weights(self, xp, dtype):
+        # Check with integer weights
+        v = xp.asarray([1, 2, 2, 4], dtype=dtype)
+        w = xp.asarray([4, 3, 2, 1], dtype=dtype)
+        wa, wb = xp.histogram(v, bins=4, weights=w)
+        testing.assert_array_equal(wa, [4, 5, 0, 1])
+        return wa, wb
+
+    @testing.for_int_dtypes(no_bool=True)
+    @testing.numpy_cupy_allclose()
+    def test_histogram_int_weights_normalized(self, xp, dtype):
+        v = xp.asarray([1, 2, 2, 4], dtype=dtype)
+        w = xp.asarray([4, 3, 2, 1], dtype=dtype)
+        wa, wb = xp.histogram(v, bins=4, weights=w, density=True)
+        testing.assert_array_almost_equal(
+            wa, xp.asarray([4, 5, 0, 1]) / 10. / 3. * 4)
+        return wb
+
+    @testing.for_int_dtypes(no_bool=True)
+    @testing.numpy_cupy_array_equal()
+    def test_histogram_int_weights_nonuniform_bins(self, xp, dtype):
+        # Check weights with non-uniform bin widths
+        a, b = xp.histogram(
+            xp.arange(9, dtype=dtype),
+            xp.asarray([0, 1, 3, 6, 10], dtype=dtype),
+            weights=xp.asarray([2, 1, 1, 1, 1, 1, 1, 1, 1], dtype=dtype),
+            density=True)
+        testing.assert_array_almost_equal(a, [.2, .1, .1, .075])
+        return a, b
+
+    @testing.for_complex_dtypes()
+    @testing.numpy_cupy_array_equal(type_check=False)
+    def test_histogram_complex_weights(self, xp, dtype):
+        values = xp.asarray([1.3, 2.5, 2.3])
+        weights = xp.asarray([1, -1, 2]) + 1j * xp.asarray([2, 1, 2])
+        weights = weights.astype(dtype)
+        a, b = xp.histogram(
+            values, bins=2, weights=weights)
+        return a, b
+
+    @testing.for_complex_dtypes()
+    @testing.numpy_cupy_array_equal(type_check=False)
+    def test_histogram_complex_weights_uneven_bins(self, xp, dtype):
+        values = xp.asarray([1.3, 2.5, 2.3])
+        weights = xp.asarray([1, -1, 2]) + 1j * xp.asarray([2, 1, 2])
+        weights = weights.astype(dtype)
+        a, b = xp.histogram(
+            values, bins=xp.asarray([0, 2, 3]), weights=weights)
+        return a, b
+
+    @testing.for_all_dtypes(no_bool=True, no_complex=True)
+    @testing.numpy_cupy_array_equal()
+    def test_histogram_empty(self, xp, dtype):
+        x = xp.array([], dtype)
+        y, bin_edges = xp.histogram(x)
+        return y, bin_edges
+
+    @testing.for_all_dtypes(no_bool=True, no_complex=True)
+    @testing.numpy_cupy_array_equal()
+    def test_histogram_int_bins(self, xp, dtype):
+        x = testing.shaped_arange((10,), xp, dtype)
+        y, bin_edges = xp.histogram(x, 4)
+        return y, bin_edges
+
+    @testing.for_all_dtypes(no_bool=True, no_complex=True)
+    @testing.numpy_cupy_array_equal()
+    def test_histogram_array_bins(self, xp, dtype):
+        x = testing.shaped_arange((10,), xp, dtype)
+        bins = testing.shaped_arange((3,), xp, dtype)
+        y, bin_edges = xp.histogram(x, bins)
+        return y, bin_edges
+
+    @testing.for_all_dtypes(no_bool=True, no_complex=True)
+    @testing.numpy_cupy_array_equal()
+    def test_histogram_numpy_bins(self, xp, dtype):
+        x = testing.shaped_arange((10,), xp, dtype)
+        bins = testing.shaped_arange((3,), numpy, dtype)
+        y, bin_edges = xp.histogram(x, bins)
+        return y, bin_edges
+
+    @testing.for_all_dtypes(no_bool=True, no_complex=True)
+    @testing.numpy_cupy_array_equal()
+    def test_histogram_list_bins(self, xp, dtype):
+        x = testing.shaped_arange((10,), xp, dtype)
+        bins = list(testing.shaped_arange((3,), numpy, dtype))
+        y, bin_edges = xp.histogram(x, bins)
+        return y, bin_edges
+
+    # numpy 1.13.1 does not check this error correctly with unsigned int.
+    @testing.for_all_dtypes(no_bool=True, no_complex=True)
+    def test_histogram_bins_not_ordered(self, dtype):
+        for xp in (numpy, cupy):
+            x = testing.shaped_arange((10,), xp, dtype)
+            bins = xp.array([1, 3, 2], dtype)
+            with pytest.raises(ValueError):
+                xp.histogram(x, bins)
+
+    @for_all_dtypes_bincount()
+    @testing.numpy_cupy_allclose(accept_error=TypeError)
+    def test_bincount(self, xp, dtype):
+        x = testing.shaped_arange((3,), xp, dtype)
+        return xp.bincount(x)
+
+    @for_all_dtypes_bincount()
+    @testing.numpy_cupy_allclose(accept_error=TypeError)
+    def test_bincount_duplicated_value(self, xp, dtype):
+        x = xp.array([1, 2, 2, 1, 2, 4], dtype)
+        return xp.bincount(x)
+
+    @for_all_dtypes_combination_bincount(names=['x_type', 'w_type'])
+    @testing.numpy_cupy_allclose(accept_error=TypeError)
+    def test_bincount_with_weight(self, xp, x_type, w_type):
+        x = testing.shaped_arange((3,), xp, x_type)
+        w = testing.shaped_arange((3,), xp, w_type)
+        return xp.bincount(x, weights=w)
+
+    @for_all_dtypes_bincount()
+    @testing.numpy_cupy_allclose(accept_error=TypeError)
+    def test_bincount_with_minlength(self, xp, dtype):
+        x = testing.shaped_arange((3,), xp, dtype)
+        return xp.bincount(x, minlength=5)
+
+    @for_all_dtypes_combination_bincount(names=['x_type', 'w_type'])
+    def test_bincount_invalid_weight_length(self, x_type, w_type):
+        for xp in (numpy, cupy):
+            x = testing.shaped_arange((1,), xp, x_type)
+            w = testing.shaped_arange((2,), xp, w_type)
+            # TODO(imanishi): Split this test into a test for ValueError and
+            # a test for TypeError.
+            with pytest.raises((ValueError, TypeError)):
+                xp.bincount(x, weights=w)
+
+    @for_signed_dtypes_bincount()
+    def test_bincount_negative(self, dtype):
+        for xp in (numpy, cupy):
+            x = testing.shaped_arange((3,), xp, dtype) - 2
+            with pytest.raises(ValueError):
+                xp.bincount(x)
+
+    @for_all_dtypes_bincount()
+    def test_bincount_too_deep(self, dtype):
+        for xp in (numpy, cupy):
+            x = xp.array([[1]], dtype)
+            with pytest.raises(ValueError):
+                xp.bincount(x)
+
+    @for_all_dtypes_bincount()
+    def test_bincount_too_small(self, dtype):
+        for xp in (numpy, cupy):
+            x = xp.zeros((), dtype)
+            with pytest.raises(ValueError):
+                xp.bincount(x)
+
+    @for_all_dtypes_bincount()
+    @testing.numpy_cupy_allclose(accept_error=TypeError)
+    def test_bincount_zero(self, xp, dtype):
+        x = testing.shaped_arange((3,), xp, dtype)
+        return xp.bincount(x, minlength=0)
+
+    @for_all_dtypes_bincount()
+    def test_bincount_too_small_minlength(self, dtype):
+        for xp in (numpy, cupy):
+            x = testing.shaped_arange((3,), xp, dtype)
+            # TODO(imanishi): Split this test into a test for ValueError and
+            # a test for TypeError.
+            with pytest.raises((ValueError, TypeError)):
+                xp.bincount(x, minlength=-1)
+
+
+# This class compares CUB results against NumPy's
+# @testing.gpu
+# @unittest.skipUnless(cupy.cuda.cub.available, 'The CUB routine is not enabled')
+# class TestCubHistogram(unittest.TestCase):
+
+    # def setUp(self):
+        # self.old_accelerators = _accelerator.get_routine_accelerators()
+        # _accelerator.set_routine_accelerators(['cub'])
+
+    # def tearDown(self):
+        # _accelerator.set_routine_accelerators(self.old_accelerators)
+
+    # @testing.for_all_dtypes(no_bool=True, no_complex=True)
+    # @testing.numpy_cupy_array_equal()
+    # def test_histogram(self, xp, dtype):
+        # x = testing.shaped_arange((10,), xp, dtype)
+
+        # if xp is numpy:
+            # return xp.histogram(x)
+
+        # # xp is cupy, first ensure we really use CUB
+        # cub_func = 'cupy._statistics.histogram.cub.device_histogram'
+        # with testing.AssertFunctionIsCalled(cub_func):
+            # xp.histogram(x)
+        # # ...then perform the actual computation
+        # return xp.histogram(x)
+
+    # @testing.for_all_dtypes(no_bool=True, no_complex=True)
+    # @testing.numpy_cupy_array_equal()
+    # def test_histogram_range_float(self, xp, dtype):
+        # a = testing.shaped_arange((10,), xp, dtype)
+        # h, b = xp.histogram(a, testing.shaped_arange((10,), xp, numpy.float64))
+        # assert int(h.sum()) == 10
+        # return h, b
+
+
+@testing.gpu
+@testing.parameterize(*testing.product(
+    {'bins': [
+        # Test monotonically increasing with in-bounds values
+        [1.5, 2.5, 4.0, 6.0],
+        # Explicit out-of-bounds for x values
+        [-1.0, 1.0, 2.5, 4.0, 20.0],
+        # Repeated values should yield right-most or left-most indexes
+        [0.0, 1.0, 1.0, 4.0, 4.0, 10.0],
+    ],
+        'increasing': [True, False],
+        'right': [True, False],
+        'shape': [(), (10,), (6, 3, 3)]})
+)
+class TestDigitize(unittest.TestCase):
+
+    @testing.for_all_dtypes(no_bool=True, no_complex=True)
+    @testing.numpy_cupy_array_equal()
+    def test_digitize(self, xp, dtype):
+        x = testing.shaped_arange(self.shape, xp, dtype)
+        bins = self.bins
+        if not self.increasing:
+            bins = bins[::-1]
+        bins = xp.array(bins)
+        y = xp.digitize(x, bins, right=self.right)
+        return y,
+
+
+@testing.gpu
+@testing.parameterize(
+    {'right': True},
+    {'right': False})
+class TestDigitizeNanInf(unittest.TestCase):
+
+    @testing.numpy_cupy_array_equal()
+    def test_digitize_nan(self, xp):
+        x = testing.shaped_arange((14,), xp, xp.float32)
+        x[5] = float('nan')
+        bins = xp.array([1.0, 3.0, 5.0, 8.0, 12.0], xp.float32)
+        y = xp.digitize(x, bins, right=self.right)
+        return y,
+
+    @testing.numpy_cupy_array_equal()
+    def test_digitize_nan_bins(self, xp):
+        x = testing.shaped_arange((14,), xp, xp.float32)
+        bins = xp.array([1.0, 3.0, 5.0, 8.0, float('nan')], xp.float32)
+        y = xp.digitize(x, bins, right=self.right)
+        return y,
+
+    @testing.numpy_cupy_array_equal()
+    def test_digitize_nan_bins_repeated(self, xp):
+        x = testing.shaped_arange((14,), xp, xp.float32)
+        x[5] = float('nan')
+        bins = [1.0, 3.0, 5.0, 8.0, float('nan'), float('nan')]
+        bins = xp.array(bins, xp.float32)
+        y = xp.digitize(x, bins, right=self.right)
+        return y,
+
+    @testing.numpy_cupy_array_equal()
+    def test_digitize_nan_bins_decreasing(self, xp):
+        x = testing.shaped_arange((14,), xp, xp.float32)
+        x[5] = float('nan')
+        bins = [float('nan'), 8.0, 5.0, 3.0, 1.0]
+        bins = xp.array(bins, xp.float32)
+        y = xp.digitize(x, bins, right=self.right)
+        return y,
+
+    @testing.numpy_cupy_array_equal()
+    def test_digitize_nan_bins_decreasing_repeated(self, xp):
+        x = testing.shaped_arange((14,), xp, xp.float32)
+        x[5] = float('nan')
+        bins = [float('nan'), float('nan'), float('nan'), 5.0, 3.0, 1.0]
+        bins = xp.array(bins, xp.float32)
+        y = xp.digitize(x, bins, right=self.right)
+        return y,
+
+    @testing.numpy_cupy_array_equal()
+    def test_digitize_all_nan_bins(self, xp):
+        x = testing.shaped_arange((14,), xp, xp.float32)
+        x[5] = float('nan')
+        bins = [float('nan'), float('nan'), float('nan'), float('nan')]
+        bins = xp.array(bins, xp.float32)
+        y = xp.digitize(x, bins, right=self.right)
+        return y,
+
+    @testing.numpy_cupy_array_equal()
+    def test_searchsorted_inf(self, xp):
+        x = testing.shaped_arange((14,), xp, xp.float64)
+        x[5] = float('inf')
+        bins = xp.array([0, 1, 2, 4, 10])
+        y = xp.digitize(x, bins, right=self.right)
+        return y,
+
+    @testing.numpy_cupy_array_equal()
+    def test_searchsorted_minf(self, xp):
+        x = testing.shaped_arange((14,), xp, xp.float64)
+        x[5] = float('-inf')
+        bins = xp.array([0, 1, 2, 4, 10])
+        y = xp.digitize(x, bins, right=self.right)
+        return y,
+
+
+@testing.gpu
+class TestDigitizeInvalid(unittest.TestCase):
+
+    def test_digitize_complex(self):
+        for xp in (numpy, cupy):
+            x = testing.shaped_arange((14,), xp, xp.complex)
+            bins = xp.array([1.0, 3.0, 5.0, 8.0, 12.0], xp.complex)
+            with pytest.raises(TypeError):
+                xp.digitize(x, bins)
+
+    def test_digitize_nd_bins(self):
+        for xp in (numpy, cupy):
+            x = testing.shaped_arange((14,), xp, xp.float64)
+            bins = xp.array([[1], [2]])
+            with pytest.raises(ValueError):
+                xp.digitize(x, bins)
diff --git a/tests/third_party/cupy/statistics_tests/test_mean.py b/tests/third_party/cupy/statics_tests/test_meanvar.py
similarity index 79%
rename from tests/third_party/cupy/statistics_tests/test_mean.py
rename to tests/third_party/cupy/statics_tests/test_meanvar.py
index d72a51c4dd15..ce6953812bbe 100644
--- a/tests/third_party/cupy/statistics_tests/test_mean.py
+++ b/tests/third_party/cupy/statics_tests/test_meanvar.py
@@ -10,6 +10,78 @@
     "ignore", category=RuntimeWarning)
 
 
+@testing.gpu
+class TestMedian(unittest.TestCase):
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_allclose()
+    def test_median_noaxis(self, xp, dtype):
+        a = testing.shaped_random((3, 4, 5), xp, dtype)
+        return xp.median(a)
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_allclose()
+    def test_median_axis1(self, xp, dtype):
+        a = testing.shaped_random((3, 4, 5), xp, dtype)
+        return xp.median(a, axis=1)
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_allclose()
+    def test_median_axis2(self, xp, dtype):
+        a = testing.shaped_random((3, 4, 5), xp, dtype)
+        return xp.median(a, axis=2)
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_allclose()
+    def test_median_overwrite_input(self, xp, dtype):
+        a = testing.shaped_random((3, 4, 5), xp, dtype)
+        return xp.median(a, overwrite_input=True)
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_allclose()
+    def test_median_keepdims_axis1(self, xp, dtype):
+        a = testing.shaped_random((3, 4, 5), xp, dtype)
+        return xp.median(a, axis=1, keepdims=True)
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_allclose()
+    def test_median_keepdims_noaxis(self, xp, dtype):
+        a = testing.shaped_random((3, 4, 5), xp, dtype)
+        return xp.median(a, keepdims=True)
+
+    def test_median_invalid_axis(self):
+        for xp in [numpy, cupy]:
+            a = testing.shaped_random((3, 4, 5), xp)
+            with pytest.raises(numpy.AxisError):
+                return xp.median(a, -a.ndim - 1, keepdims=False)
+
+            with pytest.raises(numpy.AxisError):
+                return xp.median(a, a.ndim, keepdims=False)
+
+            with pytest.raises(numpy.AxisError):
+                return xp.median(a, (-a.ndim - 1, 1), keepdims=False)
+
+            with pytest.raises(numpy.AxisError):
+                return xp.median(a, (0, a.ndim,), keepdims=False)
+
+
+@testing.parameterize(
+    *testing.product({
+        'shape': [(3, 4, 5)],
+        'axis': [(0, 1), (0, -1), (1, 2), (1,)],
+        'keepdims': [True, False]
+    })
+)
+@testing.gpu
+class TestMedianAxis(unittest.TestCase):
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_allclose()
+    def test_median_axis_sequence(self, xp, dtype):
+        a = testing.shaped_random(self.shape, xp, dtype)
+        return xp.median(a, self.axis, keepdims=self.keepdims)
+
+
 @testing.gpu
 class TestAverage(unittest.TestCase):
 
@@ -88,11 +160,24 @@ def test_external_mean_axis(self, xp, dtype):
         a = testing.shaped_arange((2, 3, 4), xp, dtype)
         return xp.mean(a, axis=1)
 
+    @testing.for_all_dtypes(no_complex=True)
     @testing.numpy_cupy_allclose()
-    def test_mean_all_dtype(self, xp):
-        a = xp.full((2, 3, 4), 123456789, dtype=numpy.int64)
+    def test_mean_all_float64_dtype(self, xp, dtype):
+        a = xp.full((2, 3, 4), 123456789, dtype=dtype)
         return xp.mean(a, dtype=numpy.float64)
 
+    @testing.for_all_dtypes(no_complex=True)
+    @testing.numpy_cupy_allclose()
+    def test_mean_all_int64_dtype(self, xp, dtype):
+        a = testing.shaped_arange((2, 3, 4), xp, dtype)
+        return xp.mean(a, dtype=numpy.int64)
+
+    @testing.for_all_dtypes()
+    @testing.numpy_cupy_allclose()
+    def test_mean_all_complex_dtype(self, xp, dtype):
+        a = testing.shaped_arange((2, 3, 4), xp, dtype)
+        return xp.mean(a, dtype=numpy.complex64)
+
     @testing.for_all_dtypes()
     @testing.numpy_cupy_allclose()
     def test_var_all(self, xp, dtype):
diff --git a/tests/third_party/cupy/statics_tests/test_order.py b/tests/third_party/cupy/statics_tests/test_order.py
new file mode 100644
index 000000000000..28785181a2b2
--- /dev/null
+++ b/tests/third_party/cupy/statics_tests/test_order.py
@@ -0,0 +1,246 @@
+import unittest
+import warnings
+
+import numpy
+import pytest
+
+import dpnp as cupy
+from tests.third_party.cupy import testing
+
+
+_all_interpolations = (
+    'lower',
+    'higher',
+    'midpoint',
+    # 'nearest', # TODO(hvy): Not implemented
+    'linear')
+
+
+def for_all_interpolations(name='interpolation'):
+    return testing.for_orders(_all_interpolations, name=name)
+
+
+@testing.gpu
+class TestOrder(unittest.TestCase):
+
+    @for_all_interpolations()
+    @testing.for_all_dtypes(no_float16=True, no_bool=True, no_complex=True)
+    @testing.numpy_cupy_allclose()
+    def test_percentile_defaults(self, xp, dtype, interpolation):
+        a = testing.shaped_random((2, 3, 8), xp, dtype)
+        q = testing.shaped_random((3,), xp, dtype=dtype, scale=100)
+        return xp.percentile(a, q, interpolation=interpolation)
+
+    @for_all_interpolations()
+    @testing.for_all_dtypes(no_float16=True, no_bool=True, no_complex=True)
+    @testing.numpy_cupy_allclose()
+    def test_percentile_q_list(self, xp, dtype, interpolation):
+        a = testing.shaped_arange((1001,), xp, dtype)
+        q = [99, 99.9]
+        return xp.percentile(a, q, interpolation=interpolation)
+
+    @for_all_interpolations()
+    @testing.for_all_dtypes(no_float16=True, no_bool=True, no_complex=True)
+    @testing.numpy_cupy_allclose(rtol=1e-6)
+    def test_percentile_no_axis(self, xp, dtype, interpolation):
+        a = testing.shaped_random((10, 2, 4, 8), xp, dtype)
+        q = testing.shaped_random((5,), xp, dtype=dtype, scale=100)
+        return xp.percentile(a, q, axis=None, interpolation=interpolation)
+
+    @for_all_interpolations()
+    @testing.for_all_dtypes(no_float16=True, no_bool=True, no_complex=True)
+    @testing.numpy_cupy_allclose(rtol=1e-6)
+    def test_percentile_neg_axis(self, xp, dtype, interpolation):
+        a = testing.shaped_random((4, 3, 10, 2, 8), xp, dtype)
+        q = testing.shaped_random((5,), xp, dtype=dtype, scale=100)
+        return xp.percentile(a, q, axis=-1, interpolation=interpolation)
+
+    @for_all_interpolations()
+    @testing.for_all_dtypes(no_float16=True, no_bool=True, no_complex=True)
+    @testing.numpy_cupy_allclose(rtol=1e-6)
+    def test_percentile_tuple_axis(self, xp, dtype, interpolation):
+        a = testing.shaped_random((1, 6, 3, 2), xp, dtype)
+        q = testing.shaped_random((5,), xp, dtype=dtype, scale=100)
+        return xp.percentile(a, q, axis=(0, 1, 2), interpolation=interpolation)
+
+    @for_all_interpolations()
+    @testing.for_all_dtypes(no_float16=True, no_bool=True, no_complex=True)
+    @testing.numpy_cupy_allclose()
+    def test_percentile_scalar_q(self, xp, dtype, interpolation):
+        a = testing.shaped_random((2, 3, 8), xp, dtype)
+        q = 13.37
+        return xp.percentile(a, q, interpolation=interpolation)
+
+    @for_all_interpolations()
+    @testing.for_all_dtypes(no_float16=True, no_bool=True, no_complex=True)
+    @testing.numpy_cupy_allclose(rtol=1e-5)
+    def test_percentile_keepdims(self, xp, dtype, interpolation):
+        a = testing.shaped_random((7, 2, 9, 2), xp, dtype)
+        q = testing.shaped_random((5,), xp, dtype=dtype, scale=100)
+        return xp.percentile(
+            a, q, axis=None, keepdims=True, interpolation=interpolation)
+
+    @for_all_interpolations()
+    @testing.for_float_dtypes(no_float16=True)  # NumPy raises error on int8
+    @testing.numpy_cupy_allclose(rtol=1e-6)
+    def test_percentile_out(self, xp, dtype, interpolation):
+        a = testing.shaped_random((10, 2, 3, 2), xp, dtype)
+        q = testing.shaped_random((5,), xp, dtype=dtype, scale=100)
+        out = testing.shaped_random((5, 10, 2, 3), xp, dtype)
+        return xp.percentile(
+            a, q, axis=-1, interpolation=interpolation, out=out)
+
+    @for_all_interpolations()
+    @testing.for_all_dtypes(no_float16=True, no_bool=True, no_complex=True)
+    def test_percentile_bad_q(self, dtype, interpolation):
+        for xp in (numpy, cupy):
+            a = testing.shaped_random((4, 2, 3, 2), xp, dtype)
+            q = testing.shaped_random((1, 2, 3), xp, dtype=dtype, scale=100)
+            with pytest.raises(ValueError):
+                xp.percentile(a, q, axis=-1, interpolation=interpolation)
+
+    @testing.for_all_dtypes(no_float16=True, no_bool=True, no_complex=True)
+    def test_percentile_uxpected_interpolation(self, dtype):
+        for xp in (numpy, cupy):
+            a = testing.shaped_random((4, 2, 3, 2), xp, dtype)
+            q = testing.shaped_random((5,), xp, dtype=dtype, scale=100)
+            with pytest.raises(ValueError):
+                xp.percentile(a, q, axis=-1, interpolation='deadbeef')
+
+    @testing.for_all_dtypes(no_complex=True)
+    @testing.numpy_cupy_allclose()
+    def test_nanmax_all(self, xp, dtype):
+        a = testing.shaped_random((2, 3), xp, dtype)
+        return xp.nanmax(a)
+
+    @testing.for_all_dtypes(no_complex=True)
+    @testing.numpy_cupy_allclose()
+    def test_nanmax_axis_large(self, xp, dtype):
+        a = testing.shaped_random((3, 1000), xp, dtype)
+        return xp.nanmax(a, axis=0)
+
+    @testing.for_all_dtypes(no_complex=True)
+    @testing.numpy_cupy_allclose()
+    def test_nanmax_axis0(self, xp, dtype):
+        a = testing.shaped_random((2, 3, 4), xp, dtype)
+        return xp.nanmax(a, axis=0)
+
+    @testing.for_all_dtypes(no_complex=True)
+    @testing.numpy_cupy_allclose()
+    def test_nanmax_axis1(self, xp, dtype):
+        a = testing.shaped_random((2, 3, 4), xp, dtype)
+        return xp.nanmax(a, axis=1)
+
+    @testing.for_all_dtypes(no_complex=True)
+    @testing.numpy_cupy_allclose()
+    def test_nanmax_axis2(self, xp, dtype):
+        a = testing.shaped_random((2, 3, 4), xp, dtype)
+        return xp.nanmax(a, axis=2)
+
+    @testing.for_float_dtypes()
+    @testing.numpy_cupy_allclose()
+    def test_nanmax_nan(self, xp, dtype):
+        a = xp.array([float('nan'), 1, -1], dtype)
+        with warnings.catch_warnings():
+            return xp.nanmax(a)
+
+    @testing.for_float_dtypes()
+    @testing.numpy_cupy_allclose()
+    def test_nanmax_all_nan(self, xp, dtype):
+        a = xp.array([float('nan'), float('nan')], dtype)
+        with warnings.catch_warnings(record=True) as w:
+            warnings.simplefilter('always')
+            m = xp.nanmax(a)
+        self.assertEqual(len(w), 1)
+        self.assertIs(w[0].category, RuntimeWarning)
+        return m
+
+    @testing.for_all_dtypes(no_complex=True)
+    @testing.numpy_cupy_allclose()
+    def test_nanmin_all(self, xp, dtype):
+        a = testing.shaped_random((2, 3), xp, dtype)
+        return xp.nanmin(a)
+
+    @testing.for_all_dtypes(no_complex=True)
+    @testing.numpy_cupy_allclose()
+    def test_nanmin_axis_large(self, xp, dtype):
+        a = testing.shaped_random((3, 1000), xp, dtype)
+        return xp.nanmin(a, axis=0)
+
+    @testing.for_all_dtypes(no_complex=True)
+    @testing.numpy_cupy_allclose()
+    def test_nanmin_axis0(self, xp, dtype):
+        a = testing.shaped_random((2, 3, 4), xp, dtype)
+        return xp.nanmin(a, axis=0)
+
+    @testing.for_all_dtypes(no_complex=True)
+    @testing.numpy_cupy_allclose()
+    def test_nanmin_axis1(self, xp, dtype):
+        a = testing.shaped_random((2, 3, 4), xp, dtype)
+        return xp.nanmin(a, axis=1)
+
+    @testing.for_all_dtypes(no_complex=True)
+    @testing.numpy_cupy_allclose()
+    def test_nanmin_axis2(self, xp, dtype):
+        a = testing.shaped_random((2, 3, 4), xp, dtype)
+        return xp.nanmin(a, axis=2)
+
+    @testing.for_float_dtypes()
+    @testing.numpy_cupy_allclose()
+    def test_nanmin_nan(self, xp, dtype):
+        a = xp.array([float('nan'), 1, -1], dtype)
+        with warnings.catch_warnings():
+            return xp.nanmin(a)
+
+    @testing.for_float_dtypes()
+    @testing.numpy_cupy_allclose()
+    def test_nanmin_all_nan(self, xp, dtype):
+        a = xp.array([float('nan'), float('nan')], dtype)
+        with warnings.catch_warnings(record=True) as w:
+            warnings.simplefilter('always')
+            m = xp.nanmin(a)
+        self.assertEqual(len(w), 1)
+        self.assertIs(w[0].category, RuntimeWarning)
+        return m
+
+    @testing.for_all_dtypes(no_bool=True)
+    @testing.numpy_cupy_allclose()
+    def test_ptp_all(self, xp, dtype):
+        a = testing.shaped_random((2, 3), xp, dtype)
+        return xp.ptp(a)
+
+    @testing.for_all_dtypes(no_bool=True)
+    @testing.numpy_cupy_allclose()
+    def test_ptp_axis_large(self, xp, dtype):
+        a = testing.shaped_random((3, 1000), xp, dtype)
+        return xp.ptp(a, axis=0)
+
+    @testing.for_all_dtypes(no_bool=True)
+    @testing.numpy_cupy_allclose()
+    def test_ptp_axis0(self, xp, dtype):
+        a = testing.shaped_random((2, 3, 4), xp, dtype)
+        return xp.ptp(a, axis=0)
+
+    @testing.for_all_dtypes(no_bool=True)
+    @testing.numpy_cupy_allclose()
+    def test_ptp_axis1(self, xp, dtype):
+        a = testing.shaped_random((2, 3, 4), xp, dtype)
+        return xp.ptp(a, axis=1)
+
+    @testing.for_all_dtypes(no_bool=True)
+    @testing.numpy_cupy_allclose()
+    def test_ptp_axis2(self, xp, dtype):
+        a = testing.shaped_random((2, 3, 4), xp, dtype)
+        return xp.ptp(a, axis=2)
+
+    @testing.for_float_dtypes()
+    @testing.numpy_cupy_allclose()
+    def test_ptp_nan(self, xp, dtype):
+        a = xp.array([float('nan'), 1, -1], dtype)
+        return xp.ptp(a)
+
+    @testing.for_float_dtypes()
+    @testing.numpy_cupy_allclose()
+    def test_ptp_all_nan(self, xp, dtype):
+        a = xp.array([float('nan'), float('nan')], dtype)
+        return xp.ptp(a)
diff --git a/tests/third_party/cupy/testing/__init__.py b/tests/third_party/cupy/testing/__init__.py
index 97e0442b93cb..2ef543c4bacc 100644
--- a/tests/third_party/cupy/testing/__init__.py
+++ b/tests/third_party/cupy/testing/__init__.py
@@ -14,7 +14,7 @@
 from tests.third_party.cupy.testing.attr import gpu
 # from tests.third_party.cupy.testing.attr import multi_gpu
 from tests.third_party.cupy.testing.attr import slow
-# from tests.third_party.cupy.testing.helper import assert_warns
+from tests.third_party.cupy.testing.helper import assert_warns
 # from tests.third_party.cupy.testing.helper import empty
 from tests.third_party.cupy.testing.helper import for_all_dtypes
 from tests.third_party.cupy.testing.helper import for_all_dtypes_combination