{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Sobol Indices"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*Sobol's method* is one of the most popular for global sensitivity analysis. It builds on the [ANOVA decomposition](anova.ipynb)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"20D TT tensor:\n",
"\n",
" 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32\n",
" | | | | | | | | | | | | | | | | | | | |\n",
" (0) (1) (2) (3) (4) (5) (6) (7) (8) (9) (10)(11)(12)(13)(14)(15)(16)(17)(18)(19)\n",
" / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\ / \\\n",
"1 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 1"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import tntorch as tn\n",
"import torch\n",
"import time\n",
"\n",
"N = 20\n",
"t = tn.rand([32]*N, ranks_tt=10)\n",
"t"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With *tntorch* we can handle *all* Sobol indices (i.e. for all subsets $\\alpha \\subseteq \\{0, \\dots, N-1\\}$) at once. We can access and aggregate them using the function `sobol()` and whatever [mask](logic.ipynb) is appropriate."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Single Variables\n",
"\n",
"### Variance Components\n",
"\n",
"The relative influence (proportion of the overall model variance) attributable to one variable $n$ only, without interactions with others, is known as its *variance component* and denoted as $S_n$. Let's compute it for the first variable $x$:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"tensor(0.1882)\n",
"This compressed tensor has 58240 parameters; computing this index took only 0.0380359s\n"
]
}
],
"source": [
"x, y, z = tn.symbols(N)[:3]\n",
"start = time.time()\n",
"print(tn.sobol(t, mask=tn.only(x)))\n",
"print('This compressed tensor has {} parameters; computing this index took only {:g}s'.format(t.numcoef(), time.time()-start))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Input parameters $x, y, \\dots$ are assumed independently distributed. By default, uniform marginal distributions are used, but you can specify others with the `marginals` argument (list of vectors). For instance, if the first variable can take one value only, then its sensitivity indices will be 0 (no matter how strong its effect on the multidimensional model is!):"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"tensor(0.)"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"marginals = [None]*N # By default, None means uniform\n",
"marginals[0] = torch.zeros(t.shape[0])\n",
"marginals[0][0] = 1 # This marginal's PMF is all zeros but the first value\n",
"tn.sobol(t, tn.only(x), marginals=marginals)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Total Sobol Indices\n",
"\n",
"The effect that also includes $x$'s interaction with other variables is called *total Sobol index* (it's always larger or equal than the corresponding variance component):"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"tensor(0.2150)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tn.sobol(t, x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Tuples of variables\n",
"\n",
"What are the indices for the first and third variables $x$ and $z$?"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"tensor(0.0005)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tn.sobol(t, tn.only(x & z)) # Variance component"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"tensor(0.2401)"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tn.sobol(t, x | z) # Total index"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What's the relative importance of $x$ with respect to the group $\\{y, z\\}$?"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"tensor(0.1638)"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tn.sobol(t, x & (y|z)) / tn.sobol(t, y|z)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Closed Sobol Indices\n",
"\n",
"For tuples of variables two additional kinds of indices exist. The *closed index* aggregates all components for tuples *included* in $\\alpha$, and for tuple $\\{x, z\\}$ it can be computed as follows:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"tensor(0.2100)"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tn.sobol(t, tn.only(x | z))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Superset Indices\n",
"\n",
"The *superset index* aggregates all components for tuples *that include* $\\alpha$:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"tensor(0.0009)"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tn.sobol(t, x & z)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Counting $k$-plets of Variables\n",
"\n",
"We can also easily count the influence of all $k$-plets of variables combined:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"tensor(0.9222)"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tn.sobol(t, tn.weight_mask(N, weight=[1]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Often, there are different ways to express the same mask. For example, these three are equivalent:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"tensor(0.2401)\n",
"tensor(0.2401)\n",
"tensor(0.2401)\n"
]
}
],
"source": [
"print(tn.sobol(t, x | z))\n",
"print(tn.sobol(t, x & ~z) + tn.sobol(t, ~x & z) + tn.sobol(t, x & z))\n",
"print(tn.sobol(t, x) + tn.sobol(t, z) - tn.sobol(t, x & z))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The Mean Dimension\n",
"\n",
"Variance components are the basis for an important advanced sensitivity metric, the [*mean dimension*](https://www.jstor.org/stable/27590729). It's defined as $D_S := \\sum_{\\alpha} |\\alpha| \\cdot S_{\\alpha}$ and computed as:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"tensor(1.0831)"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tn.mean_dimension(t)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can also compute it in one line by weighting the Sobol indices by their tuple weight (according to the definition of mean dimension):"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"tensor(1.0831)"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tn.sobol(t, tn.weight(N))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The mean dimension is always greater or equal than 1. It gives a notion of *complexity* of a multidimensional function (the lower the mean dimension, the simpler it is). For example, rounding a tensor usually results in a lower mean dimension:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEKCAYAAAA4t9PUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VdW5//HPk5OJIcyBRAIiiiAySYKi0jpUK1hvsS22VStqtdReawdvW+vt7a/z7WB77fVatVocaC21itYJUdtarVZQQAQEmRUiAYJAIEwhyfP7Y+/gMWY4IeecneH7fr3265y99nCeHDZ5stZeey1zd0RERI5URtQBiIhI+6ZEIiIiraJEIiIiraJEIiIiraJEIiIiraJEIiIiraJEIiIiraJEIiIiraJEIiIirZIZdQDp0K9fPx8yZEjUYYiItCuLFi3a7u75ze3XKRLJkCFDWLhwYdRhiIi0K2b2diL7qWlLRERaRYlERERaRYlERERaRYlERERaRYlERERaRYlERERaRYlERERapVM8R3KkHl5cylvb94IZBphBRtx7MwtesXDbe+8t3C8jfr/DZXXHvfc+w4yMjLp1I2Z2+NiMcHss4719YxnBknn4NSN4jb2/PBa3PTNmZGVkkBWrO5dF+v2KSMegRNKEJ5aW8dyqbXTUae2zYkZWLIPMDCM7M+NwssmOhUknlkF2ZgbZ4WtOZoyczLr3723LycogOxYjOzODLlkZdM3JpHtOJl2zY+FruJ4TrOdkZiiJiXQgSiRNuPuKCYffuzvu4EDt4ffha9z7WnecsMwbKYs7x+Fzhe/fW8L1Wj5Q5u5U1zg17tTUOtW1Tk1N+FrrVNfWvld+eHst1bXOoRrnUE0t1TW1HKp1DlXXldeG5X64/FBNLVU1tRysrmXPgWrera7iYHUNVTW1VFUHy8Hwtbo28Wwby7C4JBMjw4wad2prPXzlcNy14c9Yt62mNvi+iKsB1tXYqKvZxdXk6mqBWbEM8nKDhJaXm0n33CzycjPJq1vPySQvN4vuuZlheRb98rLJ755DZkwtwCJNUSJJUF3zFEAM/TVdX02YjPZX1bC3qpq9B2uoPFjNvqpq9h4M1vdWVQdl79tWQ607GRlBc14sw8iwoGkuI8OIZRA088Vtr6vNeJhggyQcJFl4L/HWlbs7VTW1VB6oZs+BarZXVrFh+14qD1az+0A1VdW1jf5cGQb983IZ0DOXwh65FPTMpbBn3WsXCnvm0r9HDjmZsbR8zyJtkRKJJEVwLyZGblaM3t2yow6nRQ5W11B5IEhyew7ULYfYXlnFlor9lFUcYMvuA6wtr+TFtdupPFj9gXP0655NQc9cBvbqwqDeXSnq3YVBfbpSFL7vlqP/atJx6eqWTi8nM0ZO9xh9u+cktP+eA4fYuvsAZRXBsqWi7v1+1pXv5fnV5Rw49P5aTp9u2Qzq3SVILH2C10G9u3BMv24M7tNV94ykXVMiEWmhvNws8nKzOK5/XoPb3Z3tlVWU7tzHpp37g9cdweuKst08u2IrVTXvJZpeXbMYW9SLcYOCZeygXvRpZ7U66dyUSESSzMzIz8shPy+Hkwb3/sD22lpn256DlO7cx5ptlSzZuIslm3Zxy5o1h3sIHt23a5BUinoxbnAvRhb2IDdL92GkbTLvqH1b45SUlLjmI5G2rvJgNctKK1iyaRdLNu1kyaZdbN19EAi6ao8s7MHYQb04c3g+Hx6Wr95kknJmtsjdS5rdT4lEpO3aUnGAJZt28tqmXSzZuItl71Swr6qGAT1ymFZcxKdLBnF0325RhykdlBJJHCUS6SgO1dTyt5XbeODVjTy/upxah4lD+/DZCYOZPKpAzV+SVEokcZRIpCMqq9jPnEWl/HlhKRt37KNHbiZTxw3kMxMGMWpgz6jDkw5AiSSOEol0ZLW1zvwN7/LAq5t4avkWqqprOfGoHnxmwiCmjh1Iz65ZUYco7ZQSSRwlEuksKvYd4tHX3+FPr2xiRdlucjIzmDKqgC+ecSwnFPaIOjxpZ5RI4iiRSGe0/J0KHnh1E3957R0qq6r5+NijuP7c43VzXhKmRBJHiUQ6s4p9h7jjhXXc89IGqmucT08YxFc/MowBPXKjDk3aOCWSOEokIrBt9wFufW4ts1/ZSIYZV5w2hGvOOLbdjY0m6aNEEkeJROQ9G9/dx6//uppHlrxD9+xMZnx4KJ+fdIwGlpQPUCKJo0Qi8kGrtuzhl8+s4tkVW+nXPZtrzzqOS04ZrCHx5TAlkjhKJCKNW7xxJzfNW8XL699lYK8ufPWcYXzypIEagkUSTiQpvVLM7G4z22ZmyxvZbmZ2i5mtNbOlZjY+btsvzOwNM1sZ7mNhebGZLQuPOVwuIkdm/ODe/PELp/CHq06hb/dsvvXQUj52y4us3ron6tCknUj1nxz3ApOb2D4FGBYuM4DbAczsNOB0YAwwCpgAnBEec3u4b91xTZ1fRBJgZkwa1o9Hrz2d2y8dz7t7q5h660s8uuSdqEOTdiClicTdXwB2NLHLVGCWB+YDvcyskGBq9FwgG8gBsoCt4bYe7v6yB21ys4ALU/kziHQmZsaU0YU8+ZVJjBrYg6/+aQnfe3R5k9MRi0TdCDoQ2BS3XgoMdPeXgeeAsnB52t1XhvuX1t8/TbGKdBoDeuTyxy9M5OpJx3Dfy2/zmTtfpqxif9RhSRsVdSJp6P6Gm9lxwAlAEUGiONvMPtzY/g2e2GyGmS00s4Xl5eVJC1iks8iKZfBfF4zktkvHs3rLHj52y4u8tHZ71GFJGxR1IikFBsWtFwGbgU8A89290t0rgaeAieH+RQ3s/wHufqe7l7h7SX5+fkqCF+kMzh9dyGPXTaJvt2wum7mA3zy3ltrajt/bUxIXdSJ5DJge9t6aCFS4exmwETjDzDLNLIvgRvvKcNseM5sY9taaDjwaWfQincSx+d35y7Wnc8GYo7jp6VV8YdZCKvYdijosaSNS3f13NvAyMNzMSs3sKjO7xsyuCXeZC6wH1gJ3Af8elj8ErAOWAa8Dr7v74+G2LwG/C49ZR1BbEZEU65aTyf9+dhw/nHoiL6wp54Jb/8nydyqiDkvaAD2QKCIttnjjTq69fzHv7q3ix1NH8ekJg5o/SNqdNvFAooh0TOMH9+aJ6yZx8pA+fGvOUr49ZykHDtVEHZZERIlERI5I3+453Pf5k/nyWcfxp1c3ccld89lzQPdNOiMlEhE5YrEM4xvnDee2S8eztLSC6Xe/wm4lk05HiUREWu380YXcesl4lpVWMH2mkklno0QiIkkxeVQBt106njc2V3DZzFeo2K9k0lkokYhI0nz0xAJuu7SYFZsrmD5zgZJJJ6FEIiJJde7IAdx+aTErynZz2cwFenCxE1AiEZGkO2fkAO74XDFvlu3hc0omHZ4SiYikxEdOGMBvLytm1ZY9XDpzPrv2VUUdkqSIEomIpMxZI/rz2+nFrN5aySV3LWDnXiWTjkiJRERS6qzh/bnzsmLWlldy6e+UTDoiJRIRSbkzh/fnruklrC2v5JLfLWCHkkmHokQiImlxxvH5/G56CevLK7nkrvlKJh2IEomIpM2Hj89n5uUT2LB9L5fcNV/NXB2EEomIpNWkYf24+4oJrCuv5Ofz3ow6HEkCJRIRSbvTj+vH9FOH8OeFm3hzy+6ow5FWUiIRkUhcd/ZxdM/J5KdzVStp75RIRCQSvbpm85WPDOP51eW8sLo86nCkFZRIRCQyl516NIP6dOG/566kprbjT/vdUSmRiEhkcjJj3DB5BG9u2cOcxaVRhyNHSIlERCL1sdGFjBvUi189s4p9VdVRhyNHQIlERCJlZnz3ghPYuvsgv/vnhqjDkSOgRCIikSs+ug/njy7gjufXsW33gajDkRZSIhGRNuFb543gUE0tN/91ddShSAspkYhImzCkXzcumziEB17dxKote6IOR1pAiURE2ozDDyk+tTLqUKQFlEhEpM3o3S2b684exj9WlfPPNXpIsb1IWSIxs7vNbJuZLW9ku5nZLWa21syWmtn4sPwsM1sStxwwswvDbfea2Ya4beNSFb+IRGP6aUdT1LsLP3lSDym2F6mskdwLTG5i+xRgWLjMAG4HcPfn3H2cu48Dzgb2Ac/EHffNuu3uviQlkYtIZOIfUnxYDym2CylLJO7+ArCjiV2mArM8MB/oZWaF9faZBjzl7vtSFaeItD0XjAkeUvzlM6vYX1UTdTjSjCjvkQwENsWtl4Zl8T4LzK5X9pOwKexmM8tJZYAiEg0z4zsfq3tIcX3U4Ugzokwk1kDZ4QbRsHYyGng6bvuNwAhgAtAHuKHRk5vNMLOFZrawvFw37UTamwlD+jD5xAJuf34d2/boIcW2LKFEYmanmdklZja9bknCZ5cCg+LWi4DNceufBh5x90N1Be5eFjaFHQTuAU5u7OTufqe7l7h7SX5+fhLCFZF0u2HKCKqqa7n52TVRhyJNaDaRmNnvgV8CkwhqAhOAkiR89mPA9LD31kSgwt3L4rZfTL1mrbp7KGZmwIVAgz3CRKRjOKZfNz438WgeeHUjq7fqIcW2KjOBfUqAke7eon54ZjYbOBPoZ2alwPeALAB3vwOYC5wPrCXomXVl3LFDCGorz9c77f1mlk/QLLYEuKYlMYlI+/OVjwxjzuJSfjp3Jfdc2WgjhEQokUSyHCgAyprbMZ67X9zMdgeubWTbW3zwxjvufnZLYhCR9q9Pt2y+fNZx/PSpN3lxzXYmDesXdUhSTyL3SPoBK8zsaTN7rG5JdWAiInUuP20IA3t14SdzV1KrhxTbnERqJN9PdRAiIk3JzYrx9XOP5xsPvs78De9y2rGqlbQlzdZI3P154E0gL1xWhmUiImnzsdGF5OVk8tBCPe3e1iTSa+vTwCvARQRdcheY2bRUByYiEq9LdowLxhYyd3kZew4cav4ASZtE7pF8B5jg7pe7+3SCZze+m9qwREQ+aFrxIA4cqmXushb1/ZEUSySRZLj7trj1dxM8TkQkqcYP7sXQ/G48tEjNW21JIglhXthj6wozuwJ4kuAZEBGRtDIzphUX8epbO9mwfW/U4UgokZvt3wTuBMYAY4E73b3RMa5ERFLpkycVkWEwR7WSNiOR7r+4+xxgTopjERFpVkHPXD40LJ85i0v5+rnHE8toaPxXSadGayRm9mL4usfMdscte8xsd/pCFBF5v2nFRZRVHOBf67ZHHYrQRCJx90nha56794hb8ty9R/pCFBF5v3NHDqBHbiYP6pmSNiGR50iOrZtAyszONLOvmFmv1IcmItKw3KwYHx93FE+/sYWK/XqmJGqJ9NqaA9SY2XHATOAY4I8pjUpEpBkXFQ/iYHUtTy7VMyVRSySR1Lp7NfAJ4Nfu/nWg/tzqIiJpNaaoJ8P6d+fBRZua31lSKpFEcsjMLgYuB54Iy7JSF5KISPPMjItKinht4y7WbquMOpxOLZFEciVwKvATd99gZscAf0htWCIizbvwpIHEMkxPukcskQcSV7j7V9x9dri+wd1/lvrQRESa1j8vlzOPz+eR10qp0TwlkUmk19bpZvasma02s/VmtsHM1qcjOBGR5kwrLmLr7oO8sKY86lA6rUSebJ8JfB1YBNSkNhwRkZb5yAkD6N01i4cWlXLW8P5Rh9MpJZJIKtz9qZRHIiJyBLIzM5g6biB/XLCRXfuq6NU1O+qQOp1EbrY/Z2Y3mdmpZja+bkl5ZCIiCZpWXERVTS2Pv7456lA6pURqJKeEryVxZQ6cnfxwRERa7sSjejCiII8HF5Vy2alDog6n02k2kbj7WekIRETkSAXPlAziR0+sYNWWPQwvyIs6pE4lkV5bA8xsppk9Fa6PNLOrUh+aiEjiLhx3FJkZxkN60j3tErlHci/wNHBUuL4a+FqqAhIRORJ9u+dw9oj+PPLaZg7V1EYdTqeSSCLp5+5/BmoBwnG31A1YRNqcacVFbK88yPOr9ExJOiWSSPaaWV+CG+yY2USgIqVRiYgcgbNG9Kdvt2wNmZJmiSSS64HHgGPN7CVgFnBdcweZ2d1mts3Mljey3czsFjNba2ZL67oUm9lZZrYkbjlgZheG244xswVmtsbMHjAzdRgXkcOyYhlceNJA/vbmVnbsrYo6nE4jkbG2FgNnAKcBXwROdPelCZz7XmByE9unAMPCZQZwe/h5z7n7OHcfR9DFeB/wTHjMz4Gb3X0YsBPQTX8ReZ9pxUUcqnEeXfJO1KF0Gon02ooB5wMfAT4KXGdm1zd3nLu/AOxoYpepwCwPzAd6mVn9eU6mAU+5+z4zM4LE8lC47T7gwubiEJHO5YTCHowa2EPT8KZRIk1bjwNXAH2BvLiltQYC8f30SsOyeJ8FZofv+wK7wpv9je0vIsJFxYNYUbabNzbrdm46JPJke5G7j0nBZ1sDZYfHgQ5rJ6MJuh43u/8HTm42g6DJjMGDBx95lCLS7nx87FH85MmVPLSolBOP6hl1OB1eIjWSp8zsoyn47FJgUNx6ERA/UM6ngUfc/VC4vp2g+Suzkf3fx93vdPcSdy/Jz89PYtgi0tb17pbNOSP78+iSzVRV65mSVEskkcwHHjGz/Wa228z2mNnuJHz2Y8D0sPfWRIJRhsvitl/Me81auLsDzxHcN4Fg6t9HkxCHiHRA04qL2LG3ir+/uS3qUDq8RBLJrwim2u3q7j3cPc/dezR3kJnNBl4GhptZqZldZWbXmNk14S5zgfXAWuAu4N/jjh1CUFt5vt5pbwCuN7O1BPdMZiYQv4h0Qh8elk9+Xg5zFuume6olco9kDbA8rBEkzN0vbma7A9c2su0tGriR7u7rgZNbEoeIdE6ZsQwuGFPI/Qs2Unmwmu45ify6kyORSI2kDPiHmd1oZtfXLakOTESktaaMKqSqulbNWymWSCLZAPwNyCa53X9FRFKq+Oje9Ouew7zlZc3vLEcskflIfpCOQEREki2WYXz0xAE8svgdDhyqITcrFnVIHVKjNRIz+3X4+riZPVZ/SV+IIiJHbsqoAvYfquH51RoROFWaqpH8Pnz9ZToCERFJhYlD+9KzSxbzlm/hvBMLog6nQ2o0kbj7ovC1fhdcEZF2IyuWwbkjB/D0G1uoqq4lOzORW8PSEo0mEjNbRhNDkKRo2BQRkaSbMqqAhxaV8tK67Zw1vH/U4XQ4TTVtXRC+1j3rUdfUdSnB0O4iIu3C6cf1o3tOJvOWbVEiSYFG63ju/ra7vw2c7u7fcvdl4fJt4Lz0hSgi0jq5WTHOGtGfZ1dupVrzuSddIo2F3cxsUt2KmZ0GdEtdSCIiyTdlVAE79lbxyltNTZMkRyKRMQOuAu42s54E90wqgM+nNCoRkSQ7c3g+uVkZzFu+hdOO7Rd1OB1KIlPtLnL3scAYYFw4De7i1IcmIpI8XbMzOeP4fOYt30JtbYuGDpRmJNwPzt13u7umGxORdmvKqEK27TnIa5t2Rh1Kh6IO1SLSaZx9Qn+yYsZTy7ZEHUqHokQiIp1Gj9wsJh3Xj3lvbKGFM2NIExIaoD/sqTUkfn93n5WimEREUmbyqAKem7OMNzbvZtRAzeeeDM3WSMzs9wTjbU0CJoRLSYrjEhFJiXNHFhDLMJ7S0PJJk0iNpAQY2dIZEkVE2qI+3bI55Zg+PLV8C9/46HDMLOqQ2r1E7pEsBzRkpoh0GFNGFbC+fC9rtlVGHUqHkEgi6QesMLOnNR+JiHQE551YgBnqvZUkiTRtfT/VQYiIpFP/HrkUD+7NU8vL+Oo5w6IOp91LZKpdzUciIh3O5FEF/PjJlby1fS9D+mn4wNZIpNfWRDN71cwqzazKzGrMbHc6ghMRSZW62RLnvaHmrdZK5B7JrcDFwBqgC3B1WCYi0m4N6tOV0QN78tRyJZLWSujJdndfC8Tcvcbd7wHOTGlUIiJpMHlUAa9v2sXmXfujDqVdSySR7DOzbGCJmf3CzL6O5iMRkQ5gyqiweUu1klZJJJFcFu73ZWAvMAj4VCqDEhFJh6H53Rk+IE+JpJUSmY/kbcCAQnf/gbtfHzZ1iYi0e+eNKuDVt3dQvudg1KG0W4n02vo3YAkwL1wfl8gDiWZ2t5ltM7PljWw3M7vFzNaa2VIzGx+3bbCZPWNmK81shZkNCcvvNbMNZrYkXMYl9mOKiDRsyqgC3OGZFaqVHKlEmra+D5wM7AJw9yUEIwE3515gchPbpwDDwmUGcHvctlnATe5+QvjZ2+K2fTOcpXFcGIuIyBEbUZDHkL5d1bzVCokkkuojmRnR3V8AdjSxy1RglgfmA73MrNDMRgKZ7v5seJ5Kd9/X0s8XEUmEmTF5VCEvr3uXXfuqog6nXUpo0EYzuwSImdkwM/s/4F9J+OyBwKa49dKw7Hhgl5k9bGavmdlNZhaL2+8nYVPYzWaW09jJzWyGmS00s4Xl5eVJCFdEOqopowqornWeXbE16lDapUQSyXXAicBBYDawG/haEj67obGbnWDYlg8B3yCY+2QocEW4/UZgRFjeB7ihsZO7+53uXuLuJfn5+UkIV0Q6qjFFPTmqZ66at45QIr229rn7d9x9QviL+TvufiAJn11K0JW4ThGwOSx/zd3Xu3s18BdgfBhLWdgUdhC4h+D+iYhIq5gZ540q4J9rtlN5sDrqcNqdRgdtbK5nlrt/vJWf/RjwZTP7E3AKUOHuZWa2DehtZvnuXg6cDSwMYyoM9zHgQoK5UkREWm3KqELueekt/v7mNj4+9qiow2lXmhr991SCexizgQU03BTVKDObTTCUSj8zKwW+B2QBuPsdwFzgfGAtsA+4MtxWY2bfAP4WJoxFwF3hae83s/wwliXANS2JSUSkMcVH96Zf9xzmLS9TImmhphJJAXAuwYCNlwBPArPd/Y1ETuzuFzez3YFrG9n2LDCmgfKzE/lsEZGWimUY5504gIcXv8P+qhq6ZMeaP0iAJu6RhAM0znP3y4GJBDWHf5jZdWmLTkQkjaaMKmT/oRrmvVEWdSjtSpMTW4Xdaz9GUCsZAtwCPJz6sERE0m/i0D6MLOzBd//yBiMLezK8IC/qkNqFRmskZnYfwfMi44EfhL22fuTu76QtOhGRNMqMZTDzihK6Zsf4/L2vsm1PMjqodnxNdf+9jODhwK8C/zKz3eGyRzMkikhHVdizCzMvn8COvVV8YdYiDhyqiTqkNq+peyQZ7p4XLj3iljx375HOIEVE0ml0UU9+/dlxLC3dxX/8+XVqaz3qkNq0hGZIFBHpbM47sYAbp4zgyWVl/M+zq6MOp01r8ma7iEhn9oUPDWXD9r3c+txahvTrxrTioqhDapOUSEREGmFm/HDqKDbu2MeNDy+lqHcXJg7tG3VYbY6atkREmpAVy+C2S4oZ3KcrX/z9ItaXV0YdUpujRCIi0oyeXbO4+4oJxDKMq+5bqHlL6lEiERFJwNF9u3HnZcW8s3M/X/z9Iqqqa6MOqc1QIhERSVDJkD78YtoYFmzYwY0PLyMYMlB0s11EpAUuPGkgG7bv5X//toah+d249qzjog4pckokIiIt9LVzhrFh+15uenoVQ/p242NjCqMOKVJq2hIRaSEz4xfTxlB8dG+u//MSXtu4M+qQIqVEIiJyBHKzYtx5WTH9e+TwhVkLeX3TrqhDiowSiYjIEerbPYd7rjiZnMwYF/32ZeYsKo06pEgokYiItMJx/bvz+HWTKB7cm/948HV+8PgbHKrpXF2DlUhERFqpT7dsZl11MleePoR7XnqL6TNfYcfezvPQohKJiEgSZMUy+N6/ncgvLxrLoo07+bf/e5E3NldEHVZaKJGIiCTRtOIiHvziqdTUOp+6/V889vrmqENKOSUSEZEkGzuoF49fN4nRA3vyldmv8dOnVlLTgSfHUiIREUmB/Lwc7r96Ip+bOJjfPr+eK+55pcMO9qhEIiKSItmZGfz4wtH89JOjmb/+Xab+5iVWbdkTdVhJp0QiIpJiF588mD/NmMi+qho+cdtLzFteFnVISaVEIiKSBsVH9+GJ6yZx/IA8rvnDYn7z3NqoQ0oaJRIRkTQZ0COXB744kanjjuKmp1fx8OKO8SR8yhKJmd1tZtvMbHkj283MbjGztWa21MzGx20bbGbPmNlKM1thZkPC8mPMbIGZrTGzB8wsO1Xxi4ikQk5mjF9eNJaJQ/vw7YeXdYgBH1NZI7kXmNzE9inAsHCZAdwet20WcJO7nwCcDGwLy38O3Ozuw4CdwFVJjllEJOWyYhncdmkxA3rkMOP3iyir2B91SK2SskTi7i8AO5rYZSowywPzgV5mVmhmI4FMd382PE+lu+8zMwPOBh4Kj78PuDBV8YuIpFKfbtnMvHwC+w5WM2PWIvZX1UQd0hGL8h7JQGBT3HppWHY8sMvMHjaz18zsJjOLAX2BXe5eXW9/EZF26fgBefzvZ09i+eYKvjVnabudujfKRGINlDnBrI0fAr4BTACGAlc0sX/DJzebYWYLzWxheXl566MVEUmBc0YO4JvnDefx1ze3255cUSaSUmBQ3HoRsDksf83d14e1j78A44HtBM1fmfX2b5C73+nuJe5ekp+fn5IfQEQkGb50xrFcOO4ofvnMap5+Y0vU4bRYlInkMWB62HtrIlDh7mXAq0BvM6v77X82sMKDOt9zwLSw/HLg0XQHLSKSbGbGzz41hrFFPfn6A0tYWbY76pBaJJXdf2cDLwPDzazUzK4ys2vM7Jpwl7nAemAtcBfw7wDuXkPQrPU3M1tG0KR1V3jMDcD1ZraW4J7JzFTFLyKSTrlZMe6cXkJebiZX37eQdysPRh1Swqy93txpiZKSEl+4cGHUYYiINGtp6S4uuuNlxhb14g9Xn0J2ZnQNR2a2yN1LmttPT7aLiLQhY4p68YtpY3jlrR38v0eXt4ueXJnN7yIiIuk0ddxAVm/dw2+eW8fwgjyuPP2YqENqkmokIiJt0H+cO5xzRw7gR0+s4J9r2vYjDEokIiJtUEaGcfNnxjGsfx7X3r+Y9eWVUYfUKCUSEZE2qntOJr+7vITMWAZXz1pIxf5DUYfUICUSEZE2bFCfrtx+6Xg2vruP7z3a4GDqkVMiERFp404Z2pdrzzqOvyzZzN/f3Bp1OB+gRCIi0g5ce9ZxDB+Qx38+vJzdB9qp17rrAAAMcUlEQVRWE5cSiYhIO5CdmcEvpo1h254D/HTum1GH8z5KJCIi7cTYQb24+kNDmf3KRv61bnvU4RymRCIi0o58/ZzjGdK3Kzc+vKzNTIalRCIi0o50yY7x80+N4e139/GrZ1ZFHQ6gRCIi0u6cMrQvn5s4mLtf2sBrG3dGHY4SiYhIe3TD5BEU9MjlWw8t5WB1tE1cSiQiIu1QXm4W//3J0azZVslv/h7tFL1KJCIi7dSZw/vzyfEDue0f61ixObpZFZVIRETase9+bCS9umZxw5ylVNfURhKDEomISDvWu1s2P5w6imXvVHDXPzdEEoMSiYhIO3f+6EImn1jAzX9dzboIhptXIhER6QB+eOGJdMmK8e05S6mtTe/0vEokIiIdQP+8XL57wUhefWsnv5//dlo/W4lERKSD+NT4gZxxfD4/n/cmm3bsS9vnKpGIiHQQZsZ/f3I0BvznI8twT08TlxKJiEgHMrBXF749ZQT/XLOdBxeVpuUzlUhERDqYS085mpOH9OHHT6xg2+4DKf88JRIRkQ4mI8P42adGM/7o3lSl4SHFzJR/goiIpN3Q/O7ce+XJafmslNVIzOxuM9tmZssb2W5mdouZrTWzpWY2Pm5bjZktCZfH4srvNbMNcdvGpSp+ERFJTCprJPcCtwKzGtk+BRgWLqcAt4evAPvdvbEk8U13fyiJcYqISCukrEbi7i8AO5rYZSowywPzgV5mVpiqeEREJDWivNk+ENgUt14algHkmtlCM5tvZhfWO+4nYVPYzWaWk5ZIRUSkUVEmEmugrO7pmcHuXgJcAvzazI4Ny28ERgATgD7ADY2e3GxGmIwWlpeXJzFsERGJF2UiKQUGxa0XAZsB3L3udT3wD+CkcL0sbAo7CNwDNNolwd3vdPcSdy/Jz89PzU8gIiKRJpLHgOlh762JQIW7l5lZ77omKzPrB5wOrAjXC8NXAy4EGuwRJiIi6ZOyXltmNhs4E+hnZqXA94AsAHe/A5gLnA+sBfYBV4aHngD81sxqCRLdz9x9RbjtfjPLJ2gWWwJck6r4RUQkMZauQb2iZGblwJGOq9wP2J7EcJJFcbWM4moZxdUyHTWuo9292XsDnSKRtIaZLQxv/LcpiqtlFFfLKK6W6exxaawtERFpFSUSERFpFSWS5t0ZdQCNUFwto7haRnG1TKeOS/dIRESkVVQjERGRVlEiAczsR+H4XUvM7BkzO6qR/S43szXhcnlcebGZLQuHxL8lfGAyGXHdZGZvhrE9Yma9GthneNyw+kvMbLeZfS3c9n0zeydu2/npiivc763we1liZgvjyvuY2bPh9/ismfVOV1xmNsjMnjOzlWb2hpl9NW5b1N/XZDNbFV5H344rP8bMFoTf1wNmlp2kuC4Kv4NaM2uwZ09E11ezcYX7pfv6SuT7iuL6SvT7St315e6dfgF6xL3/CnBHA/v0AdaHr73D973Dba8ApxI8KPkUMCVJcX0UyAzf/xz4eTP7x4AtBH2/Ab4PfCMF31dCcQFvAf0aKP8F8O3w/beb+7mSGRdQCIwP3+cBq4GRUX9f4b/dOmAokA28HhfXn4HPhu/vAL6UpLhOAIYTDENUksD+6bq+Eoorguur2bgiur4SiSul15dqJIC7745b7cZ7g0fGOw941t13uPtO4FlgsgXDtvRw95c9+JeYRTB8SzLiesbdq8PV+QTjkTXlI8A6dz/Shy9TFVd9U4H7wvf3kcbvy4Px2haH7/cAK3lv1OmUSPD7OhlY6+7r3b0K+BMwNazdng3UzcGTzO9rpbuvasEh6bq+WhpXfam6vpqNK6LrK5HvK6XXlxJJyMx+YmabgEuB/9fALo0Nez8wfF+/PNk+T1Dbacpngdn1yr4cNqncnawqfgvicuAZM1tkZjPiyge4exkE//GA/mmOCwAzG0IwIOiCuOKovq/Grq++wK64RJSq6ysRUVxfTYny+mpWRNdXY1J6fXWaRGJmfzWz5Q0sUwHc/TvuPgi4H/hyQ6dooMybKE9KXOE+3wGqw9gaO0828HHgwbji24FjgXFAGfCrNMd1uruPJ5gN81oz+3Cin5/iuDCz7sAc4GtxNdIov6/Irq8Ez5P26ysBkVxfCZ4n7ddXc6dooKzV11edVE6126a4+zkJ7vpH4EmCQSbjlRIMQlmniKBNspT3N1UcHg4/GXFZcFP/AuAjYdNZY6YAi919a9y5D783s7uAJ9IZl783HcA2M3uEoHr9ArDVzAo9GO25ENiWzrjMLIvgP/n97v5w3Lmj/L4am1ZhO8HsoZnhX41Jvb5aIK3XV4LnSPv1lYgorq8EpOT6qtNpaiRNMbNhcasfB95sYLengY9aMMx9b4IbqE+HVec9ZjYxbG+cDjyapLgmE0ze9XF339fM7hdTr9nB3j918SdI0rD7icRlZt3MLK/uPcH3Vff5jwF1vd4uJ43fV/hvNBNY6e7/U29bZN8X8CowzIIeNNkEzUiPhUnnOWBauF/Svq8WStv1lYgorq8E40r79ZWg1F5fyew90F4Xgr8elgNLgceBgWF5CfC7uP0+TzDs/VrgyrjykvD4dcCthA96JiGutQTtmkvC5Y6w/Chgbtx+XYF3gZ71jv89sCz8uR4DCtMVF0HvkNfD5Q3gO3HH9wX+BqwJX/ukMa5JBFX3pXH7nR/19xWun0/Qy2ddve9rKEHPwLUETUs5SYrrEwR/qR4EthL8YdQWrq9m44ro+kokriiur0T/HVN2fenJdhERaRU1bYmISKsokYiISKsokYiISKsokYiISKsokYiISKsokUi7YWaVrTj2yxaMeupm1i+u3CwYsXmtBUNXjI/bVmhmCT80liwWjGrbL3z/r3R//pEysz/VeyZLOgklEuksXgLOAeoPODgFGBYuMwiGsahzPXBX/ROZWdpGhHD309L1WWYWq7duZpbQ74jw2NuBb6UiNmnblEik3Ql/wd0UjjW0zMw+E5ZnmNltFszN8ISZzTWzaQDu/pq7v9XA6aYCszwwn2C4iLonkD8FzAvPfYWZPWhmjxMMFNhYDGfG12LM7FYzuyJ8/5aZ/cDMFofHjAjL+1owD85rZvZb4sY/qquFhef9h5k9ZMHcJveHT1FjZueHZS+GtasP1KLMLBbG+2pY8/pi3HmfM7M/AsvMbIgFc2ncBiwGBpnZxWG8y83s5/GxmdkPzWwBwTQK/wTOSWeilbZBiUTao08SDHw3lqCWcVP4y/+TwBBgNHA1wS+35jQ4KqqZHQPsdPeDcdtOBS5397ObiKE52z0YaPB24Bth2feAF939JIInngc3cuxJwNeAkQRPI59uZrnAbwnmwJkE5Ddy7FVAhbtPACYAXwh/RgjGqPqOu48M14cTJNeTgEMEc6icHf68E8ysbpjxbsBydz/F3V9091qCp6PHJvA9SAeiRCLt0SRgtrvXeDAQ3vMEvxwnAQ+6e627byEYQ6g5jY1+WgiU1yt/1t13NBNDc+oG8VtEkPQAPgz8AcDdnwR2NnLsK+5eGv7CXhIePwJY7+4bwn3qD/Ne56PAdDNbQjCseV+C5ry6826I2/ftsHZG+DP9w93LPRjU7/4wXoAaguGF4m0jGJpDOhFVQaU9amwq4yOZ4rixUVELgNx6++5N4LOqef8faPXPUVfDqeH9//8SGasovnZUd3yiP7MB17n70+8rNDuT9/9ckNjPCXDA3WvqleUC+xOMSToI1UikPXoB+EzY7p9P8BfyK8CLwKfCeyUDeP+w/415jOAvdTOziQTNP2UEg9sNOYIY3gZGmlmOmfUkmFUwkZ/nUgAzm0IwlXOi3gSGWjCJEsBnGtnvaeBLFgxxjpkdb8Gouc1ZAJxhZv3CG+oXE9S+GnM8wSCK0omoRiLt0SME9yteJ/hL/lvuvsXM5hD84l5OkAgWABUAZvYVgh5FBcBSM5vr7lcDcwlGRV0L7AOuBHD3vWa2zsyOc/e1icYQftafCUZ4XQO8lsDP8wNgtpktJvglvTHRL8Ld95vZvwPzzGw7QTJryO8IEuPi8CZ9OQlMqerBnB43EjQTGsFosg0OMx4m7/1hIpZORKP/SodiZt3dvdLM+hL8Uj297hf8EZzrE0Cxu/9XUoNMsrif2YDfAGvc/eYI4vg6sNvdZ6b7syVaqpFIR/OEmfUCsoEfHWkSAXD3R8KE1NZ9wYIZGLMJakC/jSiOXQRzbkgnoxqJiIi0im62i4hIqyiRiIhIqyiRiIhIqyiRiIhIqyiRiIhIqyiRiIhIq/x/NPs3dAfGT+UAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"import numpy as np\n",
"\n",
"errors = 10**np.linspace(-1, -3, 25)\n",
"mean_dimensions = []\n",
"for eps in errors:\n",
" mean_dimensions.append(tn.mean_dimension(tn.round(t, eps=eps)))\n",
" \n",
"plt.figure()\n",
"plt.plot(np.log10(errors), mean_dimensions)\n",
"plt.xlabel('log10(rounding error)')\n",
"plt.ylabel('Mean dimension')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can also compute the *restricted* mean dimension, i.e. impose certain conditions on the set of tuples that intervene. For example, we can see which of two variables tends to show up more with higher-order terms:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"tensor(1.1355)\n",
"tensor(1.3664)\n"
]
}
],
"source": [
"print(tn.mean_dimension(t, mask=x))\n",
"print(tn.mean_dimension(t, mask=y))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The Dimension Distribution\n",
"\n",
"Last, the [*dimension distribution*](http://www3.stat.sinica.edu.tw/statistica/oldpdf/A13n11.pdf) gathers the relevance of $k$-tuples of variables for each $k = 1, \\dots, N$:"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"tensor([9.2223e-01, 7.2741e-02, 4.7444e-03, 2.6979e-04, 1.3797e-05, 6.4638e-07,\n",
" 2.8097e-08, 1.1428e-09, 4.3612e-11, 1.5722e-12, 5.3506e-14, 1.7216e-15,\n",
" 5.2287e-17, 1.4966e-18, 4.0144e-20, 1.0023e-21, 2.3034e-23, 4.7723e-25,\n",
" 8.6303e-27, 1.3940e-28])\n",
"Time: 0.11500120162963867\n"
]
}
],
"source": [
"start = time.time()\n",
"dimdist = tn.dimension_distribution(t)\n",
"print(dimdist)\n",
"print('Time:', time.time() - start)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It can be viewed as a probability mass function, namely the probability of choosing a $k$-variable tuple, if tuples are chosen according to their variance components. The expected value of this random variable is the mean dimension. Naturally, the dimension distribution must sum to $1$:"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"tensor(1.0000)"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sum(dimdist)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Again, we can extract the dimension distribution with respect to any mask:"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"tensor([0.0000, 0.6815, 0.2870, 0.0293, 0.0021, 0.0001, 0.0000, 0.0000, 0.0000,\n",
" 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,\n",
" 0.0000, 0.0000])"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tn.dimension_distribution(t, mask=y&z)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We just imposed that two variables ($y$ and $z$) appear. Note how, accordingly, the relevance of $1$-tuples has become zero since they have all been discarded."
]
}
],
"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.6.8"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
