diff --git a/light-curve/light_curve/light_curve_py/features/ptp_var.py b/light-curve/light_curve/light_curve_py/features/ptp_var.py
index 2a4f248c..68a363a4 100644
--- a/light-curve/light_curve/light_curve_py/features/ptp_var.py
+++ b/light-curve/light_curve/light_curve_py/features/ptp_var.py
@@ -9,9 +9,8 @@ class PeakToPeakVar(BaseSingleBandFeature):
     $$
     \frac{(m_i - \sigma_i)_\text{max} - (m_i + \sigma_i)_\text{min}}{(m_i - \sigma_i)_\text{max} + (m_i + \sigma_i)_\text{min}}
     $$
-    For non-variable data, it should be close to zero.
-    If data is close to be variable, the index should be more or equal than 0.10-0.15.
-    It is sensitive to magnitude of error values and can be negative in overestimated errors case.
+    Input m must be non-negative (e.g. non-differential) flux density. This feature is a variability detector, higher values correspond to more variable sources. 
+    For example, if observational count is 10^2-10^4 and Signal to Noise > 7, the feature is larger than 0.1-0.15, these value are commonly used as a variability threshold.
     
     - Depends on: **flux density**, **errors**
     - Minimum number of observations: **2**
@@ -19,10 +18,13 @@ class PeakToPeakVar(BaseSingleBandFeature):
 
     Aller M.F., Aller H.D., Hughes P.A. 1992. [DOI:10.1086/171898](https://www.doi.org/10.1086/171898)
     """
+    nstd: float = 1.0
 
     def _eval_single_band(self, t, m, sigma=None):
-        a = np.max(m - sigma)
-        b = np.min(m + sigma)
+        if np.any(m < 0): 
+            raise ValueError("m must be non-negative")
+        a = np.max(m - self.nstd*sigma)
+        b = np.min(m + self.nstd*sigma)
         return (a - b) / (a + b)
 
     @property
diff --git a/light-curve/tests/light_curve_py/features/test_ptp_var.py b/light-curve/tests/light_curve_py/features/test_ptp_var.py
index 515687df..c4834cb2 100644
--- a/light-curve/tests/light_curve_py/features/test_ptp_var.py
+++ b/light-curve/tests/light_curve_py/features/test_ptp_var.py
@@ -4,71 +4,18 @@
 from light_curve.light_curve_py import PeakToPeakVar
 
 
-def test_ptpvar_const_data():
+def test_ptpvar():
     feature = PeakToPeakVar()
-    n = 100
-    t = np.arange(n)
-    m = np.ones_like(t)
-    sigma = 0.1 * np.ones_like(t)
-    actual = feature(t, m, sigma)
-    desired = -0.1
-    assert_allclose(actual, desired)
-
-
-def test_ptpvar_periodic_data():
-    feature = PeakToPeakVar()
-    n = 100
-    t = np.linspace(0, np.pi, n)
-    m = np.sin(t)
-    sigma = 0.1 * np.ones_like(t)
-    actual = feature(t, m, sigma)
-    desired = 0.8
-    assert_allclose(actual, desired, rtol=3 / np.sqrt(n))
-
-
-def test_ptpvar_norm_data_1():
     rng = np.random.default_rng(0)
     n = 100
-    t = np.linspace(0, 1, n)
-    m = np.abs(rng.normal(0, 1, n))
-    sigma = 0.2 * np.ones_like(t)
-    feature = PeakToPeakVar()
-    actual = feature(t, m, sigma)
-    desired = 1
-    assert_allclose(actual, desired, rtol=3 / np.sqrt(n))
-
-
-def test_ptpvar_norm_data_2():
-    rng = np.random.default_rng(0)
-    n = 10000
-    t = np.linspace(0, 1, n)
-    m = np.abs(rng.normal(0, 1, n))
-    sigma = 0.2 * np.ones_like(t)
-    feature = PeakToPeakVar()
-    actual = feature(t, m, sigma)
-    desired = 1
-    assert_allclose(actual, desired, rtol=3 / np.sqrt(n))
-
-
-def test_ptpvar_expon_data_1():
-    rng = np.random.default_rng(0)
-    n = 100
-    t = np.linspace(0, 1, n)
-    m = rng.exponential(2, n)
-    sigma = np.ones_like(t)
-    feature = PeakToPeakVar()
-    actual = feature(t, m, sigma)
-    desired = 1
-    assert_allclose(actual, desired, rtol=3 / np.sqrt(n))
-
-
-def test_ptpvar_expon_data_2():
-    rng = np.random.default_rng(0)
-    n = 10000
-    t = np.linspace(0, 1, n)
-    m = rng.exponential(2, n)
-    sigma = np.ones_like(t)
-    feature = PeakToPeakVar()
-    actual = feature(t, m, sigma)
-    desired = 1
-    assert_allclose(actual, desired, rtol=3 / np.sqrt(n))
\ No newline at end of file
+    m = 10000
+    lmd = 200
+    t = np.arange(n)
+    flux_list = np.random.poisson(lmd, (m, n))
+    ptp_list = [PeakToPeakVar(t, flux_list[i], np.sqrt(flux_list[i])) for i in range(m)]
+    flux = rng.poisson(lmd, n)
+    sigma = np.sqrt(flux)
+    actual = feature(t, flux, sigma)
+    desired = (np.quantile(ptp_list, 0.025) + np.quantile(ptp_list, 0.975)) / 2
+    atol = desired - np.quantile(ptp_list, 0.025)
+    assert_allclose(actual, desired, atol=atol)
\ No newline at end of file