diff --git a/source/statistics/code/ztest_pvalue.R b/source/statistics/code/ztest_pvalue.R
new file mode 100644
index 0000000..4396556
--- /dev/null
+++ b/source/statistics/code/ztest_pvalue.R
@@ -0,0 +1,10 @@
+
+
+N=rnorm(9999999,mean=0,sd=1)
+cdf=ecdf(N)
+zscore=-1.8
+
+message(paste0("Alpha approximated is ",cdf(zscore)))
+message(paste0("Alpha from built-in CDF ",pnorm(zscore,mean=0,sd=1)))
+message(paste0("Bonus: zscore from p-value ",qnorm(0.02275,mean=0,sd=1)))
+
diff --git a/source/statistics/tests/parametric/ztest.rst b/source/statistics/tests/parametric/ztest.rst
index 20dc042..85bff69 100644
--- a/source/statistics/tests/parametric/ztest.rst
+++ b/source/statistics/tests/parametric/ztest.rst
@@ -48,13 +48,31 @@ That *p-value* is computed as follow:
 
 .. image:: ../../figures/normal_law_tails.svg
    :align: center
-
+      
 If a z-test is done over one tail (left or right) it is called a **one-tailed** z-test.
 If a z-test is done over both tails (left and right) it is called a **two-tailed** z-test.
 
+The following code shows you how to obtain the *p-value* in R:
+
+.. literalinclude:: ../../code/ztest_pvalue.R
+   :language: R
+
+Output example:
+
+.. code-block:: console
+                
+   Alpha approximated is 0.0359588035958804
+   Alpha from built-in CDF 0.0359303191129258
+   Bonus: zscore from p-value -2.0000024438996
+
 If the :math:`\alpha` value given by the test is lower or equal to the *p-value* threshold chosen prior the test,
 :math:`H_0` is rejected and :math:`H_1` is considered accepted.
 
+Using the z-score from *p-value* function allows you to create a *rejection region*.
+It is an interval :math:`[\mathrm{zscore}_{min},\mathrm{zscore}_{max}]`.
+If the sample z-score is not part of the interval, then :math:`H_0` is rejected and :math:`H_1` is considered accepted.
+
+
 One-tailed vs Two-tailed
 ========================
 
diff --git a/test.R b/test.R
deleted file mode 100644
index d6ae433..0000000
--- a/test.R
+++ /dev/null
@@ -1,8 +0,0 @@
-
-
-cdf=ecdf(N)
-initial_alpha=0.05
-
-print(round(pnorm(1.01,mean=0,sd=1),digits=4))
-
-