wishart-moments-calculator

A SAGE project for calculating (symbolically and numerically) special moments of a Wishart random variable.

https://github.com/antunescarles/wishart-moments-calculator

Repository

A SAGE project for calculating (symbolically and numerically) special moments of a Wishart random variable.

README.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# WishartMoments Package\n",
    "## _A package to compute all invariant Wishart moments_\n",
    "\n",
    "\n",
    "WishartMoments is a python Package using the mathematical software SAGE for computing all the invariant moments of a Wishart distributed variable.\n",
    "\n",
    "## Features\n",
    "\n",
    "- Compute the symbolic expression for the expectation of order k of of the variable W in terms of its parameters n and \\Sigma.\n",
    "- Evaluate the moments for concrete values of k, n, and \\Sigma.\n",
    "- Compute (symbolically) the moments for the inverse W^{-1} when possible, in terms of n, \\Sigma and the size of the matrix.\n",
    "- When usign the package inside a Jupyter Notebook, you can print the results as MathJax formulas to get a nice visualisation.\n",
    "- Obtain the LaTex code that represents any of these expressions.\n",
    "- You can compute the moments using the [Wishart Moment Calculator][WM] web interface.\n",
    " \n",
    "## Installation\n",
    "You should have already installed SAGE. If not, you can get it from [here][sage].\n",
    "\n",
    "Open a SAGE shell and install WishartMoments using the SAGE pip utility by typing\n",
    "\n",
    "```sh\n",
    "sage -pip install WishartMoments\n",
    "```\n",
    "Once the installation is complete, to use the package in the command line you have to type `sage` in the SAGE shell to initiate SAGE, and once you see the prompt `sage:`, import the package (as any other Python package):\n",
    "```sh\n",
    "sage: import WishartMoments as wm\n",
    "```\n",
    "\n",
    "## How to use it\n",
    "\n",
    "We will now show  how to compute the following Wishart moments of order $ 3 $:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{array}{lcl}\n",
    "     \\mathbb{E}(W^{3}) &=&  \\left({n}^{3} + 3 \\, {n}^{2} + 4 \\, {n}\\right) \\Sigma^{3} + \\left(2 \\, {\\left({n}^{2} + {n}\\right)} {(\\mathrm{tr} \\, {\\Sigma})}\\right) \\Sigma^{2} + \\left({n} {(\\mathrm{tr} \\, {\\Sigma})}^{2} + {\\left({n}^{2} + {n}\\right)} {(\\mathrm{tr} \\, {\\Sigma}^{2})}\\right) \\Sigma \\\\\n",
    "\\mathbb{E}(W {(\\mathrm{tr} \\, W)}^{2}) &=& 8 \\, {n} \\Sigma^{3} + 4 \\, {n}^{2} {(\\mathrm{tr} \\, {\\Sigma})} \\Sigma^{2} + \\left({n}^{3} {(\\mathrm{tr} \\, {\\Sigma})}^{2} + 2 \\, {n}^{2} {(\\mathrm{tr} \\, {\\Sigma}^{2})}\\right) \\Sigma\n",
    "\\end{array}\n",
    "\\end{equation*}\n",
    "and\n",
    "\n",
    "\\begin{align*}\n",
    "\\mathbb{E}({{W^{-3}}})  &= \\frac{{\\left({n} - r - 1\\right)} {{\\Sigma}^{-3}}}{{\\left({n} - r + 1\\right)} {\\left({n} - r\\right)} {\\left({n} - r - 3\\right)} {\\left({n} - r - 5\\right)}} \\frac{2 \\, {{\\Sigma}^{-2}} {(\\mathrm{tr} \\, {\\Sigma}^{-1})}}{{\\left({n} - r + 1\\right)} {\\left({n} - r\\right)} {\\left({n} - r - 3\\right)} {\\left({n} - r - 5\\right)}} \\\\\n",
    "&\\phantom{=}+ \\frac{{\\left(2 \\, {(\\mathrm{tr} \\, {\\Sigma}^{-1})}^{2} + {n} {(\\mathrm{tr} \\, {\\Sigma}^{-2})} - r {(\\mathrm{tr} \\, {\\Sigma}^{-2})} - {(\\mathrm{tr} \\, {\\Sigma}^{-2})}\\right)} {{\\Sigma}^{-1}}}{{\\left({n} - r + 1\\right)} {\\left({n} - r\\right)} {\\left({n} - r - 1\\right)} {\\left({n} - r - 3\\right)} {\\left({n} - r - 5\\right)}}\n",
    "\\end{align*}\n",
    "\n",
    "After importing the Wishart moments package, we have to create an instance of the class `Expectations`.\n",
    "```sh\n",
    "sage: k=3\n",
    "sage: expec = wm.Expectations(k)\n",
    "```\n",
    "\n",
    "We need to know how to reference the expressions of which we want to compute their expectations. We can get a list the expressions of order $ k $ by using the method `expressions` of `Expectations`, which returns a list of 2-element lists with the index of the portrait and the the expression for expectation of the moment corresponding to it.\n",
    "```sh\n",
    "sage: expec.expressions()\n",
    "```\n",
    "```\n",
    "[0, W*tr(W, 1)^2]\n",
    "[1, 2/3*W^2*tr(W, 1) + 1/3*W*tr(W, 2)]\n",
    "[2, W^3]\n",
    "```\n",
    "\n",
    "Here `tr(A,j)` represents $ \\mathrm{tr}(A^j) $. Therefore, to get $ W {(\\mathrm{tr}\\,{W})}^{2} $ we call the method `moment` with the index `0`\n",
    "\n",
    "```\n",
    "```\n",
    "\n",
    "```sh\n",
    "sage: expec.moment(0)\n",
    "```\n",
    "\n",
    "```\n",
    "{\n",
    "\t'var': W*tr(W, 1)^2 ,\n",
    "\t'moment': 8*n*S^3 + 4*n^2*tr(S, 1)*S^2 + (n^3*tr(S, 1)^2 \n",
    "\t+ 2*n^2*tr(S, 2))*S\n",
    "}\n",
    "```\n",
    "\n",
    "Similarly we use the index `2` to get $ W^3 $.\n",
    "```sh\n",
    "sage: expec.moment(2)\n",
    "```\n",
    "```\n",
    "{\n",
    "\t'var': W^3 ,\n",
    "\t'moment': (n^3 + 3*n^2 + 4*n)*S^3 + (2*(n^2 + n)*tr(S, 1))*S^2 \n",
    "\t+ (n*tr(S, 1)^2 + (n^2 + n)*tr(S, 2))*S\n",
    "}\n",
    "```\n",
    "\n",
    "As for the moment of the inverse, we can call `moment` with the argument `inverse` set to `True`.\n",
    "```sh\n",
    "sage: expec.expressions(inverse = True)\n",
    "```\n",
    "```\n",
    "[0, inv(W, 1)*tr(W, -1)^2]\n",
    "[1, 2/3*inv(W, 2)*tr(W, -1) + 1/3*inv(W, 1)*tr(W, -2)]\n",
    "[2, inv(W, 3)]\n",
    "```\n",
    "\n",
    "Here `inv(A,j)` represents $ A^{-j} $. We use the index `2` to get $ W^{-3} $.\n",
    "```sh\n",
    "sage: expec.moment(2, inverse = True)\n",
    "```\n",
    "```\n",
    "{\n",
    "\t'var': inv(W, 3) ,\n",
    "\t'moment': (n - r - 1)*inv(S, 3)/((n - r + 1)*(n - r)\n",
    "\t  *(n - r - 3)*(n - r - 5)) \n",
    "\t+ 2*inv(S, 2)*tr(S, -1)/((n - r + 1)*(n - r)*(n - r - 3)*(n - r - 5)) \n",
    "\t+ (2*tr(S, -1)^2 + n*tr(S, -2) - r*tr(S, -2) \n",
    "\t- tr(S, -2))*inv(S, 1)/((n - r + 1)*(n - r)*(n - r - 1)\n",
    "\t  *(n - r - 3)*(n - r - 5))\n",
    "}\n",
    "```\n",
    "### Latex code of the results\n",
    "We can obtain the string with the $\\LaTeX$ code representing these expressions by using the built-in function `latex`. For instance, if we want to get the code for the variable $W {(\\mathrm{tr} \\, W)}^{2}$, we should use the following commands\n",
    "```sh\n",
    "sage: latex(expec.moment(0)['var'])\n",
    "```\n",
    "```\n",
    "W {(\\mathrm{tr} \\, W)}^{2}\n",
    "```\n",
    "and for its expectation, ${\\Sigma} {n}^{3} {(\\mathrm{tr} \\, {\\Sigma})}^{2} + 8 \\, {\\Sigma}^{3} {n}+ 2 \\, {\\left(2 \\, {\\Sigma}^{2} {(\\mathrm{tr} \\, {\\Sigma})}+ {\\Sigma} {(\\mathrm{tr} \\, {\\Sigma}^{2})}\\right)} {n}^{2}$,\n",
    "\n",
    "```sh\n",
    "sage: latex(expec.moment(0)['moment'])\n",
    "```\n",
    "```\n",
    "8 \\, {n} \\Sigma^{3} + 4 \\, {n}^{2} {(\\mathrm{tr} \\, {\\Sigma})} \\Sigma^{2}\n",
    "+ \\left({n}^{3} {(\\mathrm{tr} \\, {\\Sigma})}^{2} \n",
    "+ 2 \\, {n}^{2} {(\\mathrm{tr} \\, {\\Sigma}^{2})}\\right) \\Sigma\n",
    "```\n",
    "\n",
    "Notice that an instance of the form `wm.Expectations(3)` only permits to compute moments of order $3$. To compute a moment of a different order, say $k=4$, the user has to instantiate a new object of the class `wm.Expectations(4)`. We will continue with the examples using the same parameter as we were doing so far, that is `k=3`, so that can keep using the same object `expec`.\n",
    "\n",
    "### Evaluating the moments for specific values of the parameters\n",
    "\n",
    "Now we show how to compute the numerical value of the moment $E(W \\hbox{tr} (W^2)) $  for a Wishart distributions with parameters $n=10$ and \n",
    "\n",
    "\\begin{equation*}\n",
    "\\Sigma = \\begin{bmatrix}\n",
    "4 & 1\\\\\n",
    "1 & 3\n",
    "\\end{bmatrix}\n",
    "\\end{equation*}\n",
    "\n",
    ">We want to remark that neither the package nor the website will check if the matrix $\\Sigma$ is positive definite.\n",
    "\n",
    "We first set the matrix $\\Sigma$:\n",
    "\n",
    "```sh\n",
    "sage: Sigma = np.array([[4,1],[1,3]]);\n",
    "```\n",
    "\n",
    "To evaluate the moment we use the `evaluate_moment` where the parameters are `t` (the index of the expression in the list  `expec.expressions()` for which we require its expectation),  `n_param` (the numerical value for the parameter $n$), `Sigma` (the numerical value for the matrix $\\Sigma$) and the boolean parameter `inverse` (`False` will compute the moment of $W$ and `True` the moment for $W^{-1}$). Here we need to compute $E(W \\mathrm{tr} (W^2))$ and therefore\n",
    "the parameters are passed as follows.\n",
    "\n",
    "```sh\n",
    "sage: ev = expec.evaluate_moment(t=0, n_param=10, Sigma=Sigma, inverse=False);\n",
    "```\n",
    "\n",
    "As `ev` is a dictonary, we can retrieve the variable by using the key `'var'`\n",
    "```sh\n",
    "sage: ev['var']\n",
    "```\n",
    "```\n",
    "W*tr(W, 1)^2\n",
    "```\n",
    "\n",
    "To get the moment we use the key `'moment'`\n",
    " \n",
    "```sh\n",
    "sage: ev['moment']\n",
    "```\n",
    "```\n",
    "array([[813600., 231120.],\n",
    "       [231120., 582480.]], dtype=object)\n",
    "```\n",
    "\n",
    "## License\n",
    "\n",
    "GNU Public License v3\n",
    "\n",
    "[//]: # (These are reference links used in the body of this note and get stripped out when the markdown processor does its job. There is no need to format nicely because it shouldn't be seen. Thanks SO - http://stackoverflow.com/questions/4823468/store-comments-in-markdown-syntax)\n",
    "   [sage]: \n",
    "   [WM]: "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

Citation (CITATION.cff)

# This CITATION.cff file was generated with cffinit.
# Visit https://bit.ly/cffinit to generate yours today!

cff-version: 1.2.0
title: Wishart Moments Calculator
message: Please cite this software using these metadata.
type: software
authors:
  - given-names: Carlos Santiago
    family-names: Antunes Percíncula
    email: antunes.p.carlos@hotmail.com
    affiliation: Universidad Nacional del Litoral
repository-code: >-
  https://github.com/antunescarles/wishart-moments-calculator
url: >-
  https://antunescarles.github.io/wishart-moments-calculator/
abstract: >-
  A Python package using SAGE to compute all
  invariant Wishart moments and a web interface for
  those computations.
keywords:
  - wishart
  - moments
  - sage
  - python
  - web
  - interface

GitHub Events

Committers

Last synced: almost 3 years ago

Issues and Pull Requests

Last synced: about 2 years ago