Repository: zhengp0/limetr (language: python)

Issue:
add more ipopt options



Discussion / hints:

Fix patch:
diff --git a/.github/workflows/python-build.yml b/.github/workflows/python-build.yml
new file mode 100644
index 0000000..855040b
--- /dev/null
+++ b/.github/workflows/python-build.yml
@@ -0,0 +1,29 @@
+name: python-build
+on: [push]
+jobs:
+  build:
+
+    runs-on: ubuntu-latest
+
+    steps:
+    - uses: actions/checkout@v2
+    - name: Set up Python 3.8
+      uses: actions/setup-python@v2
+      with:
+        python-version: 3.8
+    - name: Install dependencies
+      run: python -m pip install .[dev] --upgrade pip
+    - name: Test with pytest
+      run: pytest
+    - name: Build package distribution
+      if: startsWith(github.ref, 'refs/tags')
+      run: |
+        python -m pip install build 
+        python -m build --sdist --wheel --outdir dist/ .
+    - name: Publish package distribution to PyPI
+      if: startsWith(github.ref, 'refs/tags')
+      uses: pypa/gh-action-pypi-publish@master
+      with:
+        skip_existing: true
+        user: __token__
+        password: ${{ secrets.PYPI_API_TOKEN }}
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..09aae16
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,138 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+pip-wheel-metadata/
+share/python-wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.nox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+*.py,cover
+.hypothesis/
+.pytest_cache/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+db.sqlite3-journal
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# IPython
+profile_default/
+ipython_config.py
+
+# pyenv
+.python-version
+
+# pipenv
+#   According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
+#   However, in case of collaboration, if having platform-specific dependencies or dependencies
+#   having no cross-platform support, pipenv may install dependencies that don't work, or not
+#   install all needed dependencies.
+#Pipfile.lock
+
+# PEP 582; used by e.g. github.com/David-OConnor/pyflow
+__pypackages__/
+
+# Celery stuff
+celerybeat-schedule
+celerybeat.pid
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+.dmypy.json
+dmypy.json
+
+# Pyre type checker
+.pyre/
+
+# pycharm
+.idea/
+
+# vscode
+.vscode/
+
+# examples
+examples/
\ No newline at end of file
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..199e4f8
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2020 Peng Zheng
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
\ No newline at end of file
diff --git a/Makefile b/Makefile
deleted file mode 100644
index 899b463..0000000
--- a/Makefile
+++ /dev/null
@@ -1,5 +0,0 @@
-# make file for pynlme class
-
-clean:
-	find . -name "__pycache__" -exec rm -rf "{}" \;
-	find . -name "*.so" -delete
\ No newline at end of file
diff --git a/README.rst b/README.rst
new file mode 100644
index 0000000..a10c307
--- /dev/null
+++ b/README.rst
@@ -0,0 +1,19 @@
+.. image:: https://github.com/zhengp0/limetr/workflows/python-build/badge.svg
+    :target: https://github.com/zhengp0/limetr/actions
+
+LimeTr
+======
+
+:code:`limetr`, read as "lime tree", is a Python package solves mixed effects models with linear random effects.
+Moreover the package also provides robust regression option called "trimming".
+
+Install
+-------
+
+The package can be installed through ``pip``,
+.. code-block:: bash
+    pip install limetr
+
+Or you could install from source, by clone this repository and do,
+.. code-block:: bash
+    python setup.py install
diff --git a/pynlme/Makefile b/pynlme/Makefile
deleted file mode 100644
index 1fd25ac..0000000
--- a/pynlme/Makefile
+++ /dev/null
@@ -1,4 +0,0 @@
-# make file for the fortran module
-
-futils: utils.f90
-	f2py3 -c -m futils utils.f90 -L/usr/local/lib -llapack -lblas
\ No newline at end of file
diff --git a/pynlme/futils.cpython-36m-x86_64-linux-gnu.so b/pynlme/futils.cpython-36m-x86_64-linux-gnu.so
deleted file mode 100755
index 001d94d..0000000
Binary files a/pynlme/futils.cpython-36m-x86_64-linux-gnu.so and /dev/null differ
diff --git a/pynlme/utils.f90 b/pynlme/utils.f90
deleted file mode 100644
index 1903a73..0000000
--- a/pynlme/utils.f90
+++ /dev/null
@@ -1,336 +0,0 @@
-! fortran helper functions for speed up
-
-MODULE varmat
-
-CONTAINS
-
-! public functions
-! -----------------------------------------------------------------------------
-FUNCTION dot_mv(v, Z, d, study_sizes, x, m, n, k) RESULT(y)
-    IMPLICIT NONE
-    INTEGER(8), INTENT(IN) :: m, n, k
-    REAL(8), INTENT(IN) :: v(n), Z(n, k), d(k), x(n)
-    INTEGER(8), INTENT(IN) :: study_sizes(m)
-    REAL(8) :: y(n)
-    INTEGER(8) :: a, b, i
-
-    a = 1
-    b = 0
-    DO i = 1, m
-        b = b + study_sizes(i)
-        CALL block_dot_mv(v(a:b), Z(a:b, :), d, x(a:b), y(a:b), &
-            study_sizes(i), k)
-        a = a + study_sizes(i)
-    END DO
-
-END FUNCTION dot_mv
-
-
-FUNCTION dot_mm(v, Z, d, study_sizes, X, m, n, k, l) RESULT(Y)
-    IMPLICIT NONE
-    INTEGER(8), INTENT(IN) :: m, n, k, l
-    REAL(8), INTENT(IN) :: v(n), Z(n, k), d(k), X(n, l)
-    INTEGER(8), INTENT(IN) :: study_sizes(m)
-    REAL(8) :: Y(n, l)
-    INTEGER(8) :: a, b, i
-
-    a = 1
-    b = 0
-    DO i = 1, m
-        b = b + study_sizes(i)
-        CALL block_dot_mm(v(a:b), Z(a:b, :), d, X(a:b, :), Y(a:b, :), &
-            study_sizes(i), k, l)
-        a = a + study_sizes(i)
-    END DO
-
-END FUNCTION dot_mm
-
-
-FUNCTION invdot_mv(v, Z, d, study_sizes, x, m, n, k) RESULT(y)
-    IMPLICIT NONE
-    INTEGER(8), INTENT(IN) :: m, n, k
-    REAL(8), INTENT(IN) :: v(n), Z(n, k), d(k), x(n)
-    INTEGER(8), INTENT(IN) :: study_sizes(m)
-    REAL(8) :: y(n)
-    INTEGER(8) :: a, b, i
-
-    a = 1
-    b = 0
-    DO i = 1, m
-        b = b + study_sizes(i)
-        IF (n <= k) THEN
-            CALL block_invdot_mv_n(v(a:b), Z(a:b, :), d, x(a:b), y(a:b), &
-                study_sizes(i), k)
-        ELSE
-            CALL block_invdot_mv_k(v(a:b), Z(a:b, :), d, x(a:b), y(a:b), &
-                study_sizes(i), k)
-        END IF
-        a = a + study_sizes(i)
-    END DO
-
-END FUNCTION invdot_mv
-
-
-FUNCTION invdot_mm(v, Z, d, study_sizes, X, m, n, k, l) RESULT(Y)
-    IMPLICIT NONE
-    INTEGER(8), INTENT(IN) :: m, n, k, l
-    REAL(8), INTENT(IN) :: v(n), Z(n, k), d(k), X(n, l)
-    INTEGER(8), INTENT(IN) :: study_sizes(m)
-    REAL(8) :: Y(n, l)
-    INTEGER(8) :: a, b, i
-
-    a = 1
-    b = 0
-    DO i = 1, m
-        b = b + study_sizes(i)
-        IF (n <= k) THEN
-            CALL block_invdot_mm_n(v(a:b), Z(a:b, :), d, X(a:b, :), Y(a:b, :),&
-                study_sizes(i), k, l)
-        ELSE
-            CALL block_invdot_mm_k(v(a:b), Z(a:b, :), d, X(a:b, :), Y(a:b, :),&
-                study_sizes(i), k, l)
-        END IF
-        a = a + study_sizes(i)
-    END DO
-
-END FUNCTION invdot_mm
-
-
-FUNCTION logdet(v, Z, d, study_sizes, m, n, k) RESULT(val)
-    IMPLICIT NONE
-    INTEGER(8), INTENT(IN) :: m, n, k
-    REAL(8), INTENT(IN) :: v(n), Z(n,k), d(k)
-    INTEGER(8), INTENT(IN) :: study_sizes(m)
-    REAL(8) :: val
-    INTEGER(8) :: a, b, i
-
-    val = 0.D0
-    a = 1
-    b = 0
-    DO i = 1, m
-        b = b + study_sizes(i)
-        IF (n <= k) THEN
-            val = val + block_logdet_n(v(a:b), Z(a:b, :), d, study_sizes(i), k)
-        ELSE
-            val = val + block_logdet_k(v(a:b), Z(a:b, :), d, study_sizes(i), k)
-        END IF
-        a = a + study_sizes(i)
-    END DO
-
-END FUNCTION logdet
-
-! private SUBROUTINEs
-! -----------------------------------------------------------------------------
-SUBROUTINE block_dot_mv(v, Z, d, x, y, n, k)
-    IMPLICIT NONE
-    INTEGER(8), INTENT(IN) :: n, k
-    REAL(8), INTENT(IN) :: v(n), Z(n, k), d(k), x(n)
-    REAL(8), INTENT(INOUT) :: y(n)
-
-    y = v*x + MATMUL(Z, d*MATMUL(TRANSPOSE(Z), x))
-
-END SUBROUTINE block_dot_mv
-
-
-SUBROUTINE block_dot_mm(v, Z, d, X, Y, n, k, l)
-    IMPLICIT NONE
-    INTEGER(8), INTENT(IN) :: n, k, l
-    REAL(8), INTENT(IN) :: v(n), Z(n, k), d(k), X(n, l)
-    REAL(8), INTENT(INOUT) :: Y(n, l)
-    REAL(8) :: T(n, k)
-    INTEGER(8) :: i, j
-
-    DO j = 1, l
-        Y(:, j) = X(:, j)*v
-    END DO
-
-    DO i = 1, n
-        T(i, :) = Z(i, :)*d
-    END DO
-
-    Y = Y + MATMUL(T, MATMUL(TRANSPOSE(Z), X))
-
-END SUBROUTINE block_dot_mm
-
-
-SUBROUTINE block_invdot_mv_n(v, Z, d, x, y, n, k)
-    IMPLICIT NONE
-    INTEGER(8), INTENT(IN) :: n, k
-    REAL(8), INTENT(IN) :: v(n), Z(n, k), d(k), x(n)
-    REAL(8), INTENT(INOUT) :: y(n)
-    REAL(8) :: M(n, n), T(n, k)
-    INTEGER(8) :: i, info, ipiv(n)
-
-    ! compute covariance matrix
-    DO i = 1, n
-        T(i, :) = Z(i, :)*d
-    END DO
-
-    M = MATMUL(T, TRANSPOSE(Z))
-
-    DO i = 1, n
-        M(i, i) = M(i, i) + v(i)
-    END DO
-
-    ! solve linear system
-    y = x
-    CALL DGESV(n, 1, M, n, ipiv, y, n, info)
-
-END SUBROUTINE block_invdot_mv_n
-
-
-SUBROUTINE block_invdot_mv_k(v, Z, d, x, y, n, k)
-    IMPLICIT NONE
-    INTEGER(8), INTENT(IN) :: n, k
-    REAL(8), INTENT(IN) :: v(n), Z(n, k), d(k), x(n)
-    REAL(8), INTENT(INOUT) :: y(n)
-    REAL(8) :: M(k, k), T(n, k), u(k)
-    INTEGER(8) :: j, info, ipiv(k)
-
-    ! compute covariance matrix
-    DO j = 1, k
-        T(:, j) = Z(:, j)/v
-    END DO
-
-    M = MATMUL(TRANSPOSE(Z), T)
-
-    DO j = 1, k
-        M(j, j) = M(j, j) + 1.D0/d(j)
-    END DO
-
-    ! solve linear system
-    u = MATMUL(TRANSPOSE(T), x)
-    CALL DGESV(k, 1, M, k, ipiv, u, k, info)
-
-    y = x/v - MATMUL(T, u)
-
-END SUBROUTINE block_invdot_mv_k
-
-
-SUBROUTINE block_invdot_mm_n(v, Z, d, X, Y, n, k, l)
-    IMPLICIT NONE
-    INTEGER(8), INTENT(IN) :: n, k, l
-    REAL(8), INTENT(IN) :: v(n), Z(n, k), d(k), X(n, l)
-    REAL(8), INTENT(INOUT) :: Y(n, l)
-    REAL(8) :: M(n, n), T(n, k)
-    INTEGER(8) :: i, info, ipiv(n)
-
-    ! compute covariance matrix
-    DO i = 1, n
-        T(i, :) = Z(i, :)*d
-    END DO
-
-    M = MATMUL(T, TRANSPOSE(Z))
-
-    DO i = 1, n
-        M(i, i) = M(i, i) + v(i)
-    END DO
-
-    ! solve linear system
-    Y = X
-    CALL DGESV(n, l, M, n, ipiv, Y, n, info)
-
-END SUBROUTINE block_invdot_mm_n
-
-
-SUBROUTINE block_invdot_mm_k(v, Z, d, X, Y, n, k, l)
-    IMPLICIT NONE
-    INTEGER(8), INTENT(IN) :: n, k, l
-    REAL(8), INTENT(IN) :: v(n), Z(n, k), d(k), X(n, l)
-    REAL(8), INTENT(INOUT) :: Y(n, l)
-    REAL(8) :: M(k, k), T(n, k), U(k, l)
-    INTEGER(8) :: j, info, ipiv(k)
-
-    ! compute covariance matrix
-    DO j = 1, k
-        T(:, j) = Z(:, j)/v
-    END DO
-
-    M = MATMUL(TRANSPOSE(Z), T)
-
-    DO j = 1, k
-        M(j, j) = M(j, j) + 1.D0/d(j)
-    END DO
-
-    ! solve linear system
-    U = MATMUL(TRANSPOSE(T), X)
-    CALL DGESV(k, l, M, k, ipiv, U, k, info)
-
-    Y = MATMUL(T, U)
-
-    DO j = 1, l
-        Y(:, j) = X(:, j)/v - Y(:, j)
-    END DO
-
-END SUBROUTINE block_invdot_mm_k
-
-
-FUNCTION block_logdet_n(v, Z, d, n, k) RESULT(val)
-    IMPLICIT NONE
-    INTEGER(8), INTENT(IN) :: n, k
-    REAL(8), INTENT(IN) :: v(n), Z(n, k), d(k)
-    REAL(8) :: val
-    REAL(8) :: M(n, n), T(n, k)
-    INTEGER(8) :: i, info, ipiv(n)
-
-    ! compute covariance matrix
-    DO i = 1, n
-        T(i, :) = Z(i, :)*d
-    END DO
-
-    M = MATMUL(T, TRANSPOSE(Z))
-
-    DO i = 1, n
-        M(i, i) = M(i, i) + v(i)
-    END DO
-
-    ! compute the determinant
-    CALL DGETRF(n, n, M, n, ipiv, info)
-
-    val = 0.D0
-    DO i = 1, n
-        val = val + LOG(ABS(M(i, i)))
-    END DO
-
-END FUNCTION block_logdet_n
-
-
-FUNCTION block_logdet_k(v, Z, d, n, k) RESULT(val)
-    IMPLICIT NONE
-    INTEGER(8), INTENT(IN) :: n, k
-    REAL(8), INTENT(IN) :: v(n), Z(n, k), d(k)
-    REAL(8) :: val
-    REAL(8) :: M(k, k), T(n, k)
-    INTEGER(8) :: i, j, info, ipiv(k)
-
-    ! compute covariance matrix
-    DO i = 1, n
-        T(i, :) = Z(i, :)*d
-    END DO
-
-    DO j = 1, k
-        T(:, j) = T(:, j)/v
-    END DO
-
-    M = MATMUL(TRANSPOSE(Z), T)
-
-    DO j = 1, k
-        M(j, j) = M(j, j) + 1.D0
-    END DO
-
-    ! compute the determinant
-    CALL DGETRF(k, k, M, k, ipiv, info)
-
-    val = 0.D0
-    DO i = 1, n
-        val = val + LOG(v(i))
-    END DO
-
-    DO j = 1, k
-        val = val + LOG(ABS(M(j, j)))
-    END DO
-
-END FUNCTION block_logdet_k
-
-
-END MODULE varmat
\ No newline at end of file
diff --git a/pynlme/utils.py b/pynlme/utils.py
deleted file mode 100644
index 646e0a5..0000000
--- a/pynlme/utils.py
+++ /dev/null
@@ -1,110 +0,0 @@
-# utility classes and functions
-import numpy as np
-from scipy.linalg import block_diag
-from pynlme.futils import varmat
-
-
-class VarMat:
-    def __init__(self, v, z, gamma, study_sizes):
-        # check the dimension of the input
-        assert v.ndim == 1
-        assert v.size == sum(study_sizes)
-        assert np.all(v > 0.0)
-
-        assert z.ndim == 2
-
-        assert gamma.ndim == 1
-        assert gamma.size == z.shape[1]
-        assert np.all(gamma >= 0.0)
-
-        assert study_sizes.dtype == int
-        assert np.all(study_sizes > 0)
-
-        # pass in variables
-        self.v = v
-        self.z = np.asfortranarray(z)
-        self.gamma = gamma
-        self.study_sizes = study_sizes
-
-        self.k_gamma = gamma.size
-        self.num_studies = len(study_sizes)
-        self.num_data = sum(study_sizes)
-
-    # public functions
-    # -------------------------------------------------------------------------
-    def varMat(self):
-        split_idx = np.cumsum(self.study_sizes)[:-1]
-        v_study = np.split(self.v, split_idx)
-        z_study = np.split(self.z, split_idx, axis=0)
-
-        diag_blocks = [self._blockVarMat(v_study[i], z_study[i], self.gamma)
-                       for i in range(self.num_studies)]
-
-        return block_diag(*diag_blocks)
-
-    def invVarMat(self):
-        split_idx = np.cumsum(self.study_sizes)[:-1]
-        v_study = np.split(self.v, split_idx)
-        z_study = np.split(self.z, split_idx, axis=0)
-
-        diag_blocks = [self._blockInvVarMat(v_study[i], z_study[i], self.gamma)
-                       for i in range(self.num_studies)]
-
-        return block_diag(*diag_blocks)
-
-    def dot(self, x):
-        if x.ndim == 1:
-            return varmat.dot_mv(self.v, self.z, self.gamma,
-                                 self.study_sizes, x, self.num_studies,
-                                 self.num_data, self.k_gamma)
-        if x.ndim == 2:
-            x = np.asfortranarray(x)
-            return varmat.dot_mm(self.v, self.z, self.gamma,
-                                 self.study_sizes, x, self.num_studies,
-                                 self.num_data, self.k_gamma, x.shape[1])
-
-    def invDot(self, x):
-        if x.ndim == 1:
-            return varmat.invdot_mv(self.v, self.z, self.gamma,
-                                    self.study_sizes, x, self.num_studies,
-                                    self.num_data, self.k_gamma)
-        if x.ndim == 2:
-            x = np.asfortranarray(x)
-            return varmat.invdot_mm(self.v, self.z, self.gamma,
-                                    self.study_sizes, x, self.num_studies,
-                                    self.num_data, self.k_gamma, x.shape[1])
-
-    def logDet(self):
-        return varmat.logdet(self.v, self.z, self.gamma,
-                             self.study_sizes, self.num_studies,
-                             self.num_data, self.k_gamma)
-
-    @classmethod
-    def testProblem(cls):
-        num_studies = 5
-        study_sizes = np.array([10]*num_studies)
-        num_data = np.sum(study_sizes)
-
-        k = 2
-        s = np.random.rand(num_data)*0.09 + 0.01
-        v = s**2
-        z = np.random.randn(num_data, k)
-        gamma = np.random.rand(k)*0.1
-
-        return cls(v, z, gamma, study_sizes)
-
-    # internal functions
-    # -------------------------------------------------------------------------
-    @staticmethod
-    def _blockVarMat(v, z, gamma):
-        return np.diag(v) + (z*gamma).dot(z.T)
-
-    @staticmethod
-    def _blockInvVarMat(v, z, gamma):
-        if v.size <= gamma.size:
-            return np.linalg.inv(np.diag(v) + (z*gamma).dot(z.T))
-        else:
-            w = 1.0/v
-            y = z*w.reshape(w.size, 1)
-            d = np.diag(1.0/gamma) + y.T.dot(z)
-            return np.diag(w) - y.dot(np.linalg.inv(d).dot(y.T))
diff --git a/setup.py b/setup.py
new file mode 100644
index 0000000..1ae2eb2
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,52 @@
+import sys
+from pathlib import Path
+from setuptools import setup, find_packages
+
+
+if __name__ == '__main__':
+    base_dir = Path(__file__).parent
+    src_dir = base_dir/'src'/'limetr'
+
+    sys.path.insert(0, src_dir.as_posix())
+    import __about__ as about
+
+    with (base_dir/'README.rst').open() as f:
+        long_description = f.read()
+
+    install_requirements = [
+        'numpy',
+        'scipy',
+        'spmat',
+        'matplotlib'
+    ]
+
+    test_requirements = [
+        'pytest',
+        'pytest-mock',
+    ]
+
+    doc_requirements = []
+
+    setup(name=about.__title__,
+          version=about.__version__,
+
+          description=about.__summary__,
+          long_description=long_description,
+          license=about.__license__,
+          url=about.__uri__,
+
+          author=about.__author__,
+          author_email=about.__email__,
+
+          package_dir={'': 'src'},
+          packages=find_packages(where='src'),
+          include_package_data=True,
+
+          install_requires=install_requirements,
+          tests_require=test_requirements,
+          extras_require={
+              'docs': doc_requirements,
+              'test': test_requirements,
+              'dev': doc_requirements + test_requirements
+          },
+          zip_safe=False,)
diff --git a/src/limetr/__about__.py b/src/limetr/__about__.py
new file mode 100644
index 0000000..4740afd
--- /dev/null
+++ b/src/limetr/__about__.py
@@ -0,0 +1,16 @@
+__all__ = [
+    "__title__", "__summary__", "__uri__", "__version__", "__author__",
+    "__email__", "__license__", "__copyright__",
+]
+
+__title__ = "limetr"
+__summary__ = "limetr: linear mixed effects model with trimming."
+__uri__ = "https://github.com/zhengp0/limetr"
+
+__version__ = "0.1.0"
+
+__author__ = "Peng Zheng"
+__email__ = "zhengp@uw.eud"
+
+__license__ = "MIT License"
+__copyright__ = f"Copyright 2020 {__author__}"
diff --git a/src/limetr/__init__.py b/src/limetr/__init__.py
new file mode 100644
index 0000000..824db4c
--- /dev/null
+++ b/src/limetr/__init__.py
@@ -0,0 +1,2 @@
+from .core import LimeTr
+from .data import Data
\ No newline at end of file
diff --git a/src/limetr/core.py b/src/limetr/core.py
new file mode 100644
index 0000000..557a438
--- /dev/null
+++ b/src/limetr/core.py
@@ -0,0 +1,444 @@
+"""
+    model
+    ~~~~~
+
+    Main model module.
+"""
+from typing import Tuple, Dict
+
+import numpy as np
+from numpy import ndarray
+from scipy.linalg import block_diag
+from scipy.stats import norm
+from scipy.optimize import LinearConstraint, minimize
+from spmat import BDLMat
+
+from limetr.data import Data
+from limetr.utils import get_varmat, split_by_sizes
+from limetr.variable import FeVariable, ReVariable
+
+
+class LimeTr:
+    """
+    LimeTr model class
+
+    Attributes
+    ----------
+    data : Data
+        Data object contains the observations
+    fevar : FeVariable
+        Fixed effects variable
+    revar : ReVariable
+        Random effects variable
+    inlier_pct : float
+        Inlier percentage
+
+    Methods
+    -------
+    get_vars(var)
+        Split array into variables beta and gamma.
+    get_residual(beta)
+        Compute trimming weighted residual.
+    get_femat(beta)
+        Compute trimming weighted Jacobian matrix of fixed effects.
+    get_remat()
+        Compute trimming weighted design matrix of random effects.
+    get_obsvar()
+        Compute trimming weighted observation variance.
+    get_varmat(gamma)
+        Compute trimming weighted variance covariance matrix of the likelihood
+    objective(var)
+        Objective function of the optimiation problem
+    gradient(var)
+        Gradient function of the optimization problem
+    hessian(var)
+        Hessian function of the optimization problem, approximated by the Fisher
+        information matrix.
+    detect_outliers(var)
+        Detect outlier based on the statistical model and given varaible.
+    get_model_init()
+        Compute model initialization.
+    fit_model(var=None, num_tr_steps=3, options=None)
+        Run optimization algorithm to get solution.
+    get_random_effects(var)
+        Given estimate of beta and gamma, return the estimate of random effects.
+    """
+
+    def __init__(self,
+                 data: Data,
+                 fevar: FeVariable,
+                 revar: ReVariable,
+                 inlier_pct: float = 1.0):
+        """
+        Parameters
+        ----------
+        data : Data
+            Data object contains the observations.
+        fevar : FeVariable
+            Fixed effects variable.
+        revar : ReVariable
+            Random effects variable.
+        inlier_pct : float, optional
+            Inlier percentage, by default 1
+
+        Raises
+        ------
+        ValueError
+            When fixed effects shape not matching with the data.
+        ValueError
+            When random effects shape not matching with the data.
+        ValueError
+            When the inlier percentage is outside zero to one interval.
+        """
+
+        if data.num_obs != fevar.mapping.shape[0]:
+            raise ValueError("Fixed effects shape not matching with data.")
+        if data.num_obs != revar.mapping.shape[0]:
+            raise ValueError("Random effects shape not matching with data.")
+        if inlier_pct < 0 or inlier_pct > 1:
+            raise ValueError("`inlier_pct` must be between 0 and 1.")
+
+        self.data = data
+        self.fevar = fevar
+        self.revar = revar
+        self.inlier_pct = inlier_pct
+        self.result = None
+
+    # pylint:disable=unbalanced-tuple-unpacking
+    def get_vars(self, var: ndarray) -> Tuple[ndarray, ndarray]:
+        """
+        Split array into variables beta and gamma.
+
+        Parameters
+        ----------
+        var : ndarray
+
+        Returns
+        -------
+        Tuple[ndarray, ndarray]
+            Return beta and gamma.
+        """
+        beta, gamma = tuple(split_by_sizes(var, [self.fevar.size,
+                                                 self.revar.size]))
+        gamma[gamma < 0] = 0.0
+        return beta, gamma
+
+    def get_residual(self, beta: ndarray) -> ndarray:
+        """
+        Compute trimming weighted residual
+
+        Parameters
+        ----------
+        beta : ndarray
+
+        Returns
+        -------
+        ndarray
+            Weighted residual
+        """
+        return self.data.weight*(self.data.obs -
+                                 self.fevar.mapping(beta))
+
+    def get_femat(self, beta: ndarray) -> ndarray:
+        """
+        Compute trimming weighted Jacobian matrix of fixed effects.
+
+        Parameters
+        ----------
+        beta : ndarray
+
+        Returns
+        -------
+        ndarray
+            Weighted Jacobian matrix of fixed effects.
+        """
+        return self.data.weight[:, None]*self.fevar.mapping.jac(beta)
+
+    def get_remat(self) -> ndarray:
+        """
+        Compute trimming weighted design matrix of random effects.
+
+        Returns
+        -------
+        ndarray
+            Weighted design matrix of random effects.
+        """
+        return self.data.weight[:, None]*self.revar.mapping.mat
+
+    def get_obsvar(self) -> ndarray:
+        """
+        Compute trimming weighted observation variance
+
+        Returns
+        -------
+        ndarray
+            Weighted observation variance
+        """
+        return self.data.obs_se**(2*self.data.weight)
+
+    def get_varmat(self, gamma: ndarray) -> BDLMat:
+        """
+        Compute trimming weighted variance covariance matrix of the likelihood
+
+        Parameters
+        ----------
+        gamma : ndarray
+
+        Returns
+        -------
+        BDLMat
+            Weighted variance covariance matrix of the likelihood.
+        """
+        return get_varmat(gamma,
+                          split_by_sizes(self.get_obsvar(),
+                                         self.data.group_sizes),
+                          split_by_sizes(self.get_remat(),
+                                         self.data.group_sizes))
+
+    def objective(self, var: ndarray) -> float:
+        """
+        Objective function of the optimization problem.
+
+        Parameters
+        ----------
+        var : ndarray
+
+        Returns
+        -------
+        float
+            Objective function value.
+        """
+        beta, gamma = self.get_vars(var)
+        r = self.get_residual(beta)
+        d = self.get_varmat(gamma)
+
+        val = 0.5*(d.logdet() + r.dot(d.invdot(r)))
+        val += self.fevar.prior_objective(beta)
+        val += self.revar.prior_objective(gamma)
+
+        return val
+
+    def gradient(self, var: ndarray) -> ndarray:
+        """
+        Gradient function of the optimization problem
+
+        Parameters
+        ----------
+        var : ndarray
+
+        Returns
+        -------
+        ndarray
+            Gradient at given variable.
+        """
+        beta, gamma = self.get_vars(var)
+        r = self.get_residual(beta)
+        d = self.get_varmat(gamma)
+        femat = self.get_femat(beta)
+        remat = self.get_remat()
+
+        dr = d.invdot(r)
+        split_index = np.cumsum(np.insert(self.data.group_sizes, 0, 0))[:-1]
+
+        grad_beta = -femat.T.dot(dr) + self.fevar.prior_gradient(beta)
+        grad_gamma = 0.5*(
+            np.sum(remat*(d.invdot(remat)), axis=0) -
+            np.sum(np.add.reduceat(remat.T*dr, split_index, axis=1)**2, axis=1)
+        ) + self.revar.prior_gradient(gamma)
+
+        return np.hstack([grad_beta, grad_gamma])
+
+    def hessian(self, var: ndarray) -> ndarray:
+        """
+        Hessian function of the optimization problem, approximated by the Fisher
+        information matrix.
+
+        Parameters
+        ----------
+        var : ndarray
+
+        Returns
+        -------
+        ndarray
+            Hessian at given variable.
+        """
+        beta, gamma = self.get_vars(var)
+        d = self.get_varmat(gamma)
+        femat = self.get_femat(beta)
+        remat = split_by_sizes(self.get_remat(), self.data.group_sizes)
+
+        beta_fisher = femat.T.dot(d.invdot(femat))
+        beta_fisher += self.fevar.prior_hessian(beta)
+
+        gamma_fisher = np.zeros((self.revar.size, self.revar.size))
+        for i, dlmat in enumerate(d.dlmats):
+            gamma_fisher += 0.5*(remat[i].T.dot(dlmat.invdot(remat[i])))**2
+        gamma_fisher += self.revar.prior_hessian(gamma)
+
+        return block_diag(beta_fisher, gamma_fisher)
+
+    def detect_outliers(self, var: ndarray) -> ndarray:
+        """
+        Detect outlier based on the statistical model and given varaible
+
+        Parameters
+        ----------
+        var : ndarray
+
+        Returns
+        -------
+        ndarray
+            Indices of outliers.
+        """
+        beta, gamma = self.get_vars(var)
+        r = self.data.obs - self.fevar.mapping(beta)
+        s = np.sqrt(self.data.obs_se**2 +
+                    np.sum(self.revar.mapping.mat**2*gamma, axis=1))
+        a = norm.ppf(0.5 + 0.5*self.inlier_pct)
+        return np.abs(r) > a*s
+
+    def get_model_init(self) -> ndarray:
+        """
+        Get model initializations
+
+        Returns
+        -------
+        ndarray
+            Return the initialization variables.
+        """
+        beta = np.zeros(self.fevar.size)
+        gamma = np.zeros(self.revar.size)
+        var = np.hstack([beta, gamma])
+        grad_beta = self.gradient(var)[:self.fevar.size]
+        hess_beta = self.hessian(var)[:self.fevar.size,
+                                      :self.fevar.size]
+        beta = beta - np.linalg.solve(
+            hess_beta + np.identity(self.fevar.size),
+            grad_beta
+        )
+        return np.hstack([beta, gamma])
+
+    def _fit_model(self,
+                   var: ndarray = None,
+                   options: Dict = None):
+        """
+        (Inner) Fit model function
+
+        Parameters
+        ----------
+        var : ndarray, optional
+            Initial guess of the variable, by default None
+        options : dict, optional
+            scipy optimizer options, by default None
+        """
+        var = self.get_model_init() if var is None else var.copy()
+
+        bounds = np.hstack([self.fevar.get_uprior_info(),
+                            self.revar.get_uprior_info()]).T
+        constraints_mat = block_diag(self.fevar.get_linear_upriors_mat(),
+                                     self.revar.get_linear_upriors_mat())
+        constraints_vec = np.hstack([self.fevar.get_linear_upriors_info(),
+                                     self.revar.get_linear_upriors_info()])
+        constraints = [LinearConstraint(
+            constraints_mat,
+            constraints_vec[0],
+            constraints_vec[1]
+        )] if constraints_mat.size > 0 else []
+
+        self.result = minimize(self.objective, var,
+                               method="trust-constr",
+                               jac=self.gradient,
+                               hess=self.hessian,
+                               constraints=constraints,
+                               bounds=bounds,
+                               options=options)
+
+    def fit_model(self,
+                  var: ndarray = None,
+                  trim_steps: int = 3,
+                  options: Dict = None):
+        """
+        Fit model function
+
+        Parameters
+        ----------
+        var : ndarray, optional
+            Initial guess of the variable, by default None
+        trim_steps : int, optional
+            Number of trimming steps, by default 3
+        options : Dict, optional
+            scipy optimizer options, by default None
+
+        Raises
+        ------
+        ValueError
+            When ``num_tr_steps`` is strictly less than 2.
+        """
+
+        trim_steps = int(trim_steps)
+        if trim_steps < 2:
+            raise ValueError("At least two trimming steps.")
+
+        self._fit_model(var=var, options=options)
+        if self.inlier_pct < 1.0:
+            index = self.detect_outliers(self.result.x)
+            if index.sum() > 0:
+                for weight in np.linspace(1.0, 0.0, trim_steps)[1:]:
+                    self.data.weight.fill(1.0)
+                    self.data.weight[index] = weight
+                    self._fit_model(var=self.result.x, options=options)
+                    index = self.detect_outliers(self.result.x)
+
+    def get_random_effects(self, var: ndarray) -> ndarray:
+        """
+        Estimate random effects given beta and gamma
+
+        Parameters
+        ----------
+        var : ndarray
+
+        Returns
+        -------
+        ndarray
+            An array contains random effects.
+        """
+        beta, gamma = self.get_vars(var)
+        residual = split_by_sizes(self.get_residual(beta),
+                                  self.data.group_sizes)
+        obsvar = split_by_sizes(self.get_obsvar(),
+                                self.data.group_sizes)
+        remat = split_by_sizes(self.get_remat(),
+                               self.data.group_sizes)
+        random_effects = np.vstack([
+            gamma*np.linalg.solve(
+                (remat[i].T/obsvar[i]).dot(remat[i]*gamma) + np.identity(self.revar.size),
+                (remat[i].T/obsvar[i]).dot(residual[i])
+            )
+            for i in range(self.data.num_groups)
+        ])
+
+        return random_effects
+
+    @property
+    def soln(self) -> Dict[str, ndarray]:
+        """Solution summary"""
+        if self.result is None:
+            raise ValueError("Please fit the mdoel first.")
+        beta, gamma = self.get_vars(self.result.x)
+        beta_sd, gamma_sd = self.get_vars(1.0/np.sqrt(
+            np.diag(self.hessian(self.result.x))
+        ))
+        random_effects = self.get_random_effects(self.result.x)
+        return {
+            "beta": beta,
+            "gamma": gamma,
+            "beta_sd": beta_sd,
+            "gamma_sd": gamma_sd,
+            "random_effects": random_effects
+        }
+
+    def __repr__(self) -> str:
+        return (f"LimeTr(data={self.data},\n"
+                f"       fevar={self.fevar},\n"
+                f"       revar={self.revar},\n"
+                f"       inlier_pct={self.inlier_pct})")
diff --git a/src/limetr/data.py b/src/limetr/data.py
new file mode 100644
index 0000000..d23ef3a
--- /dev/null
+++ b/src/limetr/data.py
@@ -0,0 +1,110 @@
+"""
+    data
+    ~~~~
+
+    Data module.
+"""
+import operator
+from numbers import Number
+from typing import Iterable, Union
+
+import numpy as np
+
+from limetr.utils import default_vec_factory, iterable
+
+
+# pylint: disable=too-many-instance-attributes
+class Data:
+    """
+    Data containers observations and group information.
+
+    Attributes
+    ----------
+    obs : Vector
+        Observations. Assumed to be sorted by the group id.
+    obs_se : Vector
+        Standard deviations of observation. Default is one for each observation.
+    group_sizes : Vector
+        Number of observations for each group. Default treat every observation
+        as a group by itself.
+    weight : Vector
+        Weights for each observation. Default is one for each observation.
+    """
+
+    def __init__(self,
+                 obs: Iterable,
+                 obs_se: Union[Number, Iterable] = 1.0,
+                 group_sizes: Iterable[int] = None,
+                 weight: Union[Number, Iterable] = 1.0):
+        """
+        Parameters
+        ----------
+        obs : Iterable
+            Observations. Assumed to be sorted by the group id.
+        obs_se : Union[Number, Iterable], optional
+            Standard deviations of observation. Default is one.
+        group_sizes : Iterable[int], optional
+            Number of observations for each group. Default is ``None``.
+        weight : Union[Number, Iterable], optional
+            Weights for each observation. Default is one.
+        """
+        self.obs = obs
+        self.obs_se = obs_se
+        self.weight = weight
+        self.group_sizes = group_sizes
+
+    @property
+    def num_obs(self) -> int:
+        """Number of observations"""
+        return self.obs.size
+
+    @property
+    def num_groups(self) -> int:
+        """Number of groups"""
+        return self.group_sizes.size
+
+    obs = property(operator.attrgetter("_obs"))
+
+    @obs.setter
+    def obs(self, vec: Iterable):
+        vec = np.asarray(vec)
+        if any(np.isnan(vec)):
+            raise ValueError("`obs` must not containing nan(s).")
+        self._obs = vec
+
+    obs_se = property(operator.attrgetter("_obs_se"))
+
+    @obs_se.setter
+    def obs_se(self, vec: Union[Number, Iterable]):
+        vec = default_vec_factory(vec, self.num_obs, vec_name="obs_se")
+        if any(vec <= 0.0):
+            raise ValueError("`obs_se` must be all positive.")
+        self._obs_se = vec
+
+    weight = property(operator.attrgetter("_weight"))
+
+    @weight.setter
+    def weight(self, vec: Union[Number, Iterable]):
+        vec = default_vec_factory(vec, self.num_obs, vec_name="weight")
+        if any(vec < 0) or any(vec > 1):
+            raise ValueError("`weight` must be all between 0 and 1.")
+        self._weight = vec
+
+    group_sizes = property(operator.attrgetter("_group_sizes"))
+
+    @group_sizes.setter
+    def group_sizes(self, vec: Union[None, Iterable]):
+        if vec is None:
+            vec = np.ones(self.num_obs, dtype=int)
+        elif iterable(vec):
+            vec = np.asarray(vec).astype(int)
+        else:
+            raise TypeError("Use `None` or a vector to set `group_sizes`.")
+        if any(vec <= 0.0):
+            raise ValueError("`group_sizes` must be all positive.")
+        if sum(vec) != self.num_obs:
+            raise ValueError("Sum of `group_size` must equal to `num_obs`.")
+        self._group_sizes = vec
+
+    def __repr__(self) -> str:
+        return f"Data(num_obs={self.num_obs}, num_groups={self.num_groups})"
diff --git a/src/limetr/linalg/__init__.py b/src/limetr/linalg/__init__.py
new file mode 100644
index 0000000..3cd29c9
--- /dev/null
+++ b/src/limetr/linalg/__init__.py
@@ -0,0 +1,1 @@
+from .smooth_mapping import SmoothMapping, LinearMapping
diff --git a/src/limetr/linalg/smooth_mapping.py b/src/limetr/linalg/smooth_mapping.py
new file mode 100644
index 0000000..a5c5bae
--- /dev/null
+++ b/src/limetr/linalg/smooth_mapping.py
@@ -0,0 +1,116 @@
+"""
+    smooth_mapping
+    ~~~~~~~~~~~~~~
+
+    Smooth mapping module.
+"""
+from typing import Callable, Tuple
+
+import numpy as np
+
+
+class SmoothMapping:
+    """
+    Smooth mapping class, including function and Jacobian function.
+
+    Attributes
+    ----------
+    shape: Tuple[int, int]
+        Shape of the mapping.
+    fun: Callable
+        Mapping function.
+    jac: Callable
+        Jacobian function of the mapping.
+    """
+
+    def __init__(self,
+                 shape: Tuple[int, int],
+                 fun: Callable,
+                 jac: Callable):
+        """
+        Parameters
+        ----------
+        shape : Tuple[int, int]
+            Shape of the mapping.
+        fun : Callable
+            Mapping function.
+        jac : Callable
+            Jacobian function of the mapping.
+
+        Raises
+        ------
+        AssertionError
+            If ``shape`` is not a tuple with two positive intergers.
+        AssertionError
+            If ``fun`` is not callable.
+        AssertionError
+            If ``jac`` is not callable.
+        """
+        assert (isinstance(shape, tuple) and
+                len(shape) == 2 and
+                all([isinstance(size, int) and size > 0 for size in shape])), \
+            "`shape` has to be tuple with two positive integers."
+        assert callable(fun), "`fun` must be callable."
+        assert callable(jac), "`jac` must be callable."
+        self.shape = shape
+        self.fun = self._validate_input(fun)
+        self.jac = self._validate_input(jac)
+
+    def _validate_input(self, fun: Callable) -> Callable:
+        """
+        Decorator to validate the input before pass into the function.
+
+        Parameters
+        ----------
+        fun : Callable
+            Function to be decorated
+
+        Returns
+        -------
+        Callable
+            Function that valide the input.
+        """
+        def decorated_fun(x: np.ndarray) -> np.ndarray:
+            x = np.asarray(x)
+            if x.size != self.shape[1]:
+                raise ValueError("Input size not matching with mapping shape.")
+            return fun(x)
+        return decorated_fun
+
+    def __call__(self, x: np.ndarray) -> np.ndarray:
+        return self.fun(x)
+
+    def __repr__(self) -> str:
+        return f"SmoothMapping(shape={self.shape})"
+
+
+# pylint:disable=too-few-public-methods
+class LinearMapping(SmoothMapping):
+    """
+    Linear mapping class, construct smooth mapping from a matrix.
+
+    Attributes
+    ----------
+    mat: ndarray
+        Matrix as the linear mapping.
+    """
+
+    def __init__(self, mat: np.ndarray):
+        """
+        Parameters
+        ----------
+        mat : ndarray
+            Matrix as the linear mapping.
+        """
+        mat = np.asarray(mat)
+        assert mat.ndim == 2, "`mat` must be a matrix."
+        self.mat = mat
+
+        # pylint: disable=unused-argument
+        def fun(x): return mat.dot(x)
+        def jac(x): return mat
+
+        super().__init__(mat.shape, fun, jac)
+
+    def __repr__(self) -> str:
+        return f"LinearMapping(shape={self.shape})"
diff --git a/src/limetr/stats/__init__.py b/src/limetr/stats/__init__.py
new file mode 100644
index 0000000..eba59ee
--- /dev/null
+++ b/src/limetr/stats/__init__.py
@@ -0,0 +1,6 @@
+from .prior import (Prior,
+                    GaussianPrior,
+                    UniformPrior,
+                    LinearPrior,
+                    LinearGaussianPrior,
+                    LinearUniformPrior)
diff --git a/src/limetr/stats/prior.py b/src/limetr/stats/prior.py
new file mode 100644
index 0000000..d81857e
--- /dev/null
+++ b/src/limetr/stats/prior.py
@@ -0,0 +1,299 @@
+"""
+Prior Module
+"""
+from numbers import Number
+from typing import Iterable, List, Union, Any
+
+import numpy as np
+from limetr.utils import broadcast, get_maxlen
+
+
+class Prior:
+    """
+    Generic prior class, need to be inherited.
+
+    Attributes
+    ----------
+    info : ndarray
+        Information array of the prior.
+    size : int
+        Size of the prior.
+
+    Methods
+    -------
+    objective(var)
+        Objective function of optimization interface
+    gradient(var)
+        Gradient function of optimization interface
+    hessian(var)
+        Hessian function of the optimization interface
+    """
+
+    def __init__(self,
+                 info: List[Any] = None,
+                 size: int = 0):
+        """
+        Parameters
+        ----------
+        info : List[Any]
+            Information of the prior. Length of the list is number of
+            information components. Each component can be either scalar or a
+            vector. If there are more than one vector with size more than one,
+            the size of the vectors need to match. By default ``None``.
+        size : int, optional
+            Size of the prior, by default 0.
+        """
+        size = max(int(size), get_maxlen(info))
+        info = broadcast(info, size)
+        self.info = info
+        self.size = size
+
+    @property
+    def is_empty(self) -> bool:
+        """
+        Returns
+        -------
+        bool
+            If prior has zero size.
+        """
+        return self.size == 0
+
+    # pylint:disable=unused-argument
+    # pylint:disable=no-self-use
+    def objective(self, var: np.ndarray) -> float:
+        """
+        Objective function for optimiation interface.
+
+        Parameters
+        ----------
+        var : np.ndarray
+            Variable that prior is acting on.
+
+        Returns
+        -------
+        float
+            Objective value regarding the log likelihood of the prior.
+        """
+        return 0.0
+
+    def gradient(self, var: np.ndarray) -> np.ndarray:
+        """
+        Gradient function for optimiation interface.
+
+        Parameters
+        ----------
+        var : ndarray
+            Variable that prior is acting on.
+
+        Returns
+        -------
+        ndarray
+            Gradient value regarding the log likelihood of the prior.
+        """
+        return np.zeros(len(var))
+
+    def hessian(self, var: np.ndarray) -> np.ndarray:
+        """
+        Hessian function for optimiation interface.
+
+        Parameters
+        ----------
+        var : ndarray
+            Variable that prior is acting on.
+
+        Returns
+        -------
+        ndarray
+            Hessian value regarding the log likelihood of the prior.
+        """
+        return np.zeros((len(var), len(var)))
+
+    def __repr__(self) -> str:
+        return f"Prior(size={self.size})"
+
+
+class GaussianPrior(Prior):
+    """
+    Gaussian Prior
+
+    Attributes
+    ----------
+        mean : ndarray
+            Mean vector of the prior.
+        sd : ndarray
+            Standard deviation vector of the prior.
+    """
+
+    def __init__(self,
+                 mean: Union[Number, Iterable] = 0.0,
+                 sd: Union[Number, Iterable] = np.inf,
+                 size: int = 0):
+        """
+        Parameters
+        ----------
+        mean : Union[Number, Iterable], optional
+            Mean of the Gaussian prior, by default 0.
+        sd : Union[Number, Iterable], optional
+            Standard deviation of the Gaussian prior, by default inf.
+        size : int, optional
+            Size of the prior, by default 0.
+
+        Raises
+        ------
+        ValueError
+            If any standard deviations are less or equal to zero.
+        """
+        super().__init__([mean, sd], size=size)
+        if not all(self.info[1] > 0):
+            raise ValueError("Standard deviations have to be positive numbers.")
+        self.mean = self.info[0]
+        self.sd = self.info[1]
+
+    def objective(self, var: np.ndarray) -> float:
+        return 0.5*np.sum((var - self.mean)**2/self.sd**2)
+
+    def gradient(self, var: np.ndarray) -> np.ndarray:
+        return (var - self.mean)/self.sd**2
+
+    # pylint: disable=unused-argument
+    def hessian(self, var: np.ndarray) -> np.ndarray:
+        return np.diag(1/self.sd**2)
+
+    def __repr__(self) -> str:
+        return f"GaussianPrior(mean={self.mean}, sd={self.sd})"
+
+
+class UniformPrior(Prior):
+    """
+    Uniform Prior
+
+    Attributes
+    ----------
+        lb : ndarray
+            Lower bounds of the prior.
+        ub : ndarray
+            Upper bounds of the prior.
+    """
+
+    def __init__(self,
+                 lb: Union[Number, Iterable] = -np.inf,
+                 ub: Union[Number, Iterable] = np.inf,
+                 size: int = 0):
+        """
+        Parameters
+        ----------
+        lb : Union[Number, Iterable], optional
+            Lower bounds of the prior, by default -inf
+        ub : Union[Number, Iterable], optional
+            Upper bounds of the prior, by default inf
+        size : int, optional
+            Size of the prior, by default 0
+
+        Raises
+        ------
+        ValueError
+            If any lower bounds are greater than upper bounds.
+        """
+        super().__init__([lb, ub], size=size)
+        if any(self.info[0] > self.info[1]):
+            raise ValueError("Lower bounds must be less or equal than upper bounds.")
+        self.lb = self.info[0]
+        self.ub = self.info[1]
+
+    def __repr__(self) -> str:
+        return f"UniformPrior(lb={self.lb}, ub={self.ub})"
+
+
+class LinearPrior(Prior):
+    """
+    Linear Prior
+
+    Attributes
+    ----------
+    mat: ndarray
+        Linear mapping (matrix).
+    """
+
+    def __init__(self,
+                 mat: Iterable,
+                 info: List[Any]):
+        """
+        Parameters
+        ----------
+        mat : Iterable
+            Linear mapping (matrix).
+        info : List[Any]
+            Information array of the prior.
+        """
+        mat = np.asarray(mat)
+        if mat.ndim != 2:
+            raise ValueError("`mat` has to be a matrix.")
+        Prior.__init__(self, info, mat.shape[0])
+        self.mat = mat
+
+    def __repr__(self) -> str:
+        return f"LinearPrior(shape={self.mat.shape})"
+
+
+class LinearGaussianPrior(LinearPrior, GaussianPrior):
+    """
+    Linear Gaussian Prior
+    """
+
+    def __init__(self,
+                 mat: Iterable,
+                 mean: Union[Number, Iterable] = 0.0,
+                 sd: Union[Number, Iterable] = np.inf):
+        """
+        Parameters
+        ----------
+        mat: Iterable
+            Linear mapping(matrix).
+        mean: Union[Number, Iterable], optional
+            Mean of the prior, by default tuple with zero length.
+        sd: Union[Number, Iterable], optional
+            Standard deviation of the prior, by default tuple with zero length.
+        """
+        LinearPrior.__init__(self, mat, [mean, sd])
+        GaussianPrior.__init__(self, self.info[0], self.info[1], size=self.size)
+
+    def objective(self, var: np.ndarray) -> float:
+        trans_var = self.mat.dot(var)
+        return super().objective(trans_var)
+
+    def gradient(self, var: np.ndarray) -> np.ndarray:
+        trans_var = self.mat.dot(var)
+        return self.mat.T.dot(super().gradient(trans_var))
+
+    def hessian(self, var: np.ndarray) -> np.ndarray:
+        trans_var = self.mat.dot(var)
+        return self.mat.T.dot(super().hessian(trans_var).dot(self.mat))
+
+    def __repr__(self) -> str:
+        return f"LinearGaussianPrior(mean={self.mean}, sd={self.sd}, shape={self.mat.shape})"
+
+
+class LinearUniformPrior(LinearPrior, UniformPrior):
+    """
+    Linear Uniform Prior
+    """
+
+    def __init__(self,
+                 mat: Iterable,
+                 lb: Union[Number, Iterable] = -np.inf,
+                 ub: Union[Number, Iterable] = np.inf):
+        """
+        Parameters
+        ----------
+        mat: Iterable
+            Linear mapping(matrix).
+        lb: Union[Number, Iterable], optional
+            Lower bounds of the prior, by default tuple with zero length.
+        ub: Union[Number, Iterable], optional
+            Upper bounds of the prior, by default tuple with zero length.
+        """
+        LinearPrior.__init__(self, mat, [lb, ub])
+        UniformPrior.__init__(self, self.info[0], self.info[1], size=self.size)
+
+    def __repr__(self) -> str:
+        return f"LinearUniformPrior(lb={self.lb}, ub={self.ub}, shape={self.mat.shape})"
diff --git a/src/limetr/utils.py b/src/limetr/utils.py
new file mode 100644
index 0000000..114e8ae
--- /dev/null
+++ b/src/limetr/utils.py
@@ -0,0 +1,273 @@
+"""
+    utils
+    ~~~~~
+
+    Helper functions.
+"""
+from collections.abc import Iterable
+from numbers import Number
+from typing import Any, List, Union
+
+import numpy as np
+from spmat.dlmat import BDLMat, DLMat
+
+
+def split_by_sizes(array: np.ndarray,
+                   sizes: List[int],
+                   axis: int = 0) -> List[np.ndarray]:
+    """
+    Function that split an array into a list of arrays, provided the size for
+    each sub-array size.
+
+    Parameters
+    ----------
+    array : ndarray
+        The array need to be splitted.
+    sizes : List[int]
+        A list of sizes for each sub-array.
+    axis: int, optional
+        Along which axis, array will be splitted, default is 0.
+
+    Raises
+    ------
+    AssertionError
+        If the sum of the ``sizes`` does not equal to the shape of ``array``
+        along ``axis``.
+
+    Returns
+    -------
+    List[ndarray]
+        A list of splitted array. 
+    """
+    assert array.shape[axis] == sum(sizes)
+    return np.split(array, np.cumsum(sizes)[:-1], axis=axis)
+
+
+def empty_array() -> np.ndarray:
+    """
+    Function used for 'default_factory', creates and returns empty array.
+
+    Returns
+    -------
+    ndarray
+        An empty array with ``dtype`` being ``float``.
+    """
+    return np.array([])
+
+
+def default_vec_factory(vec: Union[Number, Iterable],
+                        size: int,
+                        default_value: Any = None,
+                        vec_name: str = 'vector') -> np.ndarray:
+    """
+    Function that automatically create and fill values of a vector.
+
+    Parameters
+    ----------
+    vec : Union[Number, Iterable]
+        A vector or number that need to be checked or expand.
+    size : int
+        The desired size of the vector.
+    default_value : Any, optional
+        Default value of the vector, will be used when ``vec`` is empty.
+        Default is ``None``.
+    vec_name : str, optional
+        Name of the vector, for more informative error message.
+        Default to be ``'vector'``.
+
+    Raises
+    ------
+    AssertionError
+        If ``vec`` is empty and ``default_value`` is ``None``.
+    AssertionError
+        If ``vec`` is ``Iterable`` and length does not equal to ``size``.
+
+    Returns
+    -------
+    ndarray:
+        Final processed array.
+    """
+    if np.isscalar(vec):
+        vec = np.repeat(vec, size)
+    elif len(vec) == 0:
+        assert default_value is not None, \
+            "Must provide `default_value` when `vec` is empty."
+        vec = np.repeat(default_value, size)
+    else:
+        vec = np.asarray(vec)
+        check_size(vec, size, vec_name=vec_name)
+
+    return vec
+
+
+def check_size(vec: Iterable, size: int, vec_name: str = 'vector'):
+    """
+    Function that check the size consistency.
+
+    Parameters
+    ----------
+    vec : Iterable
+        Iterable vector which length will be checked.
+    size : int
+        Desired size of the vector.
+    vec_name : str, optional
+        Name of the vector, for more informative error message.
+        Default to be ``'vector'``.
+
+    Raises
+    ------
+    ValueError
+        If vector length does not equal to provided ``size``.
+    """
+    if len(vec) != size:
+        raise ValueError(f"{vec_name} must be length {size}.")
+
+
+def iterable(__obj: object) -> bool:
+    """
+    Function that check if an object is iterable.
+
+    Parameters
+    ----------
+
+    __obj : object
+        Object to be examed.
+
+    Returns
+    -------
+    bool
+        ``True`` if object is iterable, and ``False`` otherwise.
+    """
+    return isinstance(__obj, Iterable)
+
+
+def has_no_repeat(array: np.ndarray) -> bool:
+    """
+    Function that check if an array have no repeat values.
+
+    Parameters
+    ----------
+    array : ndarray
+        Array that need to be examed.
+
+    Returns
+    -------
+    bool
+        ``True`` if array has no repeat values, and ``False`` otherwise.
+    """
+    return array.size == np.unique(array).size
+
+
+def sizes_to_slices(sizes: Iterable) -> List[slice]:
+    """
+    Function that convert sizes of sub-arrays to corresponding slices in the
+    original array.
+
+    Parameters
+    ----------
+    sizes : Iterable[int]
+        Iterable object contains positive integers as the sizes of the arrays.
+
+    Returns
+    -------
+    List[slice]
+        A list of ``slice`` to access each sub-array in the original array.
+    """
+    ends = np.cumsum(sizes)
+    starts = np.insert(ends, 0, 0)[:-1]
+    return [slice(*pair) for pair in zip(starts, ends)]
+
+
+def get_varmat(gamma: np.ndarray,
+               obsvar: List[np.ndarray],
+               remat: List[np.ndarray]) -> BDLMat:
+    """
+    Function that compute the variance-covariance matrix for the mixed effects
+    model, used in optimization interface.
+
+    Parameters
+    ----------
+    gamma : ndarray
+        Variance vector of the random effects.
+    obsvar : List[ndarray]
+        Variance vectors of the obervation errors. Each array in the list
+        corresponding to a group.
+    remat : List[ndarray]
+        Random effects design matrices. Each matrix in the list corresponding to
+        a group.
+
+    Raises
+    ------
+    AssertionError
+        If length of ``obsvar`` does not equal to length of ``remat``. Both
+        lengths should agree with number of groups.
+
+    Returns
+    -------
+    BDLMat
+        Block diagonal matrix object.
+    """
+    assert len(obsvar) == len(remat)
+    sqrt_gamma = np.sqrt(gamma)
+    dlmats = [
+        DLMat(obsvar[i], remat[i]*sqrt_gamma)
+        for i in range(len(obsvar))
+    ]
+    return BDLMat(dlmats)
+
+
+def get_maxlen(objs: List[Any]) -> int:
+    """
+    Get the maximum len of a list of objects.
+
+    Parameters
+    ----------
+    objs : List[Any]
+        A list of objects.
+
+    Returns
+    -------
+    int
+        Maximum length among objects.
+    """
+    return max([len(obj) if iterable(obj) else 1 for obj in objs])
+
+
+def broadcast(objs: List[Any],
+              size: int) -> np.ndarray:
+    """
+    Broadcast a list of objects.
+
+    Parameters
+    ----------
+    objs : List[Any]
+        A list of objects.
+    size : int
+        Size for the broadcasting.
+
+    Raises
+    ------
+    ValueError
+        If there is iterable object in the list whose size is not 1 and not
+        agree with broadcast size.
+
+    Returns
+    -------
+    np.ndarray
+        Two dimensional array that stores the squared objects.
+    """
+    size = int(size)
+    assert size >= 0, "Size has to be a non-negative integer."
+    if size == 0:
+        vecs = np.empty(shape=(len(objs), 0))
+    else:
+        for i, obj in enumerate(objs):
+            if np.isscalar(obj) or len(obj) == 1:
+                objs[i] = np.repeat(obj, size)
+            elif len(obj) == size:
+                objs[i] = np.asarray(obj)
+            else:
+                raise ValueError("Object size not consistent with broadcast "
+                                 "size.")
+        vecs = np.vstack(objs)
+    return vecs
diff --git a/src/limetr/variable.py b/src/limetr/variable.py
new file mode 100644
index 0000000..8f2ba89
--- /dev/null
+++ b/src/limetr/variable.py
@@ -0,0 +1,352 @@
+"""
+Variable Module
+"""
+from typing import Any, Iterable, Type
+
+import numpy as np
+from numpy import ndarray
+
+from limetr.linalg import LinearMapping, SmoothMapping
+from limetr.stats import (GaussianPrior, LinearGaussianPrior, LinearPrior,
+                          LinearUniformPrior, Prior, UniformPrior)
+
+
+class Variable:
+    """
+    Variable class, contains mapping and prior information
+
+    Attributes
+    ----------
+    mapping : SmoothMapping
+        Smooth mapping, map the variable to predict the data.
+    priors : Iterable[Prior], optional
+        Priors for the variable.
+    name : Any
+        Name of the variable.
+    gprior : GaussianPrior
+        Gaussian prior of the variable.
+    uprior : UniformPrior
+        Uniform prior of the variable.
+    linear_gpriors : List[LinearGaussianPrior]
+        List of linear Gaussian priors.
+    linear_upriors : List[LinearUniformPrior]
+        List of linear Uniform priors.
+
+    Methods
+    -------
+    update_gprior(prior)
+        Update the Gaussian prior.
+    update_uprior(prior)
+        Update the Uniform prior.
+    update_linear_gpriors(prior)
+        Update linear Gaussian prior, append to the list.
+    update_linear_upriors(prior)
+        Update linear Uniform prior, append to the list.
+    update_priors(priors)
+        Update all priors providing a list of priors.
+    reset_priors()
+        Reset priors to default settings.
+    prior_objective(var)
+        Objective function from Gaussian prior and linear Gaussian priors.
+    prior_gradient(var)
+        Gradient function from Gaussian prior and linear Gaussian priors.
+    prior_hessian(var)
+        Hessian function from Gaussian prior and linear Gaussian priors.
+    get_uprior_info()
+        Return Uniform prior upper and lower bounds.
+    get_linear_upriors_mat()
+        Return linear Uniform prior linear mapping
+    get_linear_upriors_info()
+        Return linear Uniform upper and lower bounds.
+    """
+
+    def __init__(self,
+                 mapping: SmoothMapping,
+                 priors: Iterable[Prior] = (),
+                 name: Any = "unknown"):
+        """
+        Parameters
+        ----------
+        mapping : SmoothMapping
+            Smooth mapping, map the variable to predict the data.
+        priors : Iterable[Prior], optional
+            Priors for the variable, by default tuple with zero length.
+        name : Any, optional
+            Name of the variable, default by ``'unknown'``.
+
+        Raises
+        ------
+        TypeError
+            If ``mapping`` is not ``SmoothMapping``.
+        """
+        if not isinstance(mapping, SmoothMapping):
+            raise TypeError("`mapping` has to be SmoothMapping.")
+        self.mapping = mapping
+        self.name = name
+
+        self.gprior = GaussianPrior(size=self.size)
+        self.uprior = UniformPrior(size=self.size)
+        self.linear_gpriors = []
+        self.linear_upriors = []
+        self.update_priors(priors)
+
+    @property
+    def size(self) -> int:
+        """Size of the variable"""
+        return self.mapping.shape[1]
+
+    def _validate_prior(self,
+                        prior: Prior,
+                        prior_type: Type) -> Prior:
+        if not isinstance(prior, prior_type):
+            raise TypeError(f"`prior` must be type {prior_type.__name__}.")
+        size_matching = prior.mat.shape[1] == self.size \
+            if issubclass(prior_type, LinearPrior) else prior.size == self.size
+        if not size_matching:
+            raise ValueError("`prior` size not match with variable.")
+
+        return prior
+
+    def update_gprior(self, prior: GaussianPrior):
+        """
+        Update Gaussian prior.
+
+        Parameters
+        ----------
+        prior : GaussianPrior
+        """
+        self.gprior = self._validate_prior(prior, GaussianPrior)
+
+    def update_uprior(self, prior: UniformPrior):
+        """
+        Update Uniform prior
+
+        Parameters
+        ----------
+        prior : UniformPrior
+        """
+        self.uprior = self._validate_prior(prior, UniformPrior)
+
+    def update_linear_gpriors(self, prior: LinearGaussianPrior):
+        """
+        Update linear Gaussian priors
+
+        Parameters
+        ----------
+        prior : LinearGaussianPrior
+        """
+        self.linear_gpriors.append(self._validate_prior(
+            prior, LinearGaussianPrior
+        ))
+
+    def update_linear_upriors(self, prior: LinearUniformPrior):
+        """
+        Update linear Uniform priors
+
+        Parameters
+        ----------
+        prior : LinearUniformPrior
+        """
+        self.linear_upriors.append(self._validate_prior(
+            prior, LinearUniformPrior
+        ))
+
+    def update_priors(self, priors: Iterable[Prior]):
+        """
+        Update priors
+
+        Parameters
+        ----------
+        priors : Iterable[Prior]
+
+        Raises
+        ------
+        TypeError
+            When prior type is not recognizable.
+        """
+        for prior in priors:
+            if isinstance(prior, LinearGaussianPrior):
+                self.update_linear_gpriors(prior)
+            elif isinstance(prior, LinearUniformPrior):
+                self.update_linear_upriors(prior)
+            elif isinstance(prior, GaussianPrior):
+                self.update_gprior(prior)
+            elif isinstance(prior, UniformPrior):
+                self.update_uprior(prior)
+            else:
+                raise TypeError(f"Unrecognize prior type {type(prior).__name__}")
+
+    def reset_priors(self):
+        """
+        Reset prior to default settings
+        """
+        self.gprior = GaussianPrior(size=self.size)
+        self.uprior = UniformPrior(size=self.size)
+        self.linear_gpriors = []
+        self.linear_upriors = []
+
+    def prior_objective(self, var: ndarray) -> float:
+        """
+        Objective function from Gaussian prior and linear Gaussian priors.
+
+        Parameters
+        ----------
+        var : ndarray
+            Variables.
+
+        Returns
+        -------
+        float
+            Objective function value.
+        """
+        val = self.gprior.objective(var)
+        for prior in self.linear_gpriors:
+            val += prior.objective(var)
+        return val
+
+    def prior_gradient(self, var: ndarray) -> ndarray:
+        """
+        Gradient function from Gaussian prior and linear Gaussian priors.
+
+        Parameters
+        ----------
+        var : ndarray
+            Variables.
+
+        Returns
+        -------
+        ndarray
+            Gradient at given value.
+        """
+        val = self.gprior.gradient(var)
+        for prior in self.linear_gpriors:
+            val += prior.gradient(var)
+        return val
+
+    def prior_hessian(self, var: ndarray) -> ndarray:
+        """
+        Hessian function from Gaussian prior and linear Gaussian priors.
+
+        Parameters
+        ----------
+        var : ndarray
+            Variables.
+
+        Returns
+        -------
+        ndarray
+            Hessian at given value.
+        """
+        val = self.gprior.hessian(var)
+        for prior in self.linear_gpriors:
+            val += prior.hessian(var)
+        return val
+
+    def get_uprior_info(self) -> ndarray:
+        """
+        Get Uniform prior information
+
+        Returns
+        -------
+        ndarray
+            Lower and upper bounds of the prior.
+        """
+        if self.uprior is None:
+            uprior = UniformPrior(size=self.size)
+        else:
+            uprior = self.uprior
+        return uprior.info
+
+    def get_linear_upriors_mat(self) -> ndarray:
+        """
+        Get linear Uniform prior linear mapping
+
+        Returns
+        -------
+        ndarray
+            Linear mapping of the prior.
+        """
+        if len(self.linear_upriors) == 0:
+            mat = np.empty((0, self.size))
+        else:
+            mat = np.vstack([prior.mat for prior in self.linear_upriors])
+        return mat
+
+    def get_linear_upriors_info(self) -> ndarray:
+        """
+        Get linear Unform prior information
+
+        Returns
+        -------
+        ndarray
+            Lower and upper bounds of the prior.
+        """
+        if len(self.linear_upriors) == 0:
+            info = np.empty((2, 0))
+        else:
+            info = np.hstack([prior.info for prior in self.linear_upriors])
+        return info
+
+    def __repr__(self) -> str:
+        return f"Variable(name={self.name}, size={self.size})"
+
+
+class FeVariable(Variable):
+    """
+    Fixed effects variable.
+    """
+
+    def __init__(self,
+                 mapping: SmoothMapping,
+                 priors: Iterable[Prior] = (),
+                 name: Any = "fixed effects"):
+        """
+        Parameters
+        ----------
+        name : Any, optional
+            Name of the variable, by default "fixed effects"
+        """
+        super().__init__(mapping, priors, name=name)
+
+    def __repr__(self) -> str:
+        return f"FeVariable(size={self.size})"
+
+
+class ReVariable(Variable):
+    """
+    Random effects variable.
+    """
+
+    def __init__(self,
+                 mapping: SmoothMapping,
+                 priors: Iterable[Prior] = (),
+                 name: Any = "random effects"):
+        """
+        Parameters
+        ----------
+        name : Any, optional
+            Name of the variable, by default "random effects"
+        """
+        if not isinstance(mapping, LinearMapping):
+            raise TypeError("Random effect design mapping has to be linear.")
+
+        priors = list(priors)
+        if not any([isinstance(prior, UniformPrior) for prior in priors]):
+            priors.append(UniformPrior(lb=0.0, ub=np.inf, size=mapping.mat.shape[1]))
+
+        super().__init__(mapping, priors, name=name)
+
+    def update_uprior(self, prior: UniformPrior):
+        prior = self._validate_prior(prior, UniformPrior)
+        if any(prior.lb < 0):
+            raise ValueError("Random effects variance must have non-negative lower bounds.")
+        self.uprior = prior
+
+    def reset_priors(self):
+        self.gprior = GaussianPrior(size=self.size)
+        self.uprior = UniformPrior(lb=0.0, ub=np.inf, size=self.size)
+        self.linear_gpriors = []
+        self.linear_upriors = []
+
+    def __repr__(self) -> str:
+        return f"ReVariable(size={self.size})"
diff --git a/src/limetr/visual.py b/src/limetr/visual.py
new file mode 100644
index 0000000..24a38c4
--- /dev/null
+++ b/src/limetr/visual.py
@@ -0,0 +1,85 @@
+"""
+Visualization module
+"""
+from typing import Tuple
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.pyplot import Axes
+from scipy.stats import norm
+from limetr.core import LimeTr
+
+
+def funnel(model: LimeTr,
+           ax: Axes = None,
+           sd_ub: float = 1.0,
+           ui_bounds: Tuple[float, float] = None,
+           summarize: bool = True) -> Axes:
+    """
+    Create funnel plot of the model
+
+    Parameters
+    ----------
+    model : LimeTr
+        Given fitting model
+    ax : Axes, optional
+        Plot axis, by default ``None``. If ``None`` will create new axis.
+    sd_ub : float, optional
+        Quantile upper bound of total standard deviation, by default 1.
+    ui_bounds : Tuple[float, float], optional
+        Quantile bounds of uncertainty interval, by default ``None``.
+        If ``None`` it will be infered from ``model.inlier_pct``.
+    summarize : bool, optional
+        If ``summarize``, a text box with information will be printed out,
+        default by ``True``.
+
+    Returns
+    -------
+    Axes
+        Axis used to plot the funnel.
+    """
+    # inputs
+    ax = plt.subplots()[1] if ax is None else ax
+    if ui_bounds is None:
+        ui_bounds = (0.5 - 0.5*model.inlier_pct,
+                     0.5 + 0.5*model.inlier_pct)
+    color_data = "gray"
+    color_funnel = "#77C3E3"
+    # compute the residual and total standard error
+    soln = model.soln
+    r = model.data.obs - model.fevar.mapping(soln["beta"])
+    s = np.sqrt(model.data.obs_se**2 +
+                np.sum(model.revar.mapping.mat**2*soln["gamma"], axis=1))
+    # plot the data
+    trim_index = model.data.weight == 0
+    ax.scatter(r, s, color=color_data, alpha=0.5)
+    ax.scatter(r[trim_index], s[trim_index], color="red", marker="x", alpha=0.5)
+    # plot funnel
+    s_range = np.array([0.0, np.quantile(s, sd_ub)])
+    m_range = norm.ppf(ui_bounds)
+    ax.set_ylim(s_range[1], s_range[0])
+    for m in m_range:
+        ax.plot(s_range*m, s_range, color=color_funnel)
+    ax.fill_betweenx(s_range,
+                     s_range*m_range[0],
+                     s_range*m_range[1], color=color_funnel, alpha=0.3)
+    ax.axvline(0.0, linestyle="--", linewidth=1.0, color="black")
+    # add label
+    ax.set_xlabel("residual")
+    ax.set_ylabel("standard deviation")
+    # summarize info
+    if summarize:
+        props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
+        ax.text(1.05, 0.95,
+                f"num obs={model.data.num_obs}\n"
+                f"num groups={model.data.num_groups}\n"
+                f"inlier pct={model.inlier_pct:.3f}\n"
+                f"num outliers={trim_index.sum():d}\n"
+                f"beta={soln['beta']}\n"
+                f"gamma={soln['gamma']}\n"
+                f"funnel ui=({ui_bounds[0]:.3f}, {ui_bounds[1]:.3f})",
+                linespacing=1.5,
+                transform=ax.transAxes,
+                verticalalignment='top',
+                bbox=props)
+
+    return ax


Test patch:
diff --git a/tests/check_all.py b/tests/check_all.py
deleted file mode 100644
index 404ea83..0000000
--- a/tests/check_all.py
+++ /dev/null
@@ -1,38 +0,0 @@
-# test suite for xspline
-import os
-import sys
-# add current directory
-sys.path.append('./')
-
-
-def run_test(name):
-    namespace = {}
-    exec('import ' + name, namespace)
-    exec('ok = ' + name + '.' + name + '()', namespace)
-    ok = namespace['ok']
-    if ok:
-        print(name + ': OK')
-    else:
-        print(name + ': Error')
-    return ok
-
-
-fun_list = [
-	'utils_varmat_dot',
-	'utils_varmat_invDot',
-	'utils_varmat_logDet'
-]
-
-error_count = 0
-
-for name in fun_list:
-    ok = run_test(name)
-    if not ok:
-        error_count += 1
-
-if error_count > 0:
-    print('check_all: error_count =', error_count)
-    sys.exit(1)
-else:
-    print('check_all: OK')
-    sys.exit(0)
diff --git a/tests/test_core.py b/tests/test_core.py
new file mode 100644
index 0000000..4d054a4
--- /dev/null
+++ b/tests/test_core.py
@@ -0,0 +1,176 @@
+"""
+Test core module
+"""
+import pytest
+import numpy as np
+
+from limetr.data import Data
+from limetr.linalg import LinearMapping
+from limetr.variable import FeVariable, ReVariable
+from limetr.core import LimeTr
+from limetr.utils import split_by_sizes
+
+
+# pylint:disable=redefined-outer-name
+
+
+def ad_jacobian(fun, x, shape, eps=1e-10):
+    n = len(x)
+    c = x + 0j
+    g = np.zeros(shape)
+    for i in range(n):
+        c[i] += eps*1j
+        g[i] = fun(c).imag/eps
+        c[i] -= eps*1j
+    return g
+
+
+@pytest.fixture
+def data():
+    np.random.seed(123)
+    obs = np.random.randn(5)
+    obs_se = 0.1 + 0.1*np.random.rand(5)
+    group_sizes = np.array([2, 2, 1])
+    return Data(obs=obs,
+                obs_se=obs_se,
+                group_sizes=group_sizes)
+
+
+@pytest.fixture
+def fevar():
+    np.random.seed(123)
+    mat = np.random.randn(5, 2)
+    return FeVariable(LinearMapping(mat))
+
+
+@pytest.fixture
+def revar():
+    np.random.seed(123)
+    mat = np.random.randn(5, 2)
+    return ReVariable(LinearMapping(mat))
+
+
+@pytest.fixture
+def inlier_pct():
+    return 0.95
+
+
+@pytest.fixture
+def model(data, fevar, revar, inlier_pct):
+    return LimeTr(data, fevar, revar, inlier_pct=inlier_pct)
+
+
+@pytest.mark.parametrize("my_fevar",
+                         [FeVariable(LinearMapping(np.ones((4, 3))))])
+def test_validate_fevar(data, my_fevar, revar, inlier_pct):
+    with pytest.raises(ValueError):
+        LimeTr(data, my_fevar, revar, inlier_pct=inlier_pct)
+
+
+@pytest.mark.parametrize("my_revar",
+                         [ReVariable(LinearMapping(np.ones((4, 3))))])
+def test_validate_revar(data, fevar, my_revar, inlier_pct):
+    with pytest.raises(ValueError):
+        LimeTr(data, fevar, my_revar, inlier_pct=inlier_pct)
+
+
+@pytest.mark.parametrize("my_inlier_pct", [-0.1, 1.1])
+def test_validate_inlier_pct(data, fevar, revar, my_inlier_pct):
+    with pytest.raises(ValueError):
+        LimeTr(data, fevar, revar, inlier_pct=my_inlier_pct)
+
+
+@pytest.mark.parametrize("var",
+                         [np.hstack([np.ones(2), np.zeros(2)]),
+                          np.hstack([np.ones(2), np.full(2, -0.1)])])
+def test_get_vars(model, var):
+    beta, gamma = model.get_vars(var)
+    assert beta.size == model.fevar.size
+    assert gamma.size == model.revar.size
+    assert np.allclose(beta, 1.0)
+    assert np.allclose(gamma, 0.0)
+
+
+@pytest.mark.parametrize("beta", [np.ones(2)])
+def test_get_residual(model, beta):
+    my_residual = (model.data.obs - model.fevar.mapping(beta))
+    assert np.allclose(model.get_residual(beta), my_residual)
+
+
+@pytest.mark.parametrize("beta", [np.ones(2)])
+def test_get_femat(model, beta):
+    my_femat = model.fevar.mapping.jac(beta)
+    assert np.allclose(model.get_femat(beta), my_femat)
+
+
+def test_get_remat(model):
+    my_remat = model.revar.mapping.mat
+    assert np.allclose(model.get_remat(), my_remat)
+
+
+def test_get_obsvar(model):
+    my_obsvar = model.data.obs_se**2
+    assert np.allclose(model.get_obsvar(), my_obsvar)
+
+
+@pytest.mark.parametrize("gamma", [np.zeros(2)])
+def test_get_varmat(model, gamma):
+    my_varmat = np.diag(model.get_varmat(gamma).mat)
+    assert np.allclose(my_varmat, model.get_obsvar())
+
+
+@pytest.mark.parametrize("var", [np.hstack([np.ones(2), np.zeros(2)])])
+def test_objective(model, var):
+    beta, _ = model.get_vars(var)
+    my_obj = model.objective(var)
+    tr_obj = 0.5*(np.sum(model.get_residual(beta)**2/model.get_obsvar()) +
+                  np.sum(np.log(model.get_obsvar())))
+
+    assert np.isclose(my_obj, tr_obj)
+
+
+@pytest.mark.parametrize("var", [np.hstack([np.ones(2), np.zeros(2)])])
+def test_gradient(model, var):
+    my_grad = model.gradient(var)
+    tr_grad = ad_jacobian(model.objective, var, (var.size,))
+
+    assert np.allclose(my_grad, tr_grad)
+
+
+@pytest.mark.parametrize("var", [np.hstack([np.zeros(2), np.zeros(2)])])
+def test_hessian(model, var):
+    my_hess = model.hessian(var)
+    eigvals = np.linalg.eigvals(my_hess)
+
+    assert all(eigvals >= 0.0)
+
+
+def test_get_model_init(model):
+    var = model.get_model_init()
+    beta, gamma = model.get_vars(var)
+    assert beta.size == model.fevar.size
+    assert gamma.size == model.revar.size
+    assert all(gamma >= 0.0)
+
+
+@pytest.mark.parametrize("var", [np.hstack([np.zeros(2), np.zeros(2)])])
+def test_get_random_effects_zero_gamma(model, var):
+    random_effects = model.get_random_effects(var)
+    assert np.allclose(random_effects, 0.0)
+
+
+@pytest.mark.parametrize("var", [np.hstack([np.zeros(2), np.ones(2)])])
+def test_get_random_effects_nonzero_gamma(model, var):
+    beta, gamma = model.get_vars(var)
+    u = model.get_random_effects(var)
+    r = split_by_sizes(model.get_residual(beta), model.data.group_sizes)
+    v = split_by_sizes(model.get_obsvar(), model.data.group_sizes)
+    z = split_by_sizes(model.get_remat(), model.data.group_sizes)
+
+    def grad_fun(z, v, r, gamma, u):
+        return (z.T/v).dot(z.dot(u) - r) + u/gamma
+
+    assert np.allclose(np.hstack([
+        grad_fun(z[i], v[i], r[i], gamma, u[i])
+        for i in range(model.data.num_groups)
+    ]), 0.0)
diff --git a/tests/test_data.py b/tests/test_data.py
new file mode 100644
index 0000000..e5c4f4c
--- /dev/null
+++ b/tests/test_data.py
@@ -0,0 +1,87 @@
+"""
+Test Data Module
+"""
+import pytest
+import numpy as np
+from limetr.data import Data
+
+
+# pylint:disable=redefined-outer-name
+
+
+@pytest.fixture
+def obs():
+    return np.random.randn(5)
+
+
+@pytest.fixture
+def obs_se():
+    return np.ones(5)
+
+
+@pytest.fixture
+def group_sizes():
+    return np.array([2, 2, 1])
+
+
+@pytest.fixture
+def weight():
+    return np.ones(5)
+
+
+@pytest.mark.parametrize("obs_bad", [[1.0, 1.0, 1.0, 1.0, np.nan]])
+def test_obs_validate(obs_bad, obs_se, group_sizes, weight):
+    with pytest.raises(ValueError):
+        Data(obs_bad, obs_se, group_sizes, weight)
+
+
+@pytest.mark.parametrize("obs_se_good", [1.0, 2.0, [1.0, 2.0, 3.0, 4.0, 5.0]])
+def test_obs_se(obs, obs_se_good, group_sizes, weight):
+    d = Data(obs, obs_se_good, group_sizes, weight)
+    assert d.obs_se.size == d.num_obs
+    assert all(d.obs_se > 0.0)
+
+
+@pytest.mark.parametrize("obs_se_bad", [0.0,
+                                        [1.0, 2.0, 3.0, 4.0, -5.0],
+                                        [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]])
+def test_obs_se_validate(obs, obs_se_bad, group_sizes, weight):
+    with pytest.raises(ValueError):
+        Data(obs, obs_se_bad, group_sizes, weight)
+
+
+@pytest.mark.parametrize("weight_good", [1.0, 0.0, [0.5, 0.5, 0.5, 0.5, 0.5]])
+def test_weight(obs, obs_se, group_sizes, weight_good):
+    d = Data(obs, obs_se, group_sizes, weight_good)
+    assert d.weight.size == d.num_obs
+    assert all(d.weight >= 0.0) and all(d.weight <= 1.0)
+
+
+@pytest.mark.parametrize("weight_bad", [-1.0,
+                                        [1.0, 2.0, 1.0, 1.0, 1.0],
+                                        [1.0, 1.0, 1.0, 1.0, 1.0, 1.0]])
+def test_weight_validate(obs, obs_se, group_sizes, weight_bad):
+    with pytest.raises(ValueError):
+        Data(obs, obs_se, group_sizes, weight_bad)
+
+
+@pytest.mark.parametrize("group_sizes_good", [None,
+                                              [2, 2, 1],
+                                              [2.0, 2.0, 1.0]])
+def test_group_sizes(obs, obs_se, group_sizes_good, weight):
+    d = Data(obs, obs_se, group_sizes_good, weight)
+    assert d.group_sizes.size == d.num_groups
+    assert sum(d.group_sizes) == d.num_obs
+
+
+@pytest.mark.parametrize("group_sizes_bad", [[2, 2, 1, 1],
+                                             [1, -1, 2, 2, 1]])
+def test_group_sizes_validate(obs, obs_se, group_sizes_bad, weight):
+    with pytest.raises(ValueError):
+        Data(obs, obs_se, group_sizes_bad, weight)
+
+
+@pytest.mark.parametrize("group_sizes_bad", [1.0])
+def test_group_sizes_validate_type(obs, obs_se, group_sizes_bad, weight):
+    with pytest.raises(TypeError):
+        Data(obs, obs_se, group_sizes_bad, weight)
diff --git a/tests/test_linalg.py b/tests/test_linalg.py
new file mode 100644
index 0000000..27518b9
--- /dev/null
+++ b/tests/test_linalg.py
@@ -0,0 +1,44 @@
+# Test linear algebra module
+import numpy as np
+import pytest
+import limetr.linalg as linalg
+
+
+#pylint: disable=redefined-outer-name
+
+
+@pytest.fixture
+def smooth_mapping():
+    m = 5
+    n = 2
+    mat1 = np.random.randn(m, n)
+    mat2 = np.random.randn(m, n)
+
+    def fun(x):
+        return np.exp(mat1.dot(x)) - np.exp(mat2.dot(x))
+
+    def jac(x):
+        return (mat1.T*np.exp(mat1.dot(x)) - mat2.T*np.exp(mat2.dot(x))).T
+
+    return linalg.SmoothMapping((m, n), fun, jac)
+
+
+@pytest.fixture
+def linear_mapping():
+    m = 5
+    n = 2
+    mat = np.random.randn(m, n)
+    return linalg.LinearMapping(mat)
+
+
+def test_smooth_mapping(smooth_mapping):
+    x = np.zeros(smooth_mapping.shape[1])
+    assert smooth_mapping(x).size == smooth_mapping.shape[0]
+    assert smooth_mapping.jac(x).shape == smooth_mapping.shape
+
+
+def test_linear_mapping(linear_mapping):
+    x = np.zeros(linear_mapping.shape[1])
+    assert linear_mapping.mat.shape == linear_mapping.shape
+    assert linear_mapping(x).size == linear_mapping.shape[0]
+    assert linear_mapping.jac(x).shape == linear_mapping.shape
diff --git a/tests/test_stats.py b/tests/test_stats.py
new file mode 100644
index 0000000..947d40c
--- /dev/null
+++ b/tests/test_stats.py
@@ -0,0 +1,129 @@
+"""
+Test Statistic Module
+"""
+import pytest
+import numpy as np
+from limetr.stats.prior import (GaussianPrior,
+                                UniformPrior,
+                                LinearPrior,
+                                LinearGaussianPrior,
+                                LinearUniformPrior)
+
+
+def test_gprior_validate():
+    with pytest.raises(ValueError):
+        GaussianPrior(mean=1.0, sd=-1.0)
+
+
+@pytest.mark.parametrize("mean", [np.array([1.0, 2.0, 3.0])])
+@pytest.mark.parametrize("sd", [np.array([1.0, 1.0, 1.0])])
+@pytest.mark.parametrize("var", [np.array([1.0, 1.0, 1.0])])
+def test_gprior_objective(mean, sd, var):
+    prior = GaussianPrior(mean=mean, sd=sd)
+    my_obj = prior.objective(var)
+    tr_obj = 0.5*np.sum((mean - var)**2/sd**2)
+    assert my_obj == tr_obj
+
+
+@pytest.mark.parametrize("mean", [np.array([1.0, 2.0, 3.0])])
+@pytest.mark.parametrize("sd", [np.array([1.0, 1.0, 1.0])])
+@pytest.mark.parametrize("var", [np.array([1.0, 1.0, 1.0])])
+def test_gprior_gradient(mean, sd, var):
+    prior = GaussianPrior(mean=mean, sd=sd)
+    my_grad = prior.gradient(var)
+    tr_grad = (var - mean)/sd**2
+    assert np.allclose(my_grad, tr_grad)
+
+
+@pytest.mark.parametrize("mean", [np.array([1.0, 2.0, 3.0])])
+@pytest.mark.parametrize("sd", [np.array([1.0, 1.0, 1.0])])
+@pytest.mark.parametrize("var", [np.array([1.0, 1.0, 1.0])])
+def test_gprior_hessian(mean, sd, var):
+    prior = GaussianPrior(mean=mean, sd=sd)
+    my_hess = prior.hessian(var)
+    tr_hess = np.diag(1.0/sd**2)
+    assert np.allclose(my_hess, tr_hess)
+
+
+def test_uprior_validate():
+    with pytest.raises(ValueError):
+        UniformPrior(lb=1.0, ub=-1.0)
+
+
+@pytest.mark.parametrize("lb", [np.array([1.0, 1.0, 1.0])])
+@pytest.mark.parametrize("ub", [np.array([1.0, 2.0, 3.0])])
+@pytest.mark.parametrize("var", [np.array([1.0, 1.0, 1.0])])
+def test_uprior_objective(lb, ub, var):
+    prior = UniformPrior(lb=lb, ub=ub)
+    my_obj = prior.objective(var)
+    tr_obj = 0.0
+    assert my_obj == tr_obj
+
+
+@pytest.mark.parametrize("lb", [np.array([1.0, 1.0, 1.0])])
+@pytest.mark.parametrize("ub", [np.array([1.0, 2.0, 3.0])])
+@pytest.mark.parametrize("var", [np.array([1.0, 1.0, 1.0])])
+def test_uprior_gradient(lb, ub, var):
+    prior = UniformPrior(lb=lb, ub=ub)
+    my_grad = prior.gradient(var)
+    tr_grad = 0.0
+    assert np.allclose(my_grad, tr_grad)
+
+
+@pytest.mark.parametrize("lb", [np.array([1.0, 1.0, 1.0])])
+@pytest.mark.parametrize("ub", [np.array([1.0, 2.0, 3.0])])
+@pytest.mark.parametrize("var", [np.array([1.0, 1.0, 1.0])])
+def test_uprior_hessian(lb, ub, var):
+    prior = UniformPrior(lb=lb, ub=ub)
+    my_hess = prior.hessian(var)
+    tr_hess = 0.0
+    assert np.allclose(my_hess, tr_hess)
+
+
+@pytest.mark.parametrize(("mat, info"),
+                         [(np.empty(shape=(0, 2)), np.empty(shape=(2, 0))),
+                          (np.ones((3, 4)), np.ones((2, 3)))])
+def test_linear_prior(mat, info):
+    prior = LinearPrior(mat, info)
+    assert prior.mat.shape[0] == prior.size
+
+
+@pytest.mark.parametrize("mat", [np.ones((3, 4))])
+@pytest.mark.parametrize("mean", [np.array([1.0, 2.0, 3.0])])
+@pytest.mark.parametrize("sd", [np.array([1.0, 1.0, 1.0])])
+@pytest.mark.parametrize("var", [np.array([1.0, 1.0, 1.0, 1.0])])
+def test_linear_gprior_objective(mat, mean, sd, var):
+    prior = LinearGaussianPrior(mat=mat, mean=mean, sd=sd)
+    my_obj = prior.objective(var)
+    tr_obj = 0.5*np.sum((mat.dot(var) - mean)**2/sd**2)
+    assert my_obj == tr_obj
+
+
+@pytest.mark.parametrize("mat", [np.ones((3, 4))])
+@pytest.mark.parametrize("mean", [np.array([1.0, 2.0, 3.0])])
+@pytest.mark.parametrize("sd", [np.array([1.0, 1.0, 1.0])])
+@pytest.mark.parametrize("var", [np.array([1.0, 1.0, 1.0, 1.0])])
+def test_linear_gprior_gradient(mat, mean, sd, var):
+    prior = LinearGaussianPrior(mat=mat, mean=mean, sd=sd)
+    my_grad = prior.gradient(var)
+    tr_grad = mat.T.dot((mat.dot(var) - mean)/sd**2)
+    assert np.allclose(my_grad, tr_grad)
+
+
+@pytest.mark.parametrize("mat", [np.ones((3, 4))])
+@pytest.mark.parametrize("mean", [np.array([1.0, 2.0, 3.0])])
+@pytest.mark.parametrize("sd", [np.array([1.0, 1.0, 1.0])])
+@pytest.mark.parametrize("var", [np.array([1.0, 1.0, 1.0, 1.0])])
+def test_linear_gprior_hessian(mat, mean, sd, var):
+    prior = LinearGaussianPrior(mat=mat, mean=mean, sd=sd)
+    my_hess = prior.hessian(var)
+    tr_hess = mat.T.dot(np.diag(1.0/sd**2).dot(mat))
+    assert np.allclose(my_hess, tr_hess)
+
+
+@pytest.mark.parametrize("mat", [np.ones((3, 4))])
+@pytest.mark.parametrize("lb", [np.array([1.0, 1.0, 1.0])])
+@pytest.mark.parametrize("ub", [np.array([1.0, 2.0, 3.0])])
+def test_linear_uprior_hessian(mat, lb, ub):
+    prior = LinearUniformPrior(mat=mat, lb=lb, ub=ub)
+    assert mat.shape[0] == prior.size
diff --git a/tests/test_utils.py b/tests/test_utils.py
new file mode 100644
index 0000000..a56a22a
--- /dev/null
+++ b/tests/test_utils.py
@@ -0,0 +1,90 @@
+# Test helper function in utils module
+import numpy as np
+import pytest
+import limetr.utils as utils
+
+
+@pytest.mark.parametrize('vec', [np.arange(6)])
+@pytest.mark.parametrize('sizes', [[1, 2, 3], [3, 2, 1]])
+def test_split_by_sizes(vec, sizes):
+    vecs = utils.split_by_sizes(vec, sizes)
+    assert all([vecs[i].size == size for i, size in enumerate(sizes)])
+
+
+def test_empty_array():
+    array = utils.empty_array()
+    assert array.size == 0
+    assert np.issubdtype(array.dtype, float)
+
+
+@pytest.mark.parametrize('vec', [[0, 1, 2]])
+@pytest.mark.parametrize('size', [2])
+def test_check_size_validate(vec, size):
+    with pytest.raises(ValueError):
+        utils.check_size(vec, size)
+
+
+def test_check_size():
+    utils.check_size([1, 2, 3], 3)
+
+
+@pytest.mark.parametrize(('obj', 'result'),
+                         [(3, False),
+                          ([3], True)])
+def test_iterable(obj, result):
+    assert utils.iterable(obj) == result
+
+
+@pytest.mark.parametrize(('array', 'result'),
+                         [(np.array([1, 1, 2]), False),
+                          (np.array([1, 2, 3]), True)])
+def test_has_no_repeat(array, result):
+    assert utils.has_no_repeat(array) == result
+
+
+def test_sizes_to_slices():
+    sizes = [1, 2, 3]
+    slices = [slice(0, 1), slice(1, 3), slice(3, 6)]
+    result = utils.sizes_to_slices(sizes)
+    assert all([result[i] == slices[i] for i in range(len(sizes))])
+
+
+@pytest.mark.parametrize(('vec', 'size', 'default_value'),
+                         [([], 5, 1),
+                          (1, 5, None),
+                          ([1]*5, 5, None)])
+@pytest.mark.parametrize('result', [np.ones(5)])
+def test_default_vec_factory(vec, size, default_value, result):
+    my_result = utils.default_vec_factory(vec, size, default_value)
+    assert np.allclose(my_result, result)
+
+
+def test_get_varmat():
+    np.random.seed(123)
+    gamma = np.full(5, 0.1)
+    obsvar = [np.full(6, (i + 1)/10) for i in range(3)]
+    remat = [np.random.randn(6, 5) for i in range(3)]
+
+    varmat = utils.get_varmat(gamma, obsvar, remat)
+    assert np.isclose(varmat.logdet(), np.log(np.linalg.det(varmat.mat)))
+
+
+@pytest.mark.parametrize("objs", [[[1.0, 2.0, 3.0], 2.0, (2.0, 1.0)]])
+def test_get_maxlen(objs):
+    assert utils.get_maxlen(objs) == 3
+
+
+@pytest.mark.parametrize("objs", [[[1.0, 2.0, 3.0], 2.0, (1.0,)]])
+def test_broadcast(objs):
+    my_result = utils.broadcast(objs, utils.get_maxlen(objs))
+    assert np.allclose(my_result,
+                       np.array([[1.0, 2.0, 3.0],
+                                 [2.0, 2.0, 2.0],
+                                 [1.0, 1.0, 1.0]]))
+
+
+@pytest.mark.parametrize("size", [0])
+@pytest.mark.parametrize("objs", [[[1.0, 2.0, 3.0], 2.0, (1.0,)]])
+def test_broadcast_size_zero(size, objs):
+    my_result = utils.broadcast(objs, size)
+    assert my_result.shape == (len(objs), 0)
diff --git a/tests/test_variable.py b/tests/test_variable.py
new file mode 100644
index 0000000..0def5a1
--- /dev/null
+++ b/tests/test_variable.py
@@ -0,0 +1,169 @@
+"""
+Test Variable Module
+"""
+import pytest
+import numpy as np
+from limetr.linalg import LinearMapping
+from limetr.stats import GaussianPrior, UniformPrior, LinearGaussianPrior, LinearUniformPrior
+from limetr.variable import Variable, FeVariable, ReVariable
+
+
+# pylint: disable=redefined-outer-name
+
+
+@pytest.fixture
+def mapping():
+    return LinearMapping(np.ones((3, 2)))
+
+
+@pytest.fixture
+def var(mapping):
+    return Variable(mapping)
+
+
+@pytest.fixture
+def gprior():
+    return GaussianPrior(mean=-1.0, size=2)
+
+
+@pytest.fixture
+def uprior():
+    return UniformPrior(lb=-1.0, size=2)
+
+
+@pytest.fixture
+def linear_gprior():
+    return LinearGaussianPrior(mat=np.ones((5, 2)), mean=-1.0)
+
+
+@pytest.fixture
+def linear_uprior():
+    return LinearUniformPrior(mat=np.ones((5, 2)), lb=-1.0)
+
+
+@pytest.mark.parametrize("mapping", [1, "a"])
+def test_var_validate_mapping(mapping):
+    with pytest.raises(TypeError):
+        Variable(mapping)
+
+
+def test_var_default(var):
+    assert var.size == 2
+    assert var.gprior.size == 2
+    assert var.uprior.size == 2
+    assert len(var.linear_gpriors) == 0
+    assert len(var.linear_upriors) == 0
+
+
+def test_var_update_gprior(var, gprior):
+    var.update_gprior(gprior)
+    assert np.allclose(var.gprior.mean, -1.0)
+
+
+def test_var_update_uprior(var, uprior):
+    var.update_uprior(uprior)
+    assert np.allclose(var.uprior.lb, -1.0)
+
+
+def test_var_update_linear_gpriors(var, linear_gprior):
+    var.update_linear_gpriors(linear_gprior)
+    assert len(var.linear_gpriors) == 1
+    assert np.allclose(var.linear_gpriors[0].mat, 1.0)
+    assert np.allclose(var.linear_gpriors[0].mean, -1.0)
+
+
+def test_var_update_linear_upriors(var, linear_uprior):
+    var.update_linear_upriors(linear_uprior)
+    assert len(var.linear_upriors) == 1
+    assert np.allclose(var.linear_upriors[0].mat, 1.0)
+    assert np.allclose(var.linear_upriors[0].lb, -1.0)
+
+
+def test_var_update_priors(var, gprior, uprior, linear_gprior, linear_uprior):
+    var.update_priors([gprior, uprior, linear_gprior, linear_uprior])
+    assert np.allclose(var.gprior.mean, -1.0)
+    assert np.allclose(var.uprior.lb, -1.0)
+    assert len(var.linear_gpriors) == 1
+    assert len(var.linear_upriors) == 1
+
+
+def test_var_reset_priors(var, gprior, uprior, linear_gprior, linear_uprior):
+    var.update_priors([gprior, uprior, linear_gprior, linear_uprior])
+    var.reset_priors()
+    assert np.allclose(var.gprior.mean, 0.0)
+    assert all(np.isneginf(var.uprior.lb))
+    assert len(var.linear_gpriors) == 0
+    assert len(var.linear_upriors) == 0
+
+
+@pytest.mark.parametrize("x", [np.array([1.0, 2.0])])
+def test_var_prior_objective(var, gprior, linear_gprior, x):
+    var.update_priors([gprior, linear_gprior])
+    my_val = var.prior_objective(x)
+    tr_val = gprior.objective(x) + linear_gprior.objective(x)
+    assert np.isclose(my_val, tr_val)
+
+
+@pytest.mark.parametrize("x", [np.array([1.0, 2.0])])
+def test_var_prior_gradient(var, gprior, linear_gprior, x):
+    var.update_priors([gprior, linear_gprior])
+    my_val = var.prior_gradient(x)
+    tr_val = gprior.gradient(x) + linear_gprior.gradient(x)
+    assert np.allclose(my_val, tr_val)
+
+
+@pytest.mark.parametrize("x", [np.array([1.0, 2.0])])
+def test_var_prior_hessian(var, gprior, linear_gprior, x):
+    var.update_priors([gprior, linear_gprior])
+    my_val = var.prior_hessian(x)
+    tr_val = gprior.hessian(x) + linear_gprior.hessian(x)
+    assert np.allclose(my_val, tr_val)
+
+
+def test_var_get_uprior_info(var, uprior):
+    my_val = var.get_uprior_info()
+    assert np.allclose(my_val, np.array([[-np.inf]*2, [np.inf]*2]))
+    var.update_priors([uprior])
+    my_val = var.get_uprior_info()
+    assert np.allclose(my_val, uprior.info)
+
+
+def test_var_get_linear_upriors_mat(var, linear_uprior):
+    my_val = var.get_linear_upriors_mat()
+    assert my_val.shape == (0, 2)
+    var.update_priors([linear_uprior])
+    my_val = var.get_linear_upriors_mat()
+    assert np.allclose(my_val, linear_uprior.mat)
+
+
+def test_var_get_linear_upriors_info(var, linear_uprior):
+    my_val = var.get_linear_upriors_info()
+    assert my_val.shape == (2, 0)
+    var.update_priors([linear_uprior])
+    my_val = var.get_linear_upriors_info()
+    assert np.allclose(my_val, linear_uprior.info)
+
+
+def test_fevar(mapping):
+    var = FeVariable(mapping)
+    assert var.name == "fixed effects"
+
+
+def test_revar(mapping):
+    var = ReVariable(mapping)
+    assert np.allclose(var.uprior.lb, 0.0)
+
+
+def test_revar_update_uprior_validate(mapping, uprior):
+    var = ReVariable(mapping)
+    with pytest.raises(ValueError):
+        var.update_priors([uprior])
+
+
+def test_revar_reset_priors(mapping):
+    var = ReVariable(mapping)
+    uprior = UniformPrior(lb=1.0, size=2)
+    var.update_priors([uprior])
+    assert all(var.uprior.lb == 1.0)
+    var.reset_priors()
+    assert all(var.uprior.lb == 0.0)
diff --git a/tests/utils_varmat_dot.py b/tests/utils_varmat_dot.py
deleted file mode 100644
index dfbe249..0000000
--- a/tests/utils_varmat_dot.py
+++ /dev/null
@@ -1,46 +0,0 @@
-# test function dot
-
-def utils_varmat_dot():
-    import numpy as np
-    from pynlme.utils import VarMat
-
-    ok = True
-    # setup test problem
-    # -------------------------------------------------------------------------
-    var_mat = VarMat.testProblem()
-
-    tol = 1e-10
-
-    # test the dot product with vector
-    # -------------------------------------------------------------------------
-    x = np.random.randn(var_mat.num_data)
-
-    tr_y = var_mat.varMat().dot(x)
-    my_y = var_mat.dot(x)
-
-    err = np.linalg.norm(tr_y - my_y)
-    ok_vec = err < tol
-
-    if not ok_vec:
-        print('err', err)
-        print('tr_y', tr_y)
-        print('my_y', my_y)
-
-    # test the dot product with matrix
-    # -------------------------------------------------------------------------
-    x = np.random.randn(var_mat.num_data, 2)
-
-    tr_y = var_mat.varMat().dot(x)
-    my_y = var_mat.dot(x)
-
-    err = np.linalg.norm(tr_y - my_y)
-    ok_mat = err < tol
-
-    if not ok_mat:
-        print('err', err)
-        print('tr_y', tr_y)
-        print('my_y', my_y)
-
-    ok = ok and ok_vec and ok_mat
-
-    return ok
diff --git a/tests/utils_varmat_invDot.py b/tests/utils_varmat_invDot.py
deleted file mode 100644
index f835357..0000000
--- a/tests/utils_varmat_invDot.py
+++ /dev/null
@@ -1,46 +0,0 @@
-# test function invDot
-
-def utils_varmat_invDot():
-    import numpy as np
-    from pynlme.utils import VarMat
-
-    ok = True
-    # setup test problem
-    # -------------------------------------------------------------------------
-    var_mat = VarMat.testProblem()
-
-    tol = 1e-8
-
-    # test the invDot product with vector
-    # -------------------------------------------------------------------------
-    x = np.random.randn(var_mat.num_data)
-
-    tr_y = np.linalg.solve(var_mat.varMat(), x)
-    my_y = var_mat.invDot(x)
-
-    err = np.linalg.norm(tr_y - my_y)
-    ok_vec = err < tol
-
-    if not ok_vec:
-        print('err', err)
-        print('tr_y', tr_y)
-        print('my_y', my_y)
-
-    # test the invDot product with matrix
-    # -------------------------------------------------------------------------
-    x = np.random.randn(var_mat.num_data, 2)
-
-    tr_y = np.linalg.solve(var_mat.varMat(), x)
-    my_y = var_mat.invDot(x)
-
-    err = np.linalg.norm(tr_y - my_y)
-    ok_mat = err < tol
-
-    if not ok_mat:
-        print('err', err)
-        print('tr_y', tr_y)
-        print('my_y', my_y)
-
-    ok = ok and ok_vec and ok_mat
-
-    return ok
diff --git a/tests/utils_varmat_logDet.py b/tests/utils_varmat_logDet.py
deleted file mode 100644
index 85f5793..0000000
--- a/tests/utils_varmat_logDet.py
+++ /dev/null
@@ -1,27 +0,0 @@
-# test function logDet
-
-def utils_varmat_logDet():
-    import numpy as np
-    from pynlme.utils import VarMat
-
-    ok = True
-    # setup test problem
-    # -------------------------------------------------------------------------
-    var_mat = VarMat.testProblem()
-
-    tol = 1e-8
-
-    # test the logDet
-    # -------------------------------------------------------------------------
-    tr_y = np.log(np.linalg.det(var_mat.varMat()))
-    my_y = var_mat.logDet()
-
-    err = np.linalg.norm(tr_y - my_y)
-    ok = ok and err < tol
-
-    if not ok:
-        print('err', err)
-        print('tr_y', tr_y)
-        print('my_y', my_y)
-
-    return ok