From 2e14b12850c42ee34fd1fd2f3fcfe291273f9953 Mon Sep 17 00:00:00 2001
From: peisenha <peisenha@protonmail.com>
Date: Fri, 19 Jun 2020 12:29:52 +0200
Subject: [PATCH 01/45] added next steps

---
 .../tutorials/conditional_sampling.ipynb      | 89 +++++++++++++++++++
 1 file changed, 89 insertions(+)
 create mode 100644 docs/source/tutorials/conditional_sampling.ipynb

diff --git a/docs/source/tutorials/conditional_sampling.ipynb b/docs/source/tutorials/conditional_sampling.ipynb
new file mode 100644
index 00000000..3afd6afe
--- /dev/null
+++ b/docs/source/tutorials/conditional_sampling.ipynb
@@ -0,0 +1,89 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Conditional sampling"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can rely on `chaospy`'s features to construct flexible multivariate distributions for input paramters using copulas. This allows us to sample from the **joint distribution** of input paramters as is required for uncertainty propagation. However, we now also want to be able to conduct a quantitative sensitivity analysis using, for example, Sobol indicies (**@timmens**) an Shapley values (**@lindamaok899**). This remains a straightforward task where the multivariate distribution is known to be Gaussian as analytical expressions exist ([Wikipedia](https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Conditional_distributions)).\n",
+    "\n",
+    "We need the capablity to construct these **conditional distributions** in a more flexible way to, for example, compute the Sobol indices and Shapley values for the EQO model where the marginal distribution of the input parameters is uniform instead of normal.\n",
+    "\n",
+    "Again, we will pursue a copula approach and use a Gaussian copula to generate the samples."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## $1^{st}$ Step: Flexible sampler for conditional Gaussian distribution"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "There appears to be no flexible Python code out there as part of a sound library that simply samples from the conditional normal distribution (please confirm).  We do have some work here already done by Janos and Tim, please reach out to them and see what can be useful for us. \n",
+    "\n",
+    "However, there is an implementation available to construct the conditional mean and variance in $R$ as part of the [condMVNorm](https://cran.r-project.org/web/packages/condMVNorm/condMVNorm.pdf) package. Please construct an analogous Python implementation.\n",
+    "\n",
+    "Be sure to validate your implementation against the original implementation. You can probably use the function in Python as well by using [rpy2](https://rpy2.github.io/). Simply set up a loop that specifies numerous multivariate normal distributions, specifies a request for a conditional distribution, and then compares the output between the two implementations.\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## $2^{nd}$ Step: Gaussian copula with flexible marginals"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We need to extend our conditional sampler to allow for more flexible marginal distributions. This boils down to transforming the quantile from the marginal distribution into the corresponding value of the Gaussian. We use the special case that is equivalent to a simple multivariate normal distribution for testing."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## $3^{rd}$ Step: Mapping correlation matrix"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We need to map the correlation between the orginal variables to the implied correlation for the Gaussian copula. Here, Section 4 in this [paper](https://doi.org/10.1016/j.cpc.2011.12.020) is a good starting point."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

From fc2cab49844a485a09bcc284333b6ec33fccf885 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Tue, 23 Jun 2020 18:35:39 +0200
Subject: [PATCH 02/45] add condMVN function

---
 econsa/sampling.py | 101 +++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 101 insertions(+)

diff --git a/econsa/sampling.py b/econsa/sampling.py
index 6b2228cd..aa5435d9 100644
--- a/econsa/sampling.py
+++ b/econsa/sampling.py
@@ -3,7 +3,108 @@
 This module contains all we need to sample random input parameters.
 
 """
+import numpy as np
 
 
 def get_sample(*args, **kwargs):
     return
+
+
+def condMVN(mean, sigma, dependent_ind, given_ind=None, x_given=None, check_sigma=True):
+    r"""Conditional mean and variance function.
+
+    Returns conditional mean and variance of dependent variables,
+    given multivariate normal distribution and indices of independent variables.
+
+    This is a translation of the main function of R package "condMVNorm".
+
+    .. math::
+        X = (X_{\text{ind}}, X_{\text{dep}}) \sim \mathcal{N}
+
+    Parameters
+    ----------
+    mean : array_like
+           The mean vector of the multivariate normal distribution.
+
+    sigma : array_like
+            Symmetric and positive-definte covariance matrix of
+            the multivariate normal distribution.
+
+    dependent_ind : array_like
+                    The indices of dependent variables.
+
+    given_ind : array_like, optional
+                The indices of independent variables (default value is ``None``).
+                If not speficied or all values are zero, return unconditional values.
+
+    x_given : array_like, optional
+              The conditioning values (default value is ``None``).
+              Should be the same length as `given_ind`, otherwise throw an error.
+
+    check_sigma : boolean, optional
+                  Check that `sigma` is symmetric,
+                  and no eigenvalue is zero (default value is ``True``).
+
+    Returns
+    -------
+    cond_mean : array_like
+               The conditional mean of dependent variables.
+
+    cond_var : arrray_like
+              The conditional covariance matrix of dependent variables.
+
+    Examples
+    --------
+    >>> mean = np.array([1, 1, 1])
+    >>> sigma = np.array([[4.0677098, -0.9620331, 0.9897267],
+    ...                   [-0.9620331, 2.2775449, 0.7475968],
+    ...                   [0.9897267, 0.7475968, 0.7336631]])
+    >>> dependent_ind = [0,]
+    >>> given_ind = [1, 2]
+    >>> x_given = [1, -1]
+    >>> cond_mean, cond_var = condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+    >>> np.testing.assert_almost_equal(cond_mean, -4.347531, decimal=6)
+    >>> np.testing.assert_almost_equal(cond_var, 0.170718, decimal=6)
+    """
+    #
+    mean_np = np.array(mean).squeeze()
+    sigma_np = np.array(sigma).squeeze()
+
+    if given_ind is None and x_given is None:
+        condMean = np.array(mean_np[dependent_ind])
+        condVar = np.array(sigma_np[dependent_ind, :][:, dependent_ind])
+        return (condMean, condVar)
+
+    # Make sure that given_len is not empty:
+    given_ind_np = np.array(given_ind).squeeze()
+    try:
+        len(given_ind_np)
+    except TypeError:
+        condMean = np.array(mean_np[dependent_ind])
+        condVar = np.array(sigma_np[dependent_ind, :][:, dependent_ind])
+        return (condMean, condVar)
+
+    # Make sure that x_given is not empty & aligns with given_len:
+    x_given_np = np.array(x_given).squeeze()
+    try:
+        if len(x_given_np) != len(given_ind_np):
+            raise ValueError("lengths of x_given and given_ind must be the same")
+    except TypeError:
+        raise
+
+    # Check sigma is symmetric & positive-definite:
+    if check_sigma:
+        if not np.allclose(sigma_np, sigma_np.T):
+            raise ValueError("sigma is not a symmetric matrix")
+        elif not np.all(np.linalg.eigvals(sigma) > 0):
+            raise ValueError("sigma is not positive-definite")
+
+    b = sigma_np[dependent_ind, :][:, dependent_ind]
+    c = sigma_np[dependent_ind, :][:, given_ind]
+    d = sigma_np[given_ind, :][:, given_ind]
+    c_dinv = c @ np.linalg.inv(d)
+
+    cond_mean = mean_np[dependent_ind] + c_dinv @ (x_given - mean_np[given_ind])
+    cond_var = b - c_dinv @ (c.T)
+
+    return (cond_mean, cond_var)

From e45da3f964c0a6c81596ca6860f011f20403c536 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Tue, 23 Jun 2020 19:03:29 +0200
Subject: [PATCH 03/45] add condMVN test

---
 econsa/sampling.py            |  6 ++--
 econsa/tests/test_sampling.py | 68 ++++++++++++++++++++++++++++++++++-
 environment.yml               |  1 +
 3 files changed, 71 insertions(+), 4 deletions(-)

diff --git a/econsa/sampling.py b/econsa/sampling.py
index aa5435d9..f9cdc8c1 100644
--- a/econsa/sampling.py
+++ b/econsa/sampling.py
@@ -27,7 +27,7 @@ def condMVN(mean, sigma, dependent_ind, given_ind=None, x_given=None, check_sigm
            The mean vector of the multivariate normal distribution.
 
     sigma : array_like
-            Symmetric and positive-definte covariance matrix of
+            Symmetric and positive-definite covariance matrix of
             the multivariate normal distribution.
 
     dependent_ind : array_like
@@ -35,7 +35,7 @@ def condMVN(mean, sigma, dependent_ind, given_ind=None, x_given=None, check_sigm
 
     given_ind : array_like, optional
                 The indices of independent variables (default value is ``None``).
-                If not speficied or all values are zero, return unconditional values.
+                If not specified or all values are zero, return unconditional values.
 
     x_given : array_like, optional
               The conditioning values (default value is ``None``).
@@ -84,7 +84,7 @@ def condMVN(mean, sigma, dependent_ind, given_ind=None, x_given=None, check_sigm
         condVar = np.array(sigma_np[dependent_ind, :][:, dependent_ind])
         return (condMean, condVar)
 
-    # Make sure that x_given is not empty & aligns with given_len:
+    # Make sure that x_given is not empty and aligns with given_len:
     x_given_np = np.array(x_given).squeeze()
     try:
         if len(x_given_np) != len(given_ind_np):
diff --git a/econsa/tests/test_sampling.py b/econsa/tests/test_sampling.py
index bca59b89..58548241 100644
--- a/econsa/tests/test_sampling.py
+++ b/econsa/tests/test_sampling.py
@@ -1,9 +1,75 @@
 """Tests for the sampling.
 
-This module contains all tests for th sampling setup.
+This module contains all tests for the sampling setup.
 
 """
+import numpy as np
+from numpy.random import RandomState
+from rpy2 import robjects
+from rpy2.robjects import numpy2ri
+from rpy2.robjects.packages import importr
+
+from econsa.sampling import condMVN
+
+# Since NumPy v1.17: from numpy.random import default_rng
+
+
+# Import R modules
+base = importr("base")
+stats = importr("stats")
+condMVNorm = importr("condMVNorm")
+
+# Import numpy.random generator
+# Since NumPy v1.17: rng = default_rng()
+rs = RandomState()
 
 
 def test_sampling():
     pass
+
+
+def get_strategies(name):
+    if name == "condMVN":
+        n = rs.randint(low=3, high=20)
+        mean = rs.randint(low=-2, high=2, size=n)
+        sigma = rs.standard_normal(size=(n, n))
+        sigma = sigma @ sigma.T
+        dependent_n = rs.randint(low=1, high=n - 2)
+        dependent = rs.choice(range(0, n), replace=False, size=dependent_n)
+        given_ind = [x for x in range(0, n) if x not in dependent]
+        x_given = rs.randint(low=-2, high=2, size=len(given_ind))
+        strategy = (n, mean, sigma, dependent, given_ind, x_given)
+    else:
+        raise NotImplementedError
+    return strategy
+
+
+def test_condMVN():
+    """Test condMVN against the original from R package condMVNorm.
+    """
+    # Evaluate Python code
+    n, mean, sigma, dependent_ind, given_ind, x_given = get_strategies("condMVN")
+    cond_mean, cond_var = condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+
+    # Evaluate R code
+    numpy2ri.activate()
+    mean_r = robjects.FloatVector(mean)
+    sigma_r = base.matrix(sigma, n, n)
+    dependent_ind_r = robjects.IntVector([x + 1 for x in dependent_ind])
+    given_ind_r = robjects.IntVector([x + 1 for x in given_ind])
+    x_given_r = robjects.IntVector(x_given)
+    condMean_r, condVar_r = condMVNorm.condMVN(
+        mean=mean_r,
+        sigma=sigma_r,
+        dependent=dependent_ind_r,
+        given=given_ind_r,
+        X=x_given_r,
+    )
+    numpy2ri.deactivate()
+
+    condMean_r = np.array(condMean_r)
+    condVar_r = np.array(condVar_r)
+
+    # Comparison
+    np.testing.assert_allclose(cond_mean, condMean_r)
+    np.testing.assert_allclose(cond_var, condVar_r)
diff --git a/environment.yml b/environment.yml
index 38590cf1..fe795fb5 100644
--- a/environment.yml
+++ b/environment.yml
@@ -38,3 +38,4 @@ dependencies:
     - flake8-nb
     - numpydoc
     - pre-commit
+    - rpy2

From e138eeee400a71490ad99bd4a3808d96d6277bea Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Tue, 23 Jun 2020 19:11:04 +0200
Subject: [PATCH 04/45] setup R for CI

---
 .github/workflows/continuous-integration-workflow.yml | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/.github/workflows/continuous-integration-workflow.yml b/.github/workflows/continuous-integration-workflow.yml
index 2be936a8..368517a7 100644
--- a/.github/workflows/continuous-integration-workflow.yml
+++ b/.github/workflows/continuous-integration-workflow.yml
@@ -16,6 +16,10 @@ jobs:
     steps:
     - uses: actions/checkout@v2
 
+    - uses: r-lib/actions/setup-r@v1
+      with:
+        r-version: '3.5.3'
+
     - uses: goanpeca/setup-miniconda@v1.0.2
       with:
           auto-update-conda: true

From a99f286ae2f7ff0d99ca8c55e033f60639848976 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Tue, 23 Jun 2020 19:19:39 +0200
Subject: [PATCH 05/45] use newer version of R

---
 .github/workflows/continuous-integration-workflow.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/continuous-integration-workflow.yml b/.github/workflows/continuous-integration-workflow.yml
index 368517a7..44608936 100644
--- a/.github/workflows/continuous-integration-workflow.yml
+++ b/.github/workflows/continuous-integration-workflow.yml
@@ -18,7 +18,7 @@ jobs:
 
     - uses: r-lib/actions/setup-r@v1
       with:
-        r-version: '3.5.3'
+        r-version: '4.0.1'
 
     - uses: goanpeca/setup-miniconda@v1.0.2
       with:

From ad27a1e95e001642c5a8d62604666b93382a7625 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Tue, 23 Jun 2020 19:35:05 +0200
Subject: [PATCH 06/45] install r packages

---
 .github/workflows/continuous-integration-workflow.yml | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/.github/workflows/continuous-integration-workflow.yml b/.github/workflows/continuous-integration-workflow.yml
index 44608936..d5aff440 100644
--- a/.github/workflows/continuous-integration-workflow.yml
+++ b/.github/workflows/continuous-integration-workflow.yml
@@ -20,6 +20,11 @@ jobs:
       with:
         r-version: '4.0.1'
 
+    - name: Install R dependencies.
+      shell: bash -l {0}
+      run: |
+        Rscript -e 'install.packages("condMVNorm")'
+
     - uses: goanpeca/setup-miniconda@v1.0.2
       with:
           auto-update-conda: true
@@ -30,6 +35,7 @@ jobs:
     - name: Run test suite.
       shell: bash -l {0}
       run: |
+        export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:/usr/lib/R/lib
         pytest --cov=econsa --cov-report=xml
 
         # We want to make sure proper formatting of the notebooks.

From f4cfdf16a637ab70b95543ae25ecd1df62a3d3a2 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Tue, 23 Jun 2020 20:00:59 +0200
Subject: [PATCH 07/45] use integrated r for rpy2; install r package on the fly

---
 .../workflows/continuous-integration-workflow.yml   | 13 +++++--------
 econsa/tests/test_sampling.py                       | 11 +++++++----
 2 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/.github/workflows/continuous-integration-workflow.yml b/.github/workflows/continuous-integration-workflow.yml
index d5aff440..f7c5dc28 100644
--- a/.github/workflows/continuous-integration-workflow.yml
+++ b/.github/workflows/continuous-integration-workflow.yml
@@ -16,16 +16,13 @@ jobs:
     steps:
     - uses: actions/checkout@v2
 
-    - uses: r-lib/actions/setup-r@v1
-      with:
-        r-version: '4.0.1'
-
-    - name: Install R dependencies.
-      shell: bash -l {0}
-      run: |
-        Rscript -e 'install.packages("condMVNorm")'
+    # - uses: r-lib/actions/setup-r@v1
+    #   with:
+    #     r-version: '4.0.1'
 
     - uses: goanpeca/setup-miniconda@v1.0.2
+      env:
+        RPY2_CFFI_MODE: ABI
       with:
           auto-update-conda: true
           activate-environment: econsa
diff --git a/econsa/tests/test_sampling.py b/econsa/tests/test_sampling.py
index 58548241..b7826171 100644
--- a/econsa/tests/test_sampling.py
+++ b/econsa/tests/test_sampling.py
@@ -4,10 +4,10 @@
 
 """
 import numpy as np
+import rpy2.robjects.packages as rpackages
 from numpy.random import RandomState
 from rpy2 import robjects
 from rpy2.robjects import numpy2ri
-from rpy2.robjects.packages import importr
 
 from econsa.sampling import condMVN
 
@@ -15,9 +15,12 @@
 
 
 # Import R modules
-base = importr("base")
-stats = importr("stats")
-condMVNorm = importr("condMVNorm")
+base = rpackages.importr("base")
+stats = rpackages.importr("stats")
+utils = rpackages.importr("utils")
+utils.chooseCRANmirror(ind=1)
+utils.install_packages("condMVNorm")
+condMVNorm = rpackages.importr("condMVNorm")
 
 # Import numpy.random generator
 # Since NumPy v1.17: rng = default_rng()

From 00484947a7ef95bcdf03c811b9daf04b3fd1d814 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Tue, 23 Jun 2020 20:09:14 +0200
Subject: [PATCH 08/45] fix RPY2_CFFI_MODE

---
 .github/workflows/continuous-integration-workflow.yml | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/.github/workflows/continuous-integration-workflow.yml b/.github/workflows/continuous-integration-workflow.yml
index f7c5dc28..3d5a47d3 100644
--- a/.github/workflows/continuous-integration-workflow.yml
+++ b/.github/workflows/continuous-integration-workflow.yml
@@ -16,10 +16,6 @@ jobs:
     steps:
     - uses: actions/checkout@v2
 
-    # - uses: r-lib/actions/setup-r@v1
-    #   with:
-    #     r-version: '4.0.1'
-
     - uses: goanpeca/setup-miniconda@v1.0.2
       env:
         RPY2_CFFI_MODE: ABI
@@ -31,6 +27,8 @@ jobs:
 
     - name: Run test suite.
       shell: bash -l {0}
+      env:
+        RPY2_CFFI_MODE: ABI
       run: |
         export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:/usr/lib/R/lib
         pytest --cov=econsa --cov-report=xml

From bfba26e1fff676db5264fc99e92b61b8add73fd1 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Tue, 23 Jun 2020 20:26:28 +0200
Subject: [PATCH 09/45] install r & use env variable

---
 .github/workflows/continuous-integration-workflow.yml | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/.github/workflows/continuous-integration-workflow.yml b/.github/workflows/continuous-integration-workflow.yml
index 3d5a47d3..b478eb42 100644
--- a/.github/workflows/continuous-integration-workflow.yml
+++ b/.github/workflows/continuous-integration-workflow.yml
@@ -16,9 +16,11 @@ jobs:
     steps:
     - uses: actions/checkout@v2
 
+    - uses: r-lib/actions/setup-R@v1
+      with:
+          r-version: '4.0.1'
+
     - uses: goanpeca/setup-miniconda@v1.0.2
-      env:
-        RPY2_CFFI_MODE: ABI
       with:
           auto-update-conda: true
           activate-environment: econsa
@@ -28,9 +30,8 @@ jobs:
     - name: Run test suite.
       shell: bash -l {0}
       env:
-        RPY2_CFFI_MODE: ABI
+        LD_LIBRARY_PATH: ${LD_LIBRARY_PATH}:/usr/lib/R/lib
       run: |
-        export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:/usr/lib/R/lib
         pytest --cov=econsa --cov-report=xml
 
         # We want to make sure proper formatting of the notebooks.

From e0679ba59a19cbfdc9e3e1eaba24e5361215ac04 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Tue, 23 Jun 2020 20:39:30 +0200
Subject: [PATCH 10/45] fix range is zero error

---
 econsa/tests/test_sampling.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/econsa/tests/test_sampling.py b/econsa/tests/test_sampling.py
index b7826171..98d62d8b 100644
--- a/econsa/tests/test_sampling.py
+++ b/econsa/tests/test_sampling.py
@@ -33,7 +33,7 @@ def test_sampling():
 
 def get_strategies(name):
     if name == "condMVN":
-        n = rs.randint(low=3, high=20)
+        n = rs.randint(low=4, high=20)
         mean = rs.randint(low=-2, high=2, size=n)
         sigma = rs.standard_normal(size=(n, n))
         sigma = sigma @ sigma.T

From 04b7ec6cf8c1a87616c8faa3c80293060e29c4fd Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Tue, 23 Jun 2020 20:41:33 +0200
Subject: [PATCH 11/45] typo

---
 .github/workflows/continuous-integration-workflow.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/continuous-integration-workflow.yml b/.github/workflows/continuous-integration-workflow.yml
index b478eb42..f68dced7 100644
--- a/.github/workflows/continuous-integration-workflow.yml
+++ b/.github/workflows/continuous-integration-workflow.yml
@@ -16,7 +16,7 @@ jobs:
     steps:
     - uses: actions/checkout@v2
 
-    - uses: r-lib/actions/setup-R@v1
+    - uses: r-lib/actions/setup-r@v1
       with:
           r-version: '4.0.1'
 

From 4d4f948d2ae861ec41c90ceea2e207f7465f5d6f Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Tue, 23 Jun 2020 21:53:21 +0200
Subject: [PATCH 12/45] debug CI: try rpy2.situation

---
 .../continuous-integration-workflow.yml       | 26 ++++++++++++++++---
 1 file changed, 22 insertions(+), 4 deletions(-)

diff --git a/.github/workflows/continuous-integration-workflow.yml b/.github/workflows/continuous-integration-workflow.yml
index f68dced7..b82a8b33 100644
--- a/.github/workflows/continuous-integration-workflow.yml
+++ b/.github/workflows/continuous-integration-workflow.yml
@@ -16,7 +16,7 @@ jobs:
     steps:
     - uses: actions/checkout@v2
 
-    - uses: r-lib/actions/setup-r@v1
+    - uses: r-lib/actions/setup-r@master
       with:
           r-version: '4.0.1'
 
@@ -27,11 +27,29 @@ jobs:
           environment-file: environment.yml
           python-version: ${{ matrix.python-version }}
 
-    - name: Run test suite.
+    - name: Run test suite on non-Windows.
+      if: runner.os != 'Windows'
       shell: bash -l {0}
-      env:
-        LD_LIBRARY_PATH: ${LD_LIBRARY_PATH}:/usr/lib/R/lib
       run: |
+        python -m rpy2.situation
+        export LD_LIBRARY_PATH=$(python -m rpy2.situation LD_LIBRARY_PATH):${LD_LIBRARY_PATH}
+
+        pytest --cov=econsa --cov-report=xml
+
+        # We want to make sure proper formatting of the notebooks.
+        flake8-nb; [ $? -eq 0 ] || exit 1
+        black-nb .; [ $? -eq 0 ] || exit 1
+
+        # We want to pull the submodules.
+        # This can be removed once temfpy is live on conda.
+        git submodule update --init --recursive
+
+    - name: Run test suite on Windows.
+      if: runner.os == 'Windows'
+      shell: bash -l {0}
+      run: |
+        python -m rpy2.situation
+
         pytest --cov=econsa --cov-report=xml
 
         # We want to make sure proper formatting of the notebooks.

From 63ef4414fc4b684ee43726cf0ff9e3b402894ff7 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Tue, 23 Jun 2020 22:35:08 +0200
Subject: [PATCH 13/45] use r 3.6.3

---
 .github/workflows/continuous-integration-workflow.yml | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/.github/workflows/continuous-integration-workflow.yml b/.github/workflows/continuous-integration-workflow.yml
index b82a8b33..eb0722c8 100644
--- a/.github/workflows/continuous-integration-workflow.yml
+++ b/.github/workflows/continuous-integration-workflow.yml
@@ -18,7 +18,7 @@ jobs:
 
     - uses: r-lib/actions/setup-r@master
       with:
-          r-version: '4.0.1'
+          r-version: '3.6.3'
 
     - uses: goanpeca/setup-miniconda@v1.0.2
       with:
@@ -31,7 +31,6 @@ jobs:
       if: runner.os != 'Windows'
       shell: bash -l {0}
       run: |
-        python -m rpy2.situation
         export LD_LIBRARY_PATH=$(python -m rpy2.situation LD_LIBRARY_PATH):${LD_LIBRARY_PATH}
 
         pytest --cov=econsa --cov-report=xml

From 83154d255e1d1cdf86ddb7420bb0750f3cdfb504 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Wed, 24 Jun 2020 16:30:37 +0200
Subject: [PATCH 14/45] skip importing rpy2 for autodoc

---
 docs/source/conf.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/docs/source/conf.py b/docs/source/conf.py
index 2e944629..7c33813d 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -217,3 +217,6 @@ def setup(app):
 
 # Use new CDN path of MathJax
 mathjax_path = "https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"
+
+# Skip import packages
+autodoc_mock_imports = ["rpy2"]

From 19b92e13d304f9a0103b6780aa6e50ce1e52ca7f Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Wed, 24 Jun 2020 16:53:23 +0200
Subject: [PATCH 15/45] add rtd_environment

---
 docs/rtd_environment.yml | 38 ++++++++++++++++++++++++++++++++++++++
 docs/source/conf.py      |  3 ---
 2 files changed, 38 insertions(+), 3 deletions(-)
 create mode 100644 docs/rtd_environment.yml

diff --git a/docs/rtd_environment.yml b/docs/rtd_environment.yml
new file mode 100644
index 00000000..55b2d5bb
--- /dev/null
+++ b/docs/rtd_environment.yml
@@ -0,0 +1,38 @@
+channels:
+  - defaults
+  - conda-forge
+  - opensourceeconomics
+
+dependencies:
+  - python=3.7
+  - sphinx
+  - pandas>=0.25
+  - numpy=1.16
+  - pip
+  - scipy
+  - jupyter
+  - jupyterlab
+  - seaborn
+  - numba
+  - pdbpp
+  - chaospy
+  - pytest
+  - pytest-black
+  - pytest-flake8
+  - pytest-cov
+  - pytest-sugar
+  - hypothesis
+  - nbsphinx
+  - sphinx_rtd_theme
+  - sphinxcontrib-bibtex
+
+  - pip:
+    - salib==1.3.8
+    - isort
+    - pytest-stress
+    - pytest-randomly
+    - black
+    - black-nb
+    - flake8-nb
+    - numpydoc
+    - pre-commit
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 7c33813d..2e944629 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -217,6 +217,3 @@ def setup(app):
 
 # Use new CDN path of MathJax
 mathjax_path = "https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"
-
-# Skip import packages
-autodoc_mock_imports = ["rpy2"]

From 365411752a7b7a519566124e19a749756833b443 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Wed, 24 Jun 2020 17:05:07 +0200
Subject: [PATCH 16/45] change readthedocs.yml

---
 docs/rtd_environment.yml | 10 ----------
 readthedocs.yml          |  2 +-
 2 files changed, 1 insertion(+), 11 deletions(-)

diff --git a/docs/rtd_environment.yml b/docs/rtd_environment.yml
index 55b2d5bb..19a3a30a 100644
--- a/docs/rtd_environment.yml
+++ b/docs/rtd_environment.yml
@@ -1,7 +1,6 @@
 channels:
   - defaults
   - conda-forge
-  - opensourceeconomics
 
 dependencies:
   - python=3.7
@@ -14,7 +13,6 @@ dependencies:
   - jupyterlab
   - seaborn
   - numba
-  - pdbpp
   - chaospy
   - pytest
   - pytest-black
@@ -27,12 +25,4 @@ dependencies:
   - sphinxcontrib-bibtex
 
   - pip:
-    - salib==1.3.8
-    - isort
-    - pytest-stress
-    - pytest-randomly
-    - black
-    - black-nb
-    - flake8-nb
     - numpydoc
-    - pre-commit
diff --git a/readthedocs.yml b/readthedocs.yml
index acf5ea5d..beeb52ae 100644
--- a/readthedocs.yml
+++ b/readthedocs.yml
@@ -7,7 +7,7 @@ python:
   version: 3.7
 
 conda:
-  environment: environment.yml
+  environment: docs/rtd_environment.yml
 
 formats:
   - pdf

From d99b9e7410a5a0e9e92b1285f08290fa6d6e3335 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Wed, 24 Jun 2020 19:20:54 +0200
Subject: [PATCH 17/45] test condMVN exceptions

---
 econsa/tests/test_sampling.py | 32 ++++++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)

diff --git a/econsa/tests/test_sampling.py b/econsa/tests/test_sampling.py
index 98d62d8b..561dd0ad 100644
--- a/econsa/tests/test_sampling.py
+++ b/econsa/tests/test_sampling.py
@@ -4,6 +4,7 @@
 
 """
 import numpy as np
+import pytest
 import rpy2.robjects.packages as rpackages
 from numpy.random import RandomState
 from rpy2 import robjects
@@ -42,6 +43,18 @@ def get_strategies(name):
         given_ind = [x for x in range(0, n) if x not in dependent]
         x_given = rs.randint(low=-2, high=2, size=len(given_ind))
         strategy = (n, mean, sigma, dependent, given_ind, x_given)
+    elif name == "condMVN_exception":
+        n = rs.randint(low=4, high=20)
+        mean = rs.randint(low=-2, high=2, size=n)
+        sigma = rs.standard_normal(size=(n, n))
+        sigma = sigma @ sigma.T
+        dependent_n = rs.randint(low=1, high=n - 2)
+        dependent = rs.choice(range(0, n), replace=False, size=dependent_n)
+        given_ind = [x for x in range(0, n) if x not in dependent] if n % 2 == 0 else []
+        x_given = (
+            rs.randint(low=-2, high=2, size=len(given_ind) + 1) if n % 3 == 0 else []
+        )
+        strategy = (n, mean, sigma, dependent, given_ind, x_given)
     else:
         raise NotImplementedError
     return strategy
@@ -76,3 +89,22 @@ def test_condMVN():
     # Comparison
     np.testing.assert_allclose(cond_mean, condMean_r)
     np.testing.assert_allclose(cond_var, condVar_r)
+
+
+def test_condMVN_exception():
+    """Test condMVN raises exceptions when invalid variables are passed.
+    """
+    n, mean, sigma, dependent_ind, given_ind, x_given = get_strategies(
+        "condMVN_exception",
+    )
+    if n % 2 == 0 and n % 3 == 0:
+        with pytest.raises(ValueError):
+            condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+    elif n % 2 == 0:
+        with pytest.raises(ValueError):
+            condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+    elif n % 3 == 0:
+        with pytest.raises(TypeError):
+            condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+    else:
+        condMVN(mean, sigma, dependent_ind, given_ind, x_given)

From 2fe8f4ce07c194599bdfcad1895d6ff51423b47a Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Wed, 24 Jun 2020 19:37:38 +0200
Subject: [PATCH 18/45] trigger diff

---
 .github/workflows/continuous-integration-workflow.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/continuous-integration-workflow.yml b/.github/workflows/continuous-integration-workflow.yml
index eb0722c8..40567206 100644
--- a/.github/workflows/continuous-integration-workflow.yml
+++ b/.github/workflows/continuous-integration-workflow.yml
@@ -18,7 +18,7 @@ jobs:
 
     - uses: r-lib/actions/setup-r@master
       with:
-          r-version: '3.6.3'
+          r-version: '3.6'
 
     - uses: goanpeca/setup-miniconda@v1.0.2
       with:

From b43827a307aded1db881b8792222ceb86f2c2e78 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Wed, 24 Jun 2020 19:54:11 +0200
Subject: [PATCH 19/45] add test bad sigma

---
 econsa/tests/test_sampling.py | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/econsa/tests/test_sampling.py b/econsa/tests/test_sampling.py
index 561dd0ad..601e9bb5 100644
--- a/econsa/tests/test_sampling.py
+++ b/econsa/tests/test_sampling.py
@@ -47,7 +47,7 @@ def get_strategies(name):
         n = rs.randint(low=4, high=20)
         mean = rs.randint(low=-2, high=2, size=n)
         sigma = rs.standard_normal(size=(n, n))
-        sigma = sigma @ sigma.T
+        sigma = sigma @ sigma.T if n % 5 != 0 else sigma
         dependent_n = rs.randint(low=1, high=n - 2)
         dependent = rs.choice(range(0, n), replace=False, size=dependent_n)
         given_ind = [x for x in range(0, n) if x not in dependent] if n % 2 == 0 else []
@@ -97,14 +97,17 @@ def test_condMVN_exception():
     n, mean, sigma, dependent_ind, given_ind, x_given = get_strategies(
         "condMVN_exception",
     )
-    if n % 2 == 0 and n % 3 == 0:
-        with pytest.raises(ValueError):
-            condMVN(mean, sigma, dependent_ind, given_ind, x_given)
-    elif n % 2 == 0:
-        with pytest.raises(ValueError):
+    if n % 2 == 0:
+        with pytest.raises(ValueError) as e:
             condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+        assert str(e.value) == "lengths of x_given and given_ind must be the same"
     elif n % 3 == 0:
-        with pytest.raises(TypeError):
+        with pytest.raises(TypeError) as e:
+            condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+        assert str(e.value) == "len() of unsized object"
+    elif n % 5 == 0:
+        with pytest.raises(ValueError) as e:
             condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+        assert str(e.value) == "sigma is not a symmetric matrix"
     else:
         condMVN(mean, sigma, dependent_ind, given_ind, x_given)

From 989625c27b5e179dcd24b3b20460817c13c753c8 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Wed, 24 Jun 2020 20:02:16 +0200
Subject: [PATCH 20/45] use assert in instead of equal

---
 econsa/tests/test_sampling.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/econsa/tests/test_sampling.py b/econsa/tests/test_sampling.py
index 601e9bb5..c7253658 100644
--- a/econsa/tests/test_sampling.py
+++ b/econsa/tests/test_sampling.py
@@ -100,14 +100,14 @@ def test_condMVN_exception():
     if n % 2 == 0:
         with pytest.raises(ValueError) as e:
             condMVN(mean, sigma, dependent_ind, given_ind, x_given)
-        assert str(e.value) == "lengths of x_given and given_ind must be the same"
+        assert "lengths of x_given and given_ind must be the same" in str(e.value)
     elif n % 3 == 0:
         with pytest.raises(TypeError) as e:
             condMVN(mean, sigma, dependent_ind, given_ind, x_given)
-        assert str(e.value) == "len() of unsized object"
+        assert "len() of unsized object" in str(e.value)
     elif n % 5 == 0:
         with pytest.raises(ValueError) as e:
             condMVN(mean, sigma, dependent_ind, given_ind, x_given)
-        assert str(e.value) == "sigma is not a symmetric matrix"
+        assert "sigma is not a symmetric matrix" in str(e.value)
     else:
         condMVN(mean, sigma, dependent_ind, given_ind, x_given)

From d5d432844c16726e0313ff302b595cf4c1800362 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Wed, 24 Jun 2020 20:18:26 +0200
Subject: [PATCH 21/45] add test non-p-d sigma

---
 econsa/tests/test_sampling.py | 16 +++++++++++-----
 1 file changed, 11 insertions(+), 5 deletions(-)

diff --git a/econsa/tests/test_sampling.py b/econsa/tests/test_sampling.py
index c7253658..9e6b9823 100644
--- a/econsa/tests/test_sampling.py
+++ b/econsa/tests/test_sampling.py
@@ -34,7 +34,7 @@ def test_sampling():
 
 def get_strategies(name):
     if name == "condMVN":
-        n = rs.randint(low=4, high=20)
+        n = rs.randint(low=4, high=21)
         mean = rs.randint(low=-2, high=2, size=n)
         sigma = rs.standard_normal(size=(n, n))
         sigma = sigma @ sigma.T
@@ -44,13 +44,15 @@ def get_strategies(name):
         x_given = rs.randint(low=-2, high=2, size=len(given_ind))
         strategy = (n, mean, sigma, dependent, given_ind, x_given)
     elif name == "condMVN_exception":
-        n = rs.randint(low=4, high=20)
+        n = rs.randint(low=4, high=31)
         mean = rs.randint(low=-2, high=2, size=n)
-        sigma = rs.standard_normal(size=(n, n))
-        sigma = sigma @ sigma.T if n % 5 != 0 else sigma
+        sigma = rs.standard_normal(size=(n, n)) if n % 5 != 0 else np.ones(shape=(n, n))
+        sigma = sigma @ sigma.T if n % 7 != 0 else sigma
         dependent_n = rs.randint(low=1, high=n - 2)
         dependent = rs.choice(range(0, n), replace=False, size=dependent_n)
-        given_ind = [x for x in range(0, n) if x not in dependent] if n % 2 == 0 else []
+        given_ind = (
+            [x for x in range(0, n) if x not in dependent] if n % 2 == 0 else None
+        )
         x_given = (
             rs.randint(low=-2, high=2, size=len(given_ind) + 1) if n % 3 == 0 else []
         )
@@ -106,6 +108,10 @@ def test_condMVN_exception():
             condMVN(mean, sigma, dependent_ind, given_ind, x_given)
         assert "len() of unsized object" in str(e.value)
     elif n % 5 == 0:
+        with pytest.raises(ValueError) as e:
+            condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+        assert "sigma is not positive-definite" in str(e.value)
+    elif n % 7 == 0:
         with pytest.raises(ValueError) as e:
             condMVN(mean, sigma, dependent_ind, given_ind, x_given)
         assert "sigma is not a symmetric matrix" in str(e.value)

From 9607a2b834e1c0e89ec5930cd90085ae4933b1d2 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Wed, 24 Jun 2020 20:37:13 +0200
Subject: [PATCH 22/45] fix lenth error

---
 econsa/tests/test_sampling.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/econsa/tests/test_sampling.py b/econsa/tests/test_sampling.py
index 9e6b9823..4f50dc17 100644
--- a/econsa/tests/test_sampling.py
+++ b/econsa/tests/test_sampling.py
@@ -46,7 +46,7 @@ def get_strategies(name):
     elif name == "condMVN_exception":
         n = rs.randint(low=4, high=31)
         mean = rs.randint(low=-2, high=2, size=n)
-        sigma = rs.standard_normal(size=(n, n)) if n % 5 != 0 else np.ones(shape=(n, n))
+        sigma = rs.standard_normal(size=(n, n)) if n % 5 != 0 else np.diagflat([-1] * n)
         sigma = sigma @ sigma.T if n % 7 != 0 else sigma
         dependent_n = rs.randint(low=1, high=n - 2)
         dependent = rs.choice(range(0, n), replace=False, size=dependent_n)
@@ -54,7 +54,7 @@ def get_strategies(name):
             [x for x in range(0, n) if x not in dependent] if n % 2 == 0 else None
         )
         x_given = (
-            rs.randint(low=-2, high=2, size=len(given_ind) + 1) if n % 3 == 0 else []
+            rs.randint(low=-2, high=2, size=n - dependent_n + 2) if n % 3 == 0 else []
         )
         strategy = (n, mean, sigma, dependent, given_ind, x_given)
     else:
@@ -106,7 +106,7 @@ def test_condMVN_exception():
     elif n % 3 == 0:
         with pytest.raises(TypeError) as e:
             condMVN(mean, sigma, dependent_ind, given_ind, x_given)
-        assert "len() of unsized object" in str(e.value)
+        assert "len()" in str(e.value)
     elif n % 5 == 0:
         with pytest.raises(ValueError) as e:
             condMVN(mean, sigma, dependent_ind, given_ind, x_given)

From 425398db238120b0258f04b14e2bcec9040fb8e0 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Wed, 24 Jun 2020 22:43:09 +0200
Subject: [PATCH 23/45] rewrite test_condMVN_exception

---
 econsa/sampling.py            |  2 +-
 econsa/tests/test_sampling.py | 65 +++++++++++++++++++++++++++--------
 2 files changed, 51 insertions(+), 16 deletions(-)

diff --git a/econsa/sampling.py b/econsa/sampling.py
index f9cdc8c1..107b26df 100644
--- a/econsa/sampling.py
+++ b/econsa/sampling.py
@@ -96,7 +96,7 @@ def condMVN(mean, sigma, dependent_ind, given_ind=None, x_given=None, check_sigm
     if check_sigma:
         if not np.allclose(sigma_np, sigma_np.T):
             raise ValueError("sigma is not a symmetric matrix")
-        elif not np.all(np.linalg.eigvals(sigma) > 0):
+        elif np.all(np.linalg.eigvals(sigma_np) > 0) == 0:
             raise ValueError("sigma is not positive-definite")
 
     b = sigma_np[dependent_ind, :][:, dependent_ind]
diff --git a/econsa/tests/test_sampling.py b/econsa/tests/test_sampling.py
index 4f50dc17..ab3299ba 100644
--- a/econsa/tests/test_sampling.py
+++ b/econsa/tests/test_sampling.py
@@ -43,20 +43,30 @@ def get_strategies(name):
         given_ind = [x for x in range(0, n) if x not in dependent]
         x_given = rs.randint(low=-2, high=2, size=len(given_ind))
         strategy = (n, mean, sigma, dependent, given_ind, x_given)
-    elif name == "condMVN_exception":
+    elif name == "condMVN_exception_given":
         n = rs.randint(low=4, high=31)
         mean = rs.randint(low=-2, high=2, size=n)
-        sigma = rs.standard_normal(size=(n, n)) if n % 5 != 0 else np.diagflat([-1] * n)
-        sigma = sigma @ sigma.T if n % 7 != 0 else sigma
+        sigma = rs.standard_normal(size=(n, n))
+        sigma = sigma @ sigma.T
         dependent_n = rs.randint(low=1, high=n - 2)
         dependent = rs.choice(range(0, n), replace=False, size=dependent_n)
         given_ind = (
-            [x for x in range(0, n) if x not in dependent] if n % 2 == 0 else None
+            [x for x in range(0, n) if x not in dependent] if n % 3 == 0 else None
         )
         x_given = (
-            rs.randint(low=-2, high=2, size=n - dependent_n + 2) if n % 3 == 0 else []
+            rs.randint(low=-2, high=2, size=n - dependent_n + 1) if n % 5 == 0 else None
         )
         strategy = (n, mean, sigma, dependent, given_ind, x_given)
+    elif name == "test_condMVN_exception_sigma":
+        n = rs.randint(low=4, high=21)
+        mean = rs.randint(low=-2, high=2, size=n)
+        sigma = rs.standard_normal(size=(n, n)) if n % 3 == 0 else np.diagflat([-1] * n)
+        sigma = sigma @ sigma.T if n % 5 == 0 else sigma
+        dependent_n = rs.randint(low=1, high=n - 2)
+        dependent = rs.choice(range(0, n), replace=False, size=dependent_n)
+        given_ind = [x for x in range(0, n) if x not in dependent]
+        x_given = rs.randint(low=-2, high=2, size=len(given_ind))
+        strategy = (n, mean, sigma, dependent, given_ind, x_given)
     else:
         raise NotImplementedError
     return strategy
@@ -93,25 +103,50 @@ def test_condMVN():
     np.testing.assert_allclose(cond_var, condVar_r)
 
 
-def test_condMVN_exception():
-    """Test condMVN raises exceptions when invalid variables are passed.
+def test_condMVN_exception_given():
+    """Test condMVN raises exceptions when invalid `given_ind` or `x_given` is passed.
     """
     n, mean, sigma, dependent_ind, given_ind, x_given = get_strategies(
-        "condMVN_exception",
+        "condMVN_exception_given",
     )
-    if n % 2 == 0:
+    print("n = ", n)
+
+    if n % 3 != 0 and n % 5 != 0:
+        # Valid case: both `given_ind` and `x_given` are empty
+        condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+    elif n % 3 != 0:
+        # Valid case: only `given_ind` is empty
+        condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+    elif n % 5 != 0:
+        # Only `x_given` is empty
+        with pytest.raises(TypeError) as e:
+            condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+        assert "len() of unsized object" in str(e.value)
+    else:
+        # `x_given` is not empty, and does not align with `given_ind`
         with pytest.raises(ValueError) as e:
             condMVN(mean, sigma, dependent_ind, given_ind, x_given)
         assert "lengths of x_given and given_ind must be the same" in str(e.value)
-    elif n % 3 == 0:
-        with pytest.raises(TypeError) as e:
-            condMVN(mean, sigma, dependent_ind, given_ind, x_given)
-        assert "len()" in str(e.value)
-    elif n % 5 == 0:
+
+
+def test_condMVN_exception_sigma():
+    """Test condMVN raises exceptions when invalid `sigma` is passed.
+    """
+    n, mean, sigma, dependent_ind, given_ind, x_given = get_strategies(
+        "test_condMVN_exception_sigma",
+    )
+    print("n = ", n)
+
+    if n % 3 != 0 and n % 5 != 0:
+        # `sigma` is negative definite matrix
         with pytest.raises(ValueError) as e:
             condMVN(mean, sigma, dependent_ind, given_ind, x_given)
         assert "sigma is not positive-definite" in str(e.value)
-    elif n % 7 == 0:
+    elif n % 3 != 0:
+        # Valid case
+        condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+    elif n % 5 != 0:
+        # `sigma` is not symmetric
         with pytest.raises(ValueError) as e:
             condMVN(mean, sigma, dependent_ind, given_ind, x_given)
         assert "sigma is not a symmetric matrix" in str(e.value)

From 0a53a2a91b74d363a3fa92eed1b424cbf3cea066 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Wed, 24 Jun 2020 22:51:49 +0200
Subject: [PATCH 24/45] skip B101 for test files

---
 .bandit | 2 ++
 1 file changed, 2 insertions(+)
 create mode 100644 .bandit

diff --git a/.bandit b/.bandit
new file mode 100644
index 00000000..d213b8e8
--- /dev/null
+++ b/.bandit
@@ -0,0 +1,2 @@
+[bandit]
+pytest_filter: test_* #excludes B101 in files matching test_*

From 4fa2f359296e062e2109fe699bbb30a4baaad06a Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Wed, 24 Jun 2020 23:02:11 +0200
Subject: [PATCH 25/45] delete .bandit

---
 .bandit | 2 --
 1 file changed, 2 deletions(-)
 delete mode 100644 .bandit

diff --git a/.bandit b/.bandit
deleted file mode 100644
index d213b8e8..00000000
--- a/.bandit
+++ /dev/null
@@ -1,2 +0,0 @@
-[bandit]
-pytest_filter: test_* #excludes B101 in files matching test_*

From abff465825cae9447231a5bee1e455af89835042 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Thu, 25 Jun 2020 14:43:37 +0200
Subject: [PATCH 26/45] delete test codes

---
 econsa/tests/test_sampling.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/econsa/tests/test_sampling.py b/econsa/tests/test_sampling.py
index ab3299ba..e51849e1 100644
--- a/econsa/tests/test_sampling.py
+++ b/econsa/tests/test_sampling.py
@@ -109,7 +109,6 @@ def test_condMVN_exception_given():
     n, mean, sigma, dependent_ind, given_ind, x_given = get_strategies(
         "condMVN_exception_given",
     )
-    print("n = ", n)
 
     if n % 3 != 0 and n % 5 != 0:
         # Valid case: both `given_ind` and `x_given` are empty
@@ -135,7 +134,6 @@ def test_condMVN_exception_sigma():
     n, mean, sigma, dependent_ind, given_ind, x_given = get_strategies(
         "test_condMVN_exception_sigma",
     )
-    print("n = ", n)
 
     if n % 3 != 0 and n % 5 != 0:
         # `sigma` is negative definite matrix

From bd15a6aeb724e1a87d8b0f94f213654e1bfbbb83 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Thu, 25 Jun 2020 14:47:11 +0200
Subject: [PATCH 27/45] adjust sampling input validation

---
 econsa/sampling.py | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/econsa/sampling.py b/econsa/sampling.py
index 107b26df..7e7d1c99 100644
--- a/econsa/sampling.py
+++ b/econsa/sampling.py
@@ -69,6 +69,15 @@ def condMVN(mean, sigma, dependent_ind, given_ind=None, x_given=None, check_sigm
     #
     mean_np = np.array(mean).squeeze()
     sigma_np = np.array(sigma).squeeze()
+    given_ind_np = np.array(given_ind).squeeze()
+    x_given_np = np.array(x_given).squeeze()
+
+    # Check sigma is symmetric & positive-definite:
+    if check_sigma:
+        if not np.allclose(sigma_np, sigma_np.T):
+            raise ValueError("sigma is not a symmetric matrix")
+        elif np.all(np.linalg.eigvals(sigma_np) > 0) == 0:
+            raise ValueError("sigma is not positive-definite")
 
     if given_ind is None and x_given is None:
         condMean = np.array(mean_np[dependent_ind])
@@ -76,7 +85,6 @@ def condMVN(mean, sigma, dependent_ind, given_ind=None, x_given=None, check_sigm
         return (condMean, condVar)
 
     # Make sure that given_len is not empty:
-    given_ind_np = np.array(given_ind).squeeze()
     try:
         len(given_ind_np)
     except TypeError:
@@ -85,20 +93,12 @@ def condMVN(mean, sigma, dependent_ind, given_ind=None, x_given=None, check_sigm
         return (condMean, condVar)
 
     # Make sure that x_given is not empty and aligns with given_len:
-    x_given_np = np.array(x_given).squeeze()
     try:
         if len(x_given_np) != len(given_ind_np):
             raise ValueError("lengths of x_given and given_ind must be the same")
     except TypeError:
         raise
 
-    # Check sigma is symmetric & positive-definite:
-    if check_sigma:
-        if not np.allclose(sigma_np, sigma_np.T):
-            raise ValueError("sigma is not a symmetric matrix")
-        elif np.all(np.linalg.eigvals(sigma_np) > 0) == 0:
-            raise ValueError("sigma is not positive-definite")
-
     b = sigma_np[dependent_ind, :][:, dependent_ind]
     c = sigma_np[dependent_ind, :][:, given_ind]
     d = sigma_np[given_ind, :][:, given_ind]

From a07e08eb65458cc57dd70458675593ce39ed6d2d Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Thu, 25 Jun 2020 14:58:48 +0200
Subject: [PATCH 28/45] trigger diff

---
 econsa/sampling.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/econsa/sampling.py b/econsa/sampling.py
index 7e7d1c99..4f82660a 100644
--- a/econsa/sampling.py
+++ b/econsa/sampling.py
@@ -59,7 +59,7 @@ def condMVN(mean, sigma, dependent_ind, given_ind=None, x_given=None, check_sigm
     >>> sigma = np.array([[4.0677098, -0.9620331, 0.9897267],
     ...                   [-0.9620331, 2.2775449, 0.7475968],
     ...                   [0.9897267, 0.7475968, 0.7336631]])
-    >>> dependent_ind = [0,]
+    >>> dependent_ind = [0, ]
     >>> given_ind = [1, 2]
     >>> x_given = [1, -1]
     >>> cond_mean, cond_var = condMVN(mean, sigma, dependent_ind, given_ind, x_given)

From b353ef254ebb320e34c977af92fa8a5edac14df8 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Thu, 25 Jun 2020 15:05:16 +0200
Subject: [PATCH 29/45] adjust test_condMVN_exception frequency

---
 econsa/tests/test_sampling.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/econsa/tests/test_sampling.py b/econsa/tests/test_sampling.py
index e51849e1..c6efdf48 100644
--- a/econsa/tests/test_sampling.py
+++ b/econsa/tests/test_sampling.py
@@ -54,14 +54,14 @@ def get_strategies(name):
             [x for x in range(0, n) if x not in dependent] if n % 3 == 0 else None
         )
         x_given = (
-            rs.randint(low=-2, high=2, size=n - dependent_n + 1) if n % 5 == 0 else None
+            rs.randint(low=-2, high=2, size=n - dependent_n + 1) if n % 2 == 0 else None
         )
         strategy = (n, mean, sigma, dependent, given_ind, x_given)
     elif name == "test_condMVN_exception_sigma":
-        n = rs.randint(low=4, high=21)
+        n = rs.randint(low=4, high=31)
         mean = rs.randint(low=-2, high=2, size=n)
         sigma = rs.standard_normal(size=(n, n)) if n % 3 == 0 else np.diagflat([-1] * n)
-        sigma = sigma @ sigma.T if n % 5 == 0 else sigma
+        sigma = sigma @ sigma.T if n % 2 == 0 else sigma
         dependent_n = rs.randint(low=1, high=n - 2)
         dependent = rs.choice(range(0, n), replace=False, size=dependent_n)
         given_ind = [x for x in range(0, n) if x not in dependent]
@@ -110,13 +110,13 @@ def test_condMVN_exception_given():
         "condMVN_exception_given",
     )
 
-    if n % 3 != 0 and n % 5 != 0:
+    if n % 3 != 0 and n % 2 != 0:
         # Valid case: both `given_ind` and `x_given` are empty
         condMVN(mean, sigma, dependent_ind, given_ind, x_given)
     elif n % 3 != 0:
         # Valid case: only `given_ind` is empty
         condMVN(mean, sigma, dependent_ind, given_ind, x_given)
-    elif n % 5 != 0:
+    elif n % 2 != 0:
         # Only `x_given` is empty
         with pytest.raises(TypeError) as e:
             condMVN(mean, sigma, dependent_ind, given_ind, x_given)
@@ -135,7 +135,7 @@ def test_condMVN_exception_sigma():
         "test_condMVN_exception_sigma",
     )
 
-    if n % 3 != 0 and n % 5 != 0:
+    if n % 3 != 0 and n % 2 != 0:
         # `sigma` is negative definite matrix
         with pytest.raises(ValueError) as e:
             condMVN(mean, sigma, dependent_ind, given_ind, x_given)
@@ -143,7 +143,7 @@ def test_condMVN_exception_sigma():
     elif n % 3 != 0:
         # Valid case
         condMVN(mean, sigma, dependent_ind, given_ind, x_given)
-    elif n % 5 != 0:
+    elif n % 2 != 0:
         # `sigma` is not symmetric
         with pytest.raises(ValueError) as e:
             condMVN(mean, sigma, dependent_ind, given_ind, x_given)

From f28453c8f0727e440c89f559d652d9352fa5dba5 Mon Sep 17 00:00:00 2001
From: peisenha <peisenha@protonmail.com>
Date: Wed, 1 Jul 2020 19:02:17 +0200
Subject: [PATCH 30/45] moved r-base to conda

---
 .github/workflows/continuous-integration-workflow.yml | 6 ------
 environment.yml                                       | 2 +-
 2 files changed, 1 insertion(+), 7 deletions(-)

diff --git a/.github/workflows/continuous-integration-workflow.yml b/.github/workflows/continuous-integration-workflow.yml
index 40567206..40919960 100644
--- a/.github/workflows/continuous-integration-workflow.yml
+++ b/.github/workflows/continuous-integration-workflow.yml
@@ -16,10 +16,6 @@ jobs:
     steps:
     - uses: actions/checkout@v2
 
-    - uses: r-lib/actions/setup-r@master
-      with:
-          r-version: '3.6'
-
     - uses: goanpeca/setup-miniconda@v1.0.2
       with:
           auto-update-conda: true
@@ -31,7 +27,6 @@ jobs:
       if: runner.os != 'Windows'
       shell: bash -l {0}
       run: |
-        export LD_LIBRARY_PATH=$(python -m rpy2.situation LD_LIBRARY_PATH):${LD_LIBRARY_PATH}
 
         pytest --cov=econsa --cov-report=xml
 
@@ -48,7 +43,6 @@ jobs:
       shell: bash -l {0}
       run: |
         python -m rpy2.situation
-
         pytest --cov=econsa --cov-report=xml
 
         # We want to make sure proper formatting of the notebooks.
diff --git a/environment.yml b/environment.yml
index fe795fb5..49f252b2 100644
--- a/environment.yml
+++ b/environment.yml
@@ -27,7 +27,7 @@ dependencies:
   - nbsphinx
   - sphinx_rtd_theme
   - sphinxcontrib-bibtex
-
+  - r-base
   - pip:
     - salib==1.3.8
     - isort

From 94b5def030461c0f5e44d81b58b8f6ec1602e8c7 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Wed, 1 Jul 2020 20:58:02 +0200
Subject: [PATCH 31/45] no camel case

---
 econsa/sampling.py            |  43 ++++++-------
 econsa/tests/test_sampling.py | 113 +++++++++++++++++-----------------
 2 files changed, 75 insertions(+), 81 deletions(-)

diff --git a/econsa/sampling.py b/econsa/sampling.py
index 4f82660a..7e561bf7 100644
--- a/econsa/sampling.py
+++ b/econsa/sampling.py
@@ -10,7 +10,9 @@ def get_sample(*args, **kwargs):
     return
 
 
-def condMVN(mean, sigma, dependent_ind, given_ind=None, x_given=None, check_sigma=True):
+def cond_mvn(
+    mean, sigma, dependent_ind, given_ind=None, given_value=None, check_sigma=True,
+):
     r"""Conditional mean and variance function.
 
     Returns conditional mean and variance of dependent variables,
@@ -37,7 +39,7 @@ def condMVN(mean, sigma, dependent_ind, given_ind=None, x_given=None, check_sigm
                 The indices of independent variables (default value is ``None``).
                 If not specified or all values are zero, return unconditional values.
 
-    x_given : array_like, optional
+    given_value : array_like, optional
               The conditioning values (default value is ``None``).
               Should be the same length as `given_ind`, otherwise throw an error.
 
@@ -61,50 +63,41 @@ def condMVN(mean, sigma, dependent_ind, given_ind=None, x_given=None, check_sigm
     ...                   [0.9897267, 0.7475968, 0.7336631]])
     >>> dependent_ind = [0, ]
     >>> given_ind = [1, 2]
-    >>> x_given = [1, -1]
-    >>> cond_mean, cond_var = condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+    >>> given_value = [1, -1]
+    >>> cond_mean, cond_var = cond_mvn(mean, sigma, dependent_ind, given_ind, given_value)
     >>> np.testing.assert_almost_equal(cond_mean, -4.347531, decimal=6)
     >>> np.testing.assert_almost_equal(cond_var, 0.170718, decimal=6)
     """
     #
-    mean_np = np.array(mean).squeeze()
-    sigma_np = np.array(sigma).squeeze()
-    given_ind_np = np.array(given_ind).squeeze()
-    x_given_np = np.array(x_given).squeeze()
+    mean_np = np.array(mean)
+    sigma_np = np.array(sigma)
+    given_ind_np = np.array(given_ind, ndmin=1)
+    given_value_np = np.array(given_value, ndmin=1)
 
-    # Check sigma is symmetric & positive-definite:
+    # Check `sigma` is symmetric & positive-definite:
     if check_sigma:
         if not np.allclose(sigma_np, sigma_np.T):
             raise ValueError("sigma is not a symmetric matrix")
         elif np.all(np.linalg.eigvals(sigma_np) > 0) == 0:
             raise ValueError("sigma is not positive-definite")
 
-    if given_ind is None and x_given is None:
+    # When `given_ind` is empty, return mean and variances of dependent values:
+    if np.all(np.array(None, ndmin=1) == given_ind):
         condMean = np.array(mean_np[dependent_ind])
         condVar = np.array(sigma_np[dependent_ind, :][:, dependent_ind])
         return (condMean, condVar)
 
-    # Make sure that given_len is not empty:
-    try:
-        len(given_ind_np)
-    except TypeError:
-        condMean = np.array(mean_np[dependent_ind])
-        condVar = np.array(sigma_np[dependent_ind, :][:, dependent_ind])
-        return (condMean, condVar)
-
-    # Make sure that x_given is not empty and aligns with given_len:
-    try:
-        if len(x_given_np) != len(given_ind_np):
-            raise ValueError("lengths of x_given and given_ind must be the same")
-    except TypeError:
-        raise
+    # Make sure that `given_value` aligns with `given_len`:
+    # This includes the case that `given_value` is empty.
+    if len(given_value_np) != len(given_ind_np):
+        raise ValueError("lengths of given_value and given_ind must be the same")
 
     b = sigma_np[dependent_ind, :][:, dependent_ind]
     c = sigma_np[dependent_ind, :][:, given_ind]
     d = sigma_np[given_ind, :][:, given_ind]
     c_dinv = c @ np.linalg.inv(d)
 
-    cond_mean = mean_np[dependent_ind] + c_dinv @ (x_given - mean_np[given_ind])
+    cond_mean = mean_np[dependent_ind] + c_dinv @ (given_value - mean_np[given_ind])
     cond_var = b - c_dinv @ (c.T)
 
     return (cond_mean, cond_var)
diff --git a/econsa/tests/test_sampling.py b/econsa/tests/test_sampling.py
index c6efdf48..30adaf20 100644
--- a/econsa/tests/test_sampling.py
+++ b/econsa/tests/test_sampling.py
@@ -10,18 +10,18 @@
 from rpy2 import robjects
 from rpy2.robjects import numpy2ri
 
-from econsa.sampling import condMVN
+from econsa.sampling import cond_mvn
 
 # Since NumPy v1.17: from numpy.random import default_rng
 
 
 # Import R modules
-base = rpackages.importr("base")
-stats = rpackages.importr("stats")
-utils = rpackages.importr("utils")
-utils.chooseCRANmirror(ind=1)
-utils.install_packages("condMVNorm")
-condMVNorm = rpackages.importr("condMVNorm")
+r_base = rpackages.importr("base")
+r_stats = rpackages.importr("stats")
+r_utils = rpackages.importr("utils")
+r_utils.chooseCRANmirror(ind=1)
+r_utils.install_packages("condMVNorm")
+r_cond_mvnorm = rpackages.importr("condMVNorm")
 
 # Import numpy.random generator
 # Since NumPy v1.17: rng = default_rng()
@@ -33,7 +33,7 @@ def test_sampling():
 
 
 def get_strategies(name):
-    if name == "condMVN":
+    if name == "cond_mvn":
         n = rs.randint(low=4, high=21)
         mean = rs.randint(low=-2, high=2, size=n)
         sigma = rs.standard_normal(size=(n, n))
@@ -41,9 +41,9 @@ def get_strategies(name):
         dependent_n = rs.randint(low=1, high=n - 2)
         dependent = rs.choice(range(0, n), replace=False, size=dependent_n)
         given_ind = [x for x in range(0, n) if x not in dependent]
-        x_given = rs.randint(low=-2, high=2, size=len(given_ind))
-        strategy = (n, mean, sigma, dependent, given_ind, x_given)
-    elif name == "condMVN_exception_given":
+        given_value = rs.randint(low=-2, high=2, size=len(given_ind))
+        strategy = (n, mean, sigma, dependent, given_ind, given_value)
+    elif name == "cond_mvn_exception_given":
         n = rs.randint(low=4, high=31)
         mean = rs.randint(low=-2, high=2, size=n)
         sigma = rs.standard_normal(size=(n, n))
@@ -53,11 +53,11 @@ def get_strategies(name):
         given_ind = (
             [x for x in range(0, n) if x not in dependent] if n % 3 == 0 else None
         )
-        x_given = (
+        given_value = (
             rs.randint(low=-2, high=2, size=n - dependent_n + 1) if n % 2 == 0 else None
         )
-        strategy = (n, mean, sigma, dependent, given_ind, x_given)
-    elif name == "test_condMVN_exception_sigma":
+        strategy = (n, mean, sigma, dependent, given_ind, given_value)
+    elif name == "test_cond_mvn_exception_sigma":
         n = rs.randint(low=4, high=31)
         mean = rs.randint(low=-2, high=2, size=n)
         sigma = rs.standard_normal(size=(n, n)) if n % 3 == 0 else np.diagflat([-1] * n)
@@ -65,88 +65,89 @@ def get_strategies(name):
         dependent_n = rs.randint(low=1, high=n - 2)
         dependent = rs.choice(range(0, n), replace=False, size=dependent_n)
         given_ind = [x for x in range(0, n) if x not in dependent]
-        x_given = rs.randint(low=-2, high=2, size=len(given_ind))
-        strategy = (n, mean, sigma, dependent, given_ind, x_given)
+        given_value = rs.randint(low=-2, high=2, size=len(given_ind))
+        strategy = (n, mean, sigma, dependent, given_ind, given_value)
     else:
         raise NotImplementedError
     return strategy
 
 
-def test_condMVN():
-    """Test condMVN against the original from R package condMVNorm.
+def test_cond_mvn():
+    """Test cond_mvn against the original from R package cond_mvnorm.
     """
     # Evaluate Python code
-    n, mean, sigma, dependent_ind, given_ind, x_given = get_strategies("condMVN")
-    cond_mean, cond_var = condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+    n, mean, sigma, dependent_ind, given_ind, given_value = get_strategies("cond_mvn")
+    cond_mean, cond_var = cond_mvn(mean, sigma, dependent_ind, given_ind, given_value)
 
     # Evaluate R code
     numpy2ri.activate()
-    mean_r = robjects.FloatVector(mean)
-    sigma_r = base.matrix(sigma, n, n)
-    dependent_ind_r = robjects.IntVector([x + 1 for x in dependent_ind])
-    given_ind_r = robjects.IntVector([x + 1 for x in given_ind])
-    x_given_r = robjects.IntVector(x_given)
-    condMean_r, condVar_r = condMVNorm.condMVN(
-        mean=mean_r,
-        sigma=sigma_r,
-        dependent=dependent_ind_r,
-        given=given_ind_r,
-        X=x_given_r,
+    r_mean = robjects.FloatVector(mean)
+    r_sigma = r_base.matrix(sigma, n, n)
+    r_dependent_ind = robjects.IntVector([x + 1 for x in dependent_ind])
+    r_given_ind = robjects.IntVector([x + 1 for x in given_ind])
+    r_given_value = robjects.IntVector(given_value)
+    r_cond_mean, r_cond_var = r_cond_mvnorm.condMVN(
+        mean=r_mean,
+        sigma=r_sigma,
+        dependent=r_dependent_ind,
+        given=r_given_ind,
+        X=r_given_value,
     )
-    numpy2ri.deactivate()
 
-    condMean_r = np.array(condMean_r)
-    condVar_r = np.array(condVar_r)
+    r_cond_mean = np.array(r_cond_mean)
+    r_cond_var = np.array(r_cond_var)
+
+    numpy2ri.deactivate()
 
     # Comparison
-    np.testing.assert_allclose(cond_mean, condMean_r)
-    np.testing.assert_allclose(cond_var, condVar_r)
+    np.testing.assert_allclose(cond_mean, r_cond_mean)
+    np.testing.assert_allclose(cond_var, r_cond_var)
 
 
-def test_condMVN_exception_given():
-    """Test condMVN raises exceptions when invalid `given_ind` or `x_given` is passed.
+def test_cond_mvn_exception_given():
+    """Test cond_mvn raises exceptions when invalid `given_ind` or `given_value` is passed.
     """
-    n, mean, sigma, dependent_ind, given_ind, x_given = get_strategies(
-        "condMVN_exception_given",
+    n, mean, sigma, dependent_ind, given_ind, given_value = get_strategies(
+        "cond_mvn_exception_given",
     )
 
     if n % 3 != 0 and n % 2 != 0:
-        # Valid case: both `given_ind` and `x_given` are empty
-        condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+        # Valid case: both `given_ind` and `given_value` are empty
+        cond_mvn(mean, sigma, dependent_ind, given_ind, given_value)
     elif n % 3 != 0:
         # Valid case: only `given_ind` is empty
-        condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+        cond_mvn(mean, sigma, dependent_ind, given_ind, given_value)
     elif n % 2 != 0:
-        # Only `x_given` is empty
+        # Only `given_value` is empty
         with pytest.raises(TypeError) as e:
-            condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+            cond_mvn(mean, sigma, dependent_ind, given_ind, given_value)
         assert "len() of unsized object" in str(e.value)
     else:
-        # `x_given` is not empty, and does not align with `given_ind`
+        # `given_value` is not empty, and does not align with `given_ind`
         with pytest.raises(ValueError) as e:
-            condMVN(mean, sigma, dependent_ind, given_ind, x_given)
-        assert "lengths of x_given and given_ind must be the same" in str(e.value)
+            cond_mvn(mean, sigma, dependent_ind, given_ind, given_value)
+        assert "lengths of given_value and given_ind must be the same" in str(e.value)
 
 
-def test_condMVN_exception_sigma():
-    """Test condMVN raises exceptions when invalid `sigma` is passed.
+def test_cond_mvn_exception_sigma():
+    """Test cond_mvn raises exceptions when invalid `sigma` is passed.
     """
-    n, mean, sigma, dependent_ind, given_ind, x_given = get_strategies(
-        "test_condMVN_exception_sigma",
+    n, mean, sigma, dependent_ind, given_ind, given_value = get_strategies(
+        "test_cond_mvn_exception_sigma",
     )
 
     if n % 3 != 0 and n % 2 != 0:
         # `sigma` is negative definite matrix
         with pytest.raises(ValueError) as e:
-            condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+            cond_mvn(mean, sigma, dependent_ind, given_ind, given_value)
         assert "sigma is not positive-definite" in str(e.value)
     elif n % 3 != 0:
         # Valid case
-        condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+        cond_mvn(mean, sigma, dependent_ind, given_ind, given_value)
     elif n % 2 != 0:
         # `sigma` is not symmetric
         with pytest.raises(ValueError) as e:
-            condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+            cond_mvn(mean, sigma, dependent_ind, given_ind, given_value)
         assert "sigma is not a symmetric matrix" in str(e.value)
     else:
-        condMVN(mean, sigma, dependent_ind, given_ind, x_given)
+        cond_mvn(mean, sigma, dependent_ind, given_ind, given_value)

From 09cbf8625b7a957cd340350104748f85a08e7195 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Wed, 1 Jul 2020 21:03:44 +0200
Subject: [PATCH 32/45] move some strategies out of if

---
 econsa/sampling.py            |  4 ----
 econsa/tests/test_sampling.py | 21 +++++----------------
 2 files changed, 5 insertions(+), 20 deletions(-)

diff --git a/econsa/sampling.py b/econsa/sampling.py
index 7e561bf7..7a7f2339 100644
--- a/econsa/sampling.py
+++ b/econsa/sampling.py
@@ -6,10 +6,6 @@
 import numpy as np
 
 
-def get_sample(*args, **kwargs):
-    return
-
-
 def cond_mvn(
     mean, sigma, dependent_ind, given_ind=None, given_value=None, check_sigma=True,
 ):
diff --git a/econsa/tests/test_sampling.py b/econsa/tests/test_sampling.py
index 30adaf20..b26a0864 100644
--- a/econsa/tests/test_sampling.py
+++ b/econsa/tests/test_sampling.py
@@ -28,28 +28,21 @@
 rs = RandomState()
 
 
-def test_sampling():
-    pass
-
-
 def get_strategies(name):
+    n = rs.randint(low=4, high=21)
+    mean = rs.randint(low=-2, high=2, size=n)
+    dependent_n = rs.randint(low=1, high=n - 2)
+    dependent = rs.choice(range(0, n), replace=False, size=dependent_n)
+
     if name == "cond_mvn":
-        n = rs.randint(low=4, high=21)
-        mean = rs.randint(low=-2, high=2, size=n)
         sigma = rs.standard_normal(size=(n, n))
         sigma = sigma @ sigma.T
-        dependent_n = rs.randint(low=1, high=n - 2)
-        dependent = rs.choice(range(0, n), replace=False, size=dependent_n)
         given_ind = [x for x in range(0, n) if x not in dependent]
         given_value = rs.randint(low=-2, high=2, size=len(given_ind))
         strategy = (n, mean, sigma, dependent, given_ind, given_value)
     elif name == "cond_mvn_exception_given":
-        n = rs.randint(low=4, high=31)
-        mean = rs.randint(low=-2, high=2, size=n)
         sigma = rs.standard_normal(size=(n, n))
         sigma = sigma @ sigma.T
-        dependent_n = rs.randint(low=1, high=n - 2)
-        dependent = rs.choice(range(0, n), replace=False, size=dependent_n)
         given_ind = (
             [x for x in range(0, n) if x not in dependent] if n % 3 == 0 else None
         )
@@ -58,12 +51,8 @@ def get_strategies(name):
         )
         strategy = (n, mean, sigma, dependent, given_ind, given_value)
     elif name == "test_cond_mvn_exception_sigma":
-        n = rs.randint(low=4, high=31)
-        mean = rs.randint(low=-2, high=2, size=n)
         sigma = rs.standard_normal(size=(n, n)) if n % 3 == 0 else np.diagflat([-1] * n)
         sigma = sigma @ sigma.T if n % 2 == 0 else sigma
-        dependent_n = rs.randint(low=1, high=n - 2)
-        dependent = rs.choice(range(0, n), replace=False, size=dependent_n)
         given_ind = [x for x in range(0, n) if x not in dependent]
         given_value = rs.randint(low=-2, high=2, size=len(given_ind))
         strategy = (n, mean, sigma, dependent, given_ind, given_value)

From bea487571921cdadbdcc7b3f2d946882df8a8318 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Wed, 1 Jul 2020 21:47:29 +0200
Subject: [PATCH 33/45] upgrade numpy

---
 docs/source/methods/index.rst    |  1 +
 docs/source/methods/sampling.rst |  5 +++++
 econsa/sampling.py               |  8 +++++---
 econsa/tests/test_sampling.py    | 31 ++++++++++++++++---------------
 environment.yml                  |  2 +-
 5 files changed, 28 insertions(+), 19 deletions(-)
 create mode 100644 docs/source/methods/sampling.rst

diff --git a/docs/source/methods/index.rst b/docs/source/methods/index.rst
index 21819e3d..1eaf337a 100644
--- a/docs/source/methods/index.rst
+++ b/docs/source/methods/index.rst
@@ -9,3 +9,4 @@ Here we explain and document in detail, the methods we implement in the ``econsa
 
    qualitative-analysis
    quantitative-analysis
+   sampling
diff --git a/docs/source/methods/sampling.rst b/docs/source/methods/sampling.rst
new file mode 100644
index 00000000..a4304c7a
--- /dev/null
+++ b/docs/source/methods/sampling.rst
@@ -0,0 +1,5 @@
+Sampling methods
+================
+
+.. automodule:: econsa.sampling
+    :members:
diff --git a/econsa/sampling.py b/econsa/sampling.py
index 7a7f2339..874307cb 100644
--- a/econsa/sampling.py
+++ b/econsa/sampling.py
@@ -14,7 +14,9 @@ def cond_mvn(
     Returns conditional mean and variance of dependent variables,
     given multivariate normal distribution and indices of independent variables.
 
-    This is a translation of the main function of R package "condMVNorm".
+    This is a translation of the main function of R package condMVNorm_.
+
+    .. _condMVNorm: https://cran.r-project.org/web/packages/condMVNorm/index.html
 
     .. math::
         X = (X_{\text{ind}}, X_{\text{dep}}) \sim \mathcal{N}
@@ -46,10 +48,10 @@ def cond_mvn(
     Returns
     -------
     cond_mean : array_like
-               The conditional mean of dependent variables.
+                The conditional mean of dependent variables.
 
     cond_var : arrray_like
-              The conditional covariance matrix of dependent variables.
+               The conditional covariance matrix of dependent variables.
 
     Examples
     --------
diff --git a/econsa/tests/test_sampling.py b/econsa/tests/test_sampling.py
index b26a0864..870bf9d1 100644
--- a/econsa/tests/test_sampling.py
+++ b/econsa/tests/test_sampling.py
@@ -6,14 +6,12 @@
 import numpy as np
 import pytest
 import rpy2.robjects.packages as rpackages
-from numpy.random import RandomState
+from numpy.random import default_rng
 from rpy2 import robjects
 from rpy2.robjects import numpy2ri
 
 from econsa.sampling import cond_mvn
 
-# Since NumPy v1.17: from numpy.random import default_rng
-
 
 # Import R modules
 r_base = rpackages.importr("base")
@@ -24,37 +22,40 @@
 r_cond_mvnorm = rpackages.importr("condMVNorm")
 
 # Import numpy.random generator
-# Since NumPy v1.17: rng = default_rng()
-rs = RandomState()
+rng = default_rng()
 
 
 def get_strategies(name):
-    n = rs.randint(low=4, high=21)
-    mean = rs.randint(low=-2, high=2, size=n)
-    dependent_n = rs.randint(low=1, high=n - 2)
-    dependent = rs.choice(range(0, n), replace=False, size=dependent_n)
+    n = rng.integers(low=4, high=21)
+    mean = rng.integers(low=-2, high=2, size=n)
+    dependent_n = rng.integers(low=1, high=n - 2)
+    dependent = rng.choice(range(0, n), replace=False, size=dependent_n)
 
     if name == "cond_mvn":
-        sigma = rs.standard_normal(size=(n, n))
+        sigma = rng.standard_normal(size=(n, n))
         sigma = sigma @ sigma.T
         given_ind = [x for x in range(0, n) if x not in dependent]
-        given_value = rs.randint(low=-2, high=2, size=len(given_ind))
+        given_value = rng.integers(low=-2, high=2, size=len(given_ind))
         strategy = (n, mean, sigma, dependent, given_ind, given_value)
     elif name == "cond_mvn_exception_given":
-        sigma = rs.standard_normal(size=(n, n))
+        sigma = rng.standard_normal(size=(n, n))
         sigma = sigma @ sigma.T
         given_ind = (
             [x for x in range(0, n) if x not in dependent] if n % 3 == 0 else None
         )
         given_value = (
-            rs.randint(low=-2, high=2, size=n - dependent_n + 1) if n % 2 == 0 else None
+            rng.integers(low=-2, high=2, size=n - dependent_n + 1)
+            if n % 2 == 0
+            else None
         )
         strategy = (n, mean, sigma, dependent, given_ind, given_value)
     elif name == "test_cond_mvn_exception_sigma":
-        sigma = rs.standard_normal(size=(n, n)) if n % 3 == 0 else np.diagflat([-1] * n)
+        sigma = (
+            rng.standard_normal(size=(n, n)) if n % 3 == 0 else np.diagflat([-1] * n)
+        )
         sigma = sigma @ sigma.T if n % 2 == 0 else sigma
         given_ind = [x for x in range(0, n) if x not in dependent]
-        given_value = rs.randint(low=-2, high=2, size=len(given_ind))
+        given_value = rng.integers(low=-2, high=2, size=len(given_ind))
         strategy = (n, mean, sigma, dependent, given_ind, given_value)
     else:
         raise NotImplementedError
diff --git a/environment.yml b/environment.yml
index fe795fb5..5387189b 100644
--- a/environment.yml
+++ b/environment.yml
@@ -9,7 +9,7 @@ dependencies:
   - python=3.7
   - sphinx
   - pandas>=0.25
-  - numpy=1.16
+  - numpy>=1.17
   - pip
   - scipy
   - jupyter

From eee93ec368d787b3b1124df8708c800bde535972 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Wed, 1 Jul 2020 23:00:04 +0200
Subject: [PATCH 34/45] move test cases

---
 docs/rtd_environment.yml      |  2 +-
 econsa/tests/test_sampling.py | 16 +++-------------
 2 files changed, 4 insertions(+), 14 deletions(-)

diff --git a/docs/rtd_environment.yml b/docs/rtd_environment.yml
index 19a3a30a..29141a9e 100644
--- a/docs/rtd_environment.yml
+++ b/docs/rtd_environment.yml
@@ -6,7 +6,7 @@ dependencies:
   - python=3.7
   - sphinx
   - pandas>=0.25
-  - numpy=1.16
+  - numpy
   - pip
   - scipy
   - jupyter
diff --git a/econsa/tests/test_sampling.py b/econsa/tests/test_sampling.py
index 870bf9d1..48c3ce4a 100644
--- a/econsa/tests/test_sampling.py
+++ b/econsa/tests/test_sampling.py
@@ -101,19 +101,12 @@ def test_cond_mvn_exception_given():
         "cond_mvn_exception_given",
     )
 
-    if n % 3 != 0 and n % 2 != 0:
-        # Valid case: both `given_ind` and `given_value` are empty
-        cond_mvn(mean, sigma, dependent_ind, given_ind, given_value)
-    elif n % 3 != 0:
+    if n % 3 != 0:
         # Valid case: only `given_ind` is empty
+        # or both `given_ind` and `given_value` are empty
         cond_mvn(mean, sigma, dependent_ind, given_ind, given_value)
-    elif n % 2 != 0:
-        # Only `given_value` is empty
-        with pytest.raises(TypeError) as e:
-            cond_mvn(mean, sigma, dependent_ind, given_ind, given_value)
-        assert "len() of unsized object" in str(e.value)
     else:
-        # `given_value` is not empty, and does not align with `given_ind`
+        # `given_value` is empty or does not align with `given_ind`
         with pytest.raises(ValueError) as e:
             cond_mvn(mean, sigma, dependent_ind, given_ind, given_value)
         assert "lengths of given_value and given_ind must be the same" in str(e.value)
@@ -131,9 +124,6 @@ def test_cond_mvn_exception_sigma():
         with pytest.raises(ValueError) as e:
             cond_mvn(mean, sigma, dependent_ind, given_ind, given_value)
         assert "sigma is not positive-definite" in str(e.value)
-    elif n % 3 != 0:
-        # Valid case
-        cond_mvn(mean, sigma, dependent_ind, given_ind, given_value)
     elif n % 2 != 0:
         # `sigma` is not symmetric
         with pytest.raises(ValueError) as e:

From 205326b0602e053bcef022be5c7531e6a6e9651a Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Thu, 2 Jul 2020 18:12:27 +0200
Subject: [PATCH 35/45] fix documentation types

---
 docs/source/conf.py | 17 +++++++++++++----
 econsa/sampling.py  | 10 +++++-----
 2 files changed, 18 insertions(+), 9 deletions(-)

diff --git a/docs/source/conf.py b/docs/source/conf.py
index 2e944629..31063a58 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -58,7 +58,7 @@ def setup(app):
     "sphinx.ext.extlinks",
     "sphinx.ext.intersphinx",
     "sphinx.ext.viewcode",
-    "numpydoc",
+    "sphinx.ext.napoleon",
     "nbsphinx",
     "sphinx.ext.mathjax",
     "sphinxcontrib.bibtex",
@@ -92,9 +92,18 @@ def setup(app):
 # The name of the Pygments (syntax highlighting) style to use.
 pygments_style = None
 
-# Configuration for numpydoc
-numpydoc_xref_param_type = True
-numpydoc_xref_ignore = {"type", "optional", "default"}
+# Napoleon settings
+napoleon_google_docstring = True
+napoleon_numpy_docstring = True
+napoleon_include_init_with_doc = False
+napoleon_include_private_with_doc = False
+napoleon_include_special_with_doc = True
+napoleon_use_admonition_for_examples = False
+napoleon_use_admonition_for_notes = False
+napoleon_use_admonition_for_references = False
+napoleon_use_ivar = True
+napoleon_use_param = True
+napoleon_use_rtype = True
 
 # -- Options for HTML output -------------------------------------------------
 
diff --git a/econsa/sampling.py b/econsa/sampling.py
index 874307cb..1f5b3c75 100644
--- a/econsa/sampling.py
+++ b/econsa/sampling.py
@@ -1,6 +1,6 @@
 """Capabilities for sampling of random input parameters.
 
-This module contains all we need to sample random input parameters.
+This module contains functions to sample random input parameters.
 
 """
 import numpy as np
@@ -16,7 +16,7 @@ def cond_mvn(
 
     This is a translation of the main function of R package condMVNorm_.
 
-    .. _condMVNorm: https://cran.r-project.org/web/packages/condMVNorm/index.html
+    .. _condMVNorm: https://cran.r-project.org/package=condMVNorm
 
     .. math::
         X = (X_{\text{ind}}, X_{\text{dep}}) \sim \mathcal{N}
@@ -41,7 +41,7 @@ def cond_mvn(
               The conditioning values (default value is ``None``).
               Should be the same length as `given_ind`, otherwise throw an error.
 
-    check_sigma : boolean, optional
+    check_sigma : bool, optional
                   Check that `sigma` is symmetric,
                   and no eigenvalue is zero (default value is ``True``).
 
@@ -83,7 +83,7 @@ def cond_mvn(
     if np.all(np.array(None, ndmin=1) == given_ind):
         condMean = np.array(mean_np[dependent_ind])
         condVar = np.array(sigma_np[dependent_ind, :][:, dependent_ind])
-        return (condMean, condVar)
+        return condMean, condVar
 
     # Make sure that `given_value` aligns with `given_len`:
     # This includes the case that `given_value` is empty.
@@ -98,4 +98,4 @@ def cond_mvn(
     cond_mean = mean_np[dependent_ind] + c_dinv @ (given_value - mean_np[given_ind])
     cond_var = b - c_dinv @ (c.T)
 
-    return (cond_mean, cond_var)
+    return cond_mean, cond_var

From 195edc1a12cf4fcf9a026bceacd179a7e0190cea Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Thu, 2 Jul 2020 18:24:16 +0200
Subject: [PATCH 36/45] add some math

---
 econsa/sampling.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/econsa/sampling.py b/econsa/sampling.py
index 1f5b3c75..88f5d6c6 100644
--- a/econsa/sampling.py
+++ b/econsa/sampling.py
@@ -19,7 +19,9 @@ def cond_mvn(
     .. _condMVNorm: https://cran.r-project.org/package=condMVNorm
 
     .. math::
-        X = (X_{\text{ind}}, X_{\text{dep}}) \sim \mathcal{N}
+        X = (X_{\text{given}}, X_{\text{dep}}) \sim \mathcal{N}
+
+    where the average of :math:`X` is `mean` and the covariance matrix of :math:`X` is `sigma`.
 
     Parameters
     ----------

From 6f0ecda6180cfbe8e819f38167e6cbe87378221a Mon Sep 17 00:00:00 2001
From: peisenha <peisenha@protonmail.com>
Date: Fri, 3 Jul 2020 10:50:37 +0200
Subject: [PATCH 37/45] refactoring R wrapper

---
 econsa/sampling.py            | 46 ++++++++++++---------------
 econsa/tests/test_sampling.py | 60 +++++++++--------------------------
 econsa/tests/wrapper_r.py     | 32 +++++++++++++++++++
 environment.yml               |  2 ++
 4 files changed, 70 insertions(+), 70 deletions(-)
 create mode 100644 econsa/tests/wrapper_r.py

diff --git a/econsa/sampling.py b/econsa/sampling.py
index 88f5d6c6..02b8a2a9 100644
--- a/econsa/sampling.py
+++ b/econsa/sampling.py
@@ -11,17 +11,14 @@ def cond_mvn(
 ):
     r"""Conditional mean and variance function.
 
-    Returns conditional mean and variance of dependent variables,
-    given multivariate normal distribution and indices of independent variables.
+    This function provides the conditional mean and variance-covariance matrix of [Y given X],
+    where Z = (X,Y) is the fully-joint multivariate normal distribution with mean equal to
+    ``mean`` and covariance matrix ``sigma``.
 
     This is a translation of the main function of R package condMVNorm_.
 
     .. _condMVNorm: https://cran.r-project.org/package=condMVNorm
 
-    .. math::
-        X = (X_{\text{given}}, X_{\text{dep}}) \sim \mathcal{N}
-
-    where the average of :math:`X` is `mean` and the covariance matrix of :math:`X` is `sigma`.
 
     Parameters
     ----------
@@ -37,22 +34,22 @@ def cond_mvn(
 
     given_ind : array_like, optional
                 The indices of independent variables (default value is ``None``).
-                If not specified or all values are zero, return unconditional values.
+                If not specified return unconditional values.
 
     given_value : array_like, optional
-              The conditioning values (default value is ``None``).
-              Should be the same length as `given_ind`, otherwise throw an error.
+              The conditioning values (default value is ``None``). Should be the same length as
+              ``given_ind``, otherwise throw an error.
 
     check_sigma : bool, optional
-                  Check that `sigma` is symmetric,
-                  and no eigenvalue is zero (default value is ``True``).
+                  Check that ``sigma`` is symmetric, and no eigenvalue is zero (default value is
+                  ``True``).
 
     Returns
     -------
     cond_mean : array_like
                 The conditional mean of dependent variables.
 
-    cond_var : arrray_like
+    cond_cov : arrray_like
                The conditional covariance matrix of dependent variables.
 
     Examples
@@ -64,31 +61,30 @@ def cond_mvn(
     >>> dependent_ind = [0, ]
     >>> given_ind = [1, 2]
     >>> given_value = [1, -1]
-    >>> cond_mean, cond_var = cond_mvn(mean, sigma, dependent_ind, given_ind, given_value)
+    >>> cond_mean, cond_cov = cond_mvn(mean, sigma, dependent_ind, given_ind, given_value)
     >>> np.testing.assert_almost_equal(cond_mean, -4.347531, decimal=6)
-    >>> np.testing.assert_almost_equal(cond_var, 0.170718, decimal=6)
+    >>> np.testing.assert_almost_equal(cond_cov, 0.170718, decimal=6)
     """
-    #
     mean_np = np.array(mean)
     sigma_np = np.array(sigma)
     given_ind_np = np.array(given_ind, ndmin=1)
     given_value_np = np.array(given_value, ndmin=1)
 
-    # Check `sigma` is symmetric & positive-definite:
+    # Check `sigma` is symmetric and positive-definite:
     if check_sigma:
         if not np.allclose(sigma_np, sigma_np.T):
             raise ValueError("sigma is not a symmetric matrix")
         elif np.all(np.linalg.eigvals(sigma_np) > 0) == 0:
             raise ValueError("sigma is not positive-definite")
 
-    # When `given_ind` is empty, return mean and variances of dependent values:
-    if np.all(np.array(None, ndmin=1) == given_ind):
-        condMean = np.array(mean_np[dependent_ind])
-        condVar = np.array(sigma_np[dependent_ind, :][:, dependent_ind])
-        return condMean, condVar
+    # When `given_ind` is None, return mean and variances of dependent values:
+    if given_ind is None:
+        cond_mean = np.array(mean_np[dependent_ind])
+        cond_cov = np.array(sigma_np[dependent_ind, :][:, dependent_ind])
+        return cond_mean, cond_cov
 
-    # Make sure that `given_value` aligns with `given_len`:
-    # This includes the case that `given_value` is empty.
+    # Make sure that `given_value` aligns with `given_len`. This includes the case that
+    # `given_value` is empty.
     if len(given_value_np) != len(given_ind_np):
         raise ValueError("lengths of given_value and given_ind must be the same")
 
@@ -98,6 +94,6 @@ def cond_mvn(
     c_dinv = c @ np.linalg.inv(d)
 
     cond_mean = mean_np[dependent_ind] + c_dinv @ (given_value - mean_np[given_ind])
-    cond_var = b - c_dinv @ (c.T)
+    cond_cov = b - c_dinv @ (c.T)
 
-    return cond_mean, cond_var
+    return cond_mean, cond_cov
diff --git a/econsa/tests/test_sampling.py b/econsa/tests/test_sampling.py
index 48c3ce4a..71d1359f 100644
--- a/econsa/tests/test_sampling.py
+++ b/econsa/tests/test_sampling.py
@@ -3,25 +3,13 @@
 This module contains all tests for the sampling setup.
 
 """
+from numpy.random import default_rng
 import numpy as np
 import pytest
-import rpy2.robjects.packages as rpackages
-from numpy.random import default_rng
-from rpy2 import robjects
-from rpy2.robjects import numpy2ri
 
+from econsa.tests.wrapper_r import r_cond_mvn
 from econsa.sampling import cond_mvn
 
-
-# Import R modules
-r_base = rpackages.importr("base")
-r_stats = rpackages.importr("stats")
-r_utils = rpackages.importr("utils")
-r_utils.chooseCRANmirror(ind=1)
-r_utils.install_packages("condMVNorm")
-r_cond_mvnorm = rpackages.importr("condMVNorm")
-
-# Import numpy.random generator
 rng = default_rng()
 
 
@@ -36,7 +24,6 @@ def get_strategies(name):
         sigma = sigma @ sigma.T
         given_ind = [x for x in range(0, n) if x not in dependent]
         given_value = rng.integers(low=-2, high=2, size=len(given_ind))
-        strategy = (n, mean, sigma, dependent, given_ind, given_value)
     elif name == "cond_mvn_exception_given":
         sigma = rng.standard_normal(size=(n, n))
         sigma = sigma @ sigma.T
@@ -48,7 +35,6 @@ def get_strategies(name):
             if n % 2 == 0
             else None
         )
-        strategy = (n, mean, sigma, dependent, given_ind, given_value)
     elif name == "test_cond_mvn_exception_sigma":
         sigma = (
             rng.standard_normal(size=(n, n)) if n % 3 == 0 else np.diagflat([-1] * n)
@@ -56,54 +42,36 @@ def get_strategies(name):
         sigma = sigma @ sigma.T if n % 2 == 0 else sigma
         given_ind = [x for x in range(0, n) if x not in dependent]
         given_value = rng.integers(low=-2, high=2, size=len(given_ind))
-        strategy = (n, mean, sigma, dependent, given_ind, given_value)
     else:
         raise NotImplementedError
+
+    strategy = (mean, sigma, dependent, given_ind, given_value)
+
     return strategy
 
 
 def test_cond_mvn():
     """Test cond_mvn against the original from R package cond_mvnorm.
     """
-    # Evaluate Python code
-    n, mean, sigma, dependent_ind, given_ind, given_value = get_strategies("cond_mvn")
-    cond_mean, cond_var = cond_mvn(mean, sigma, dependent_ind, given_ind, given_value)
-
-    # Evaluate R code
-    numpy2ri.activate()
-    r_mean = robjects.FloatVector(mean)
-    r_sigma = r_base.matrix(sigma, n, n)
-    r_dependent_ind = robjects.IntVector([x + 1 for x in dependent_ind])
-    r_given_ind = robjects.IntVector([x + 1 for x in given_ind])
-    r_given_value = robjects.IntVector(given_value)
-    r_cond_mean, r_cond_var = r_cond_mvnorm.condMVN(
-        mean=r_mean,
-        sigma=r_sigma,
-        dependent=r_dependent_ind,
-        given=r_given_ind,
-        X=r_given_value,
-    )
+    request = get_strategies("cond_mvn")
 
-    r_cond_mean = np.array(r_cond_mean)
-    r_cond_var = np.array(r_cond_var)
+    r_cond_mean, r_cond_cov = r_cond_mvn(*request)
+    cond_mean, cond_cov = cond_mvn(*request)
 
-    numpy2ri.deactivate()
-
-    # Comparison
     np.testing.assert_allclose(cond_mean, r_cond_mean)
-    np.testing.assert_allclose(cond_var, r_cond_var)
+    np.testing.assert_allclose(cond_cov, r_cond_cov)
 
 
 def test_cond_mvn_exception_given():
     """Test cond_mvn raises exceptions when invalid `given_ind` or `given_value` is passed.
     """
-    n, mean, sigma, dependent_ind, given_ind, given_value = get_strategies(
+    mean, sigma, dependent_ind, given_ind, given_value = get_strategies(
         "cond_mvn_exception_given",
     )
 
+    n = sigma.shape[0]
     if n % 3 != 0:
-        # Valid case: only `given_ind` is empty
-        # or both `given_ind` and `given_value` are empty
+        # Valid case: only `given_ind` is empty or both `given_ind` and `given_value` are empty
         cond_mvn(mean, sigma, dependent_ind, given_ind, given_value)
     else:
         # `given_value` is empty or does not align with `given_ind`
@@ -115,10 +83,12 @@ def test_cond_mvn_exception_given():
 def test_cond_mvn_exception_sigma():
     """Test cond_mvn raises exceptions when invalid `sigma` is passed.
     """
-    n, mean, sigma, dependent_ind, given_ind, given_value = get_strategies(
+    mean, sigma, dependent_ind, given_ind, given_value = get_strategies(
         "test_cond_mvn_exception_sigma",
     )
 
+    n = sigma.shape[0]
+
     if n % 3 != 0 and n % 2 != 0:
         # `sigma` is negative definite matrix
         with pytest.raises(ValueError) as e:
diff --git a/econsa/tests/wrapper_r.py b/econsa/tests/wrapper_r.py
new file mode 100644
index 00000000..4f415d18
--- /dev/null
+++ b/econsa/tests/wrapper_r.py
@@ -0,0 +1,32 @@
+"""Wrapping R
+
+This module contains all functionality related to the use of functions from R for testing purposes.
+
+"""
+import numpy as np
+
+import rpy2.robjects.packages as rpackages
+from rpy2.robjects import numpy2ri
+from rpy2 import robjects
+
+r_package_cond_mvnorm = rpackages.importr("condMVNorm")
+r_package_base = rpackages.importr("base")
+
+
+def r_cond_mvn(mean, sigma, dependent_ind, given_ind, given_value):
+    numpy2ri.activate()
+    r_mean = robjects.FloatVector(mean)
+    n = sigma.shape[0]
+    r_sigma = r_package_base.matrix(sigma, n, n)
+    r_dependent_ind = robjects.IntVector([x + 1 for x in dependent_ind])
+    r_given_ind = robjects.IntVector([x + 1 for x in given_ind])
+    r_given_value = robjects.IntVector(given_value)
+
+    args = (r_mean, r_sigma, r_dependent_ind, r_given_ind, r_given_value)
+    r_cond_mean, r_cond_cov = r_package_cond_mvnorm.condMVN(*args)
+
+    r_cond_mean, r_cond_cov = np.array(r_cond_mean), np.array(r_cond_cov)
+
+    numpy2ri.deactivate()
+
+    return r_cond_mean, r_cond_cov
diff --git a/environment.yml b/environment.yml
index 077553da..7561a3c6 100644
--- a/environment.yml
+++ b/environment.yml
@@ -28,6 +28,8 @@ dependencies:
   - sphinx_rtd_theme
   - sphinxcontrib-bibtex
   - r-base
+  - r-r.utils
+  - r-condmvnorm
   - pip:
     - salib==1.3.8
     - isort

From 109ff0903be78d484ee137cdf44c9f157dd6473a Mon Sep 17 00:00:00 2001
From: peisenha <peisenha@protonmail.com>
Date: Fri, 3 Jul 2020 11:08:35 +0200
Subject: [PATCH 38/45] further cleanup

---
 econsa/tests/wrapper_r.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/econsa/tests/wrapper_r.py b/econsa/tests/wrapper_r.py
index 4f415d18..21ced36a 100644
--- a/econsa/tests/wrapper_r.py
+++ b/econsa/tests/wrapper_r.py
@@ -10,14 +10,13 @@
 from rpy2 import robjects
 
 r_package_cond_mvnorm = rpackages.importr("condMVNorm")
-r_package_base = rpackages.importr("base")
 
 
 def r_cond_mvn(mean, sigma, dependent_ind, given_ind, given_value):
     numpy2ri.activate()
     r_mean = robjects.FloatVector(mean)
     n = sigma.shape[0]
-    r_sigma = r_package_base.matrix(sigma, n, n)
+    r_sigma = robjects.r.matrix(sigma, n, n)
     r_dependent_ind = robjects.IntVector([x + 1 for x in dependent_ind])
     r_given_ind = robjects.IntVector([x + 1 for x in given_ind])
     r_given_value = robjects.IntVector(given_value)

From 3404720d32e19c9c3f89a10d3e33a8bc1a0b9386 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Fri, 3 Jul 2020 17:04:03 +0200
Subject: [PATCH 39/45] fix pre-commit

---
 econsa/tests/test_sampling.py | 4 ++--
 econsa/tests/wrapper_r.py     | 3 +--
 2 files changed, 3 insertions(+), 4 deletions(-)

diff --git a/econsa/tests/test_sampling.py b/econsa/tests/test_sampling.py
index 71d1359f..16775a53 100644
--- a/econsa/tests/test_sampling.py
+++ b/econsa/tests/test_sampling.py
@@ -3,12 +3,12 @@
 This module contains all tests for the sampling setup.
 
 """
-from numpy.random import default_rng
 import numpy as np
 import pytest
+from numpy.random import default_rng
 
-from econsa.tests.wrapper_r import r_cond_mvn
 from econsa.sampling import cond_mvn
+from econsa.tests.wrapper_r import r_cond_mvn
 
 rng = default_rng()
 
diff --git a/econsa/tests/wrapper_r.py b/econsa/tests/wrapper_r.py
index 21ced36a..4e3492c1 100644
--- a/econsa/tests/wrapper_r.py
+++ b/econsa/tests/wrapper_r.py
@@ -4,10 +4,9 @@
 
 """
 import numpy as np
-
 import rpy2.robjects.packages as rpackages
-from rpy2.robjects import numpy2ri
 from rpy2 import robjects
+from rpy2.robjects import numpy2ri
 
 r_package_cond_mvnorm = rpackages.importr("condMVNorm")
 

From 3ba48800dc158e016d541e0e96db8e113e8a1a9b Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Sat, 4 Jul 2020 16:49:35 +0200
Subject: [PATCH 40/45] test rhome on windows

---
 .github/workflows/continuous-integration-workflow.yml | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/.github/workflows/continuous-integration-workflow.yml b/.github/workflows/continuous-integration-workflow.yml
index 40919960..1c35b3e3 100644
--- a/.github/workflows/continuous-integration-workflow.yml
+++ b/.github/workflows/continuous-integration-workflow.yml
@@ -42,6 +42,8 @@ jobs:
       if: runner.os == 'Windows'
       shell: bash -l {0}
       run: |
+        R RHOME
+        which R
         python -m rpy2.situation
         pytest --cov=econsa --cov-report=xml
 

From 10edee3580a8e796319a2dbb8059ce05419a2e07 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Sat, 4 Jul 2020 16:59:46 +0200
Subject: [PATCH 41/45] add rhome to windows

---
 .github/workflows/continuous-integration-workflow.yml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/continuous-integration-workflow.yml b/.github/workflows/continuous-integration-workflow.yml
index 1c35b3e3..d2694d91 100644
--- a/.github/workflows/continuous-integration-workflow.yml
+++ b/.github/workflows/continuous-integration-workflow.yml
@@ -42,8 +42,8 @@ jobs:
       if: runner.os == 'Windows'
       shell: bash -l {0}
       run: |
-        R RHOME
-        which R
+        # Set R_HOME for rpy2
+        setx R_HOME "C:\Miniconda\envs\econsa\lib\R"
         python -m rpy2.situation
         pytest --cov=econsa --cov-report=xml
 

From 23efdbca70e4beafea9e786fe97c063fc5125049 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Sat, 4 Jul 2020 17:08:29 +0200
Subject: [PATCH 42/45] refresh env

---
 .github/workflows/continuous-integration-workflow.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.github/workflows/continuous-integration-workflow.yml b/.github/workflows/continuous-integration-workflow.yml
index d2694d91..14114241 100644
--- a/.github/workflows/continuous-integration-workflow.yml
+++ b/.github/workflows/continuous-integration-workflow.yml
@@ -44,6 +44,7 @@ jobs:
       run: |
         # Set R_HOME for rpy2
         setx R_HOME "C:\Miniconda\envs\econsa\lib\R"
+        refreshenv
         python -m rpy2.situation
         pytest --cov=econsa --cov-report=xml
 

From 24caa2670574e3d24b853ed1902108aaaaadc5bd Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Sat, 4 Jul 2020 17:23:35 +0200
Subject: [PATCH 43/45] set & setx

---
 .github/workflows/continuous-integration-workflow.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/continuous-integration-workflow.yml b/.github/workflows/continuous-integration-workflow.yml
index 14114241..baaa54ad 100644
--- a/.github/workflows/continuous-integration-workflow.yml
+++ b/.github/workflows/continuous-integration-workflow.yml
@@ -43,8 +43,8 @@ jobs:
       shell: bash -l {0}
       run: |
         # Set R_HOME for rpy2
+        set R_HOME=C:\Miniconda\envs\econsa\lib\R
         setx R_HOME "C:\Miniconda\envs\econsa\lib\R"
-        refreshenv
         python -m rpy2.situation
         pytest --cov=econsa --cov-report=xml
 

From 25f7bfe5fda7cf0e8710800fcdd95828f9e10675 Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Sat, 4 Jul 2020 17:38:46 +0200
Subject: [PATCH 44/45] try https://bitbucket.org/rpy2/rpy2/issues/487/

---
 .github/workflows/continuous-integration-workflow.yml | 3 ---
 econsa/tests/wrapper_r.py                             | 8 +++++++-
 2 files changed, 7 insertions(+), 4 deletions(-)

diff --git a/.github/workflows/continuous-integration-workflow.yml b/.github/workflows/continuous-integration-workflow.yml
index baaa54ad..40919960 100644
--- a/.github/workflows/continuous-integration-workflow.yml
+++ b/.github/workflows/continuous-integration-workflow.yml
@@ -42,9 +42,6 @@ jobs:
       if: runner.os == 'Windows'
       shell: bash -l {0}
       run: |
-        # Set R_HOME for rpy2
-        set R_HOME=C:\Miniconda\envs\econsa\lib\R
-        setx R_HOME "C:\Miniconda\envs\econsa\lib\R"
         python -m rpy2.situation
         pytest --cov=econsa --cov-report=xml
 
diff --git a/econsa/tests/wrapper_r.py b/econsa/tests/wrapper_r.py
index 4e3492c1..fdfc0ae0 100644
--- a/econsa/tests/wrapper_r.py
+++ b/econsa/tests/wrapper_r.py
@@ -1,9 +1,15 @@
-"""Wrapping R
+"""Wrapping R.
 
 This module contains all functionality related to the use of functions from R for testing purposes.
 
 """
+import os
+
 import numpy as np
+
+if os.name == "nt":
+    os.environ["R_HOME"] = "C:/Miniconda/envs/econsa/lib/R"
+
 import rpy2.robjects.packages as rpackages
 from rpy2 import robjects
 from rpy2.robjects import numpy2ri

From 323a6472633479b60ef4337251b2fddb87ca75be Mon Sep 17 00:00:00 2001
From: loikein <wanleiqiong@gmail.com>
Date: Sat, 4 Jul 2020 17:47:18 +0200
Subject: [PATCH 45/45] rollback

---
 econsa/tests/wrapper_r.py | 6 ------
 1 file changed, 6 deletions(-)

diff --git a/econsa/tests/wrapper_r.py b/econsa/tests/wrapper_r.py
index fdfc0ae0..9b2b4f14 100644
--- a/econsa/tests/wrapper_r.py
+++ b/econsa/tests/wrapper_r.py
@@ -3,13 +3,7 @@
 This module contains all functionality related to the use of functions from R for testing purposes.
 
 """
-import os
-
 import numpy as np
-
-if os.name == "nt":
-    os.environ["R_HOME"] = "C:/Miniconda/envs/econsa/lib/R"
-
 import rpy2.robjects.packages as rpackages
 from rpy2 import robjects
 from rpy2.robjects import numpy2ri