# Statistical Inference of Minimally Complex Models

Clélia de Mulatier
Institute for Theoretical Physics and Informatics Institute
University of Amsterdam, Amsterdam, the Netherlands
Department of Physics and Astronomy
University of Pennsylvania, Philadelphia, USA
<EMAIL_ADDRESS>
&Paolo P. Mazza
Institute for Theoretical Physics
University of Tübingen
Auf der Morgenstelle 14, 72076 Tübingen, Germany
<EMAIL_ADDRESS>
&Matteo Marsili
The Abdus Salam International Centre for Theoretical Physics (ICTP)
Strada Costiera 11, I-34014 Trieste, Italy
<EMAIL_ADDRESS>

(September 2021)

###### Abstract

Finding the model that best describes a high dimensional dataset is a daunting
task. For binary data, we show that this becomes feasible when restricting the
search to a family of simple models, that we call Minimally Complex Models
(MCMs). These are spin models, with interactions of arbitrary order, that are
composed of independent components of minimal complexity (Beretta et al.,
2018). They tend to be simple in information theoretic terms, which means that
they are well-fitted to specific types of data, and are therefore easy to
falsify. We show that Bayesian model selection restricted to these models is
computationally feasible and has many other advantages. First, their evidence,
which trades off goodness-of-fit against model complexity, can be computed
easily without any parameter fitting. This allows selecting the best MCM among
all, even though the number of models is astronomically large. Furthermore,
MCMs can be inferred and sampled from without any computational effort.
Finally, model selection among MCMs is invariant with respect to changes in
the representation of the data. MCMs portray the structure of dependencies
among variables in a simple way, as illustrated in several examples, and thus
provide robust predictions on dependencies in the data. MCMs contain
interactions of any order between variables, and thus may reveal the presence
of interactions of order higher than pairwise.

_Keywords_ Statistical inference $\cdot$ Spin models $\cdot$ Information
theory

##### Significance statement.

Finding the model that best captures patterns hidden within noisy data is
fundamental for science but challenging. Here we propose a novel approach for
extracting the simplest explanation from binary data. To do so we introduce a
new family of models called Minimally Complex Models (MCMs). In these models,
variables are connected with interactions of arbitrary order that can identify
any patterns of data. Yet, MCMs are simple in information theoretic terms, and
provide robust predictions on dependencies between variables, that are easy to
falsify. Our approach contrasts with usual methods by focusing on simple
representation of data, opening new avenues of research in statistical
inference of high dimensional data.

##  1 Introduction

“All models are wrong, but some models are useful” [1]. This statement is
particularly appropriate in the context of statistical inference of high
dimensional data. Recent spectacular advances in machine learning have shown
that very complex models, such as deep neural networks [2], can be very
“useful” in learning hidden features in high dimensional data, making it
possible to generalize from examples. Models that encode the Laws of Nature,
refer to a different notion of “usefulness". As argued by Wigner [3], they
describe regularities – such as how bodies fall under the effect of gravity –
in simple forms involving only few variables and parameters. These
regularities occur in ways that are independent of many conditions that could
affect them. Simple models, such as Newton’s law, tell us more about
independence than about dependence. Their simplicity reflects specific
principles – such as invariances, symmetries and conservation laws – that are
easy to falsify. In this paper, we are interested in performing statistical
inference for binary data with models that display similar simplicity.

In information theory, the simplicity of a statistical model is defined
unambiguously in terms of the Minimum Description Length (MDL) complexity,
which is the minimal number of bits required to encode the model [4]. This
complexity is a measure of the number of different datasets that can be
described by the same parametric model [5], and therefore, of the flexibility
of the model in fitting broad types of data. A simple model thus fits well few
types of data and, as a consequence, is easy to falsify. In Bayesian
inference, the MDL complexity predicts how much a model would be penalised due
to its flexibility to data. In the context of binary data, Beretta et al. [6]
have studied the MDL complexity of spin models, which are maximum entropy
models with interactions of arbitrary order111Spin models include the so-
called Ising model as a special case, when only pairwise interactions are
considered.. They argue that, with the same number of parameters, the simplest
models are those for which statistical dependencies concentrates on the
smallest subset of variables. By contrast, pairwise spin models – which have
been used to model a variety of systems from neuronal activity [7, 8] to
voting outcomes [9] – turn out to be very complex. Interestingly, the simplest
spin models are also very easy to infer [6], while statistical inference of
pairwise spin models is known to be computationally challenging [10, 11, 12,
13, 14].

Figure 1: (colors online) Examples of Minimally Complex Models for a $6$-spin
system. Models are represented by diagrams: single spin variables are dots,
full in presence of a local field, empty otherwise; pairwise interactions are
blue lines; 3-spin interactions are orange triangles; and 4-spin interactions
are light blue polygons including four spins. Note that all three models are
displayed in a special basis.

In this paper, we show that this latter property of the simplest spin models
make Bayesian Model Selection (BMS) possible within a remarkably broad family
of simple models for binary data (for which maximum likelihood and model
evidence are easy to compute). These models, which we call Minimally Complex
Models (MCMs), are spin models formed of independent components of minimal
complexity. In some representations an MCM corresponds to a partition of the
variables in independent parts (or “communities”), with no interactions
between the parts and with all possible interactions within each part (see
Fig. 1). The best MCM for a given dataset thus captures a global structure of
dependencies and independencies in the data, providing robust predictions on
dependencies between variables. Besides being simple and easy to infer, MCMs
also enjoy properties that make statistical inference invariant under changes
in the representation of data. Finally, inferred MCMs can also be sampled from
with minimal computational effort.

In what follows, we formally define MCMs and discuss how to compare them in
the context of BMS. As the set of potential MCMs to compare is still huge (see
Fig. 2), we then propose different heuristics to find the best MCMs for data.
Our results are finally illustrated on several test cases, where one can
appreciate how our approach uncovers global data structures of dependencies
that are consistent with the discussed phenomena.

Figure 2: (colors online) Number of spin models as a function of the system
size $n$ for different families (see SI Sec. S2 for details on enumerations):
all spin models (green), all Minimally Complex Models (MCMs, violet), all
Independent Models (IM, orange), and all models with a Single Complete
Component (SCM, dark blue). MCM∗ and IM∗ indicate respectively the number of
MCM and IM that share the same preferred basis. For comparison, we also
reported the number of Pairwise Models (PM∗). The number of IMs and of MCMs
both grow exponentially with $n$, roughly as $2^{n^{2}}$, whereas the number
of PM∗ grows as $2^{n^{2}/2}$. The number of MCM∗ grows slower than $n^{n}$.
Note that the $y$-axis reports the logarithm base 10 of the number of models.
For example at $n=9$, there are of the order of $10^{153}$ models, but only
$10^{20}$ MCMs, which include $10^{18}$ IM and $10^{5}$ MCM∗; there are
$10^{13}$ PM∗.

## Bayesian Model Selection of Spin Models

A binary dataset
$\boldsymbol{\hat{s}}=(\boldsymbol{s}^{(1)},\ldots,\boldsymbol{s}^{(N)})$ is a
set of $N$ independent observations of spin configurations
$\boldsymbol{s}=(s_{1},\ldots,s_{n})$, where $s_{i}=\pm 1$. In the following,
we assume that each configuration $\boldsymbol{s}^{(i)}$ is independently
drawn from an unknown distribution, which we aim to infer. To do so, we
consider the following parametric family of probability distributions:

$p(\boldsymbol{s}\,|\,\boldsymbol{g},\mathcal{M})=\frac{1}{Z_{\mathcal{M}}(\boldsymbol{g})}\,\exp\left({\;\sum_{\boldsymbol{\mu}\in\mathcal{M}}g^{\boldsymbol{\mu}}\phi^{\boldsymbol{\mu}}(\boldsymbol{s})}\right)\,,$
(1)

where $\mathcal{M}=\\{\boldsymbol{\mu}_{1},\,\ldots,\boldsymbol{\mu}_{K}\\}$
is a set of $K$ interactions of arbitrary order, and the
$\boldsymbol{\mu}_{i}$’s are binary vectors encoding which spins are involved
in each interaction. Any choice of the set $\mathcal{M}$ of interactions
defines a spin model. The operators222The product is over all index $i$ such
that the $i$-th element of the binary vector $\boldsymbol{\mu}$ is $1$.
$\phi^{\boldsymbol{\mu}}(\boldsymbol{s})=\prod_{i\in\boldsymbol{\mu}}s_{i}\,$
are the product spin operators associated to the interactions
$\boldsymbol{\mu}$ of $\mathcal{M}$, while the conjugate parameters
$g^{\boldsymbol{\mu}}$ modulate the strength of the interactions. For example,
pairwise spin models contain interactions
$\phi^{\boldsymbol{\mu}}(\boldsymbol{s})=s_{i}s_{j}$ between pairs of spins.
Finally, the partition function $Z_{\mathcal{M}}(\boldsymbol{g})$ ensures
normalisation (see SI Sec. 1 for more details on spin models).

Equation (1) defines a complete family of models, capable of describing all
possible patterns of binary data with an appropriate choice of the set
$\mathcal{M}$ of operators [15]. In absence of prior knowledge on the system,
one must ideally compare the performance of all spin models to find which one
best describes a dataset $\boldsymbol{\hat{s}}$. In the Bayesian approach, the
best model $\mathcal{M}$ is the one achieving the largest posterior
probability $P(\mathcal{M}\,|\,\boldsymbol{\hat{s}})$, which is obtained with
Bayes’ theorem after computing the evidence [16]:

$P(\boldsymbol{\hat{s}}\,|\,\mathcal{M})=\int_{\mathbb{R}^{M}}\\!d\boldsymbol{g}\;\prod_{i=1}^{N}\,p\big{(}\boldsymbol{s}^{(i)}\,|\,\boldsymbol{g},\mathcal{M}\big{)}\,p_{0}(\boldsymbol{g}\,|\,\mathcal{M})\,,$
(2)

where $p_{0}(\boldsymbol{g}\,|\,\mathcal{M})$ is a prior distribution over the
parameters. We assume a uniform prior over the models $\mathcal{M}$, such that
the best model for data is also the one that maximizes the evidence.
Furthermore, we assume that $p_{0}(\boldsymbol{g}\,|\,\mathcal{M})$ takes the
form of Jeffreys’ prior [17]. With this choice, it has been shown [5] that in
the limit of large datasets ($N\to\infty$) the model $\mathcal{M}$ maximizing
the evidence in Eq. (2) is also the one providing the most succinct
description of the data in the Minimum Description Length (MDL) principle [18,
19, 4]. Indeed, in the large $N$ limit, one finds [5] that:

$\log P(\boldsymbol{\hat{s}}\,|\,\mathcal{M})\underset{N\to\infty}{\simeq}\log
P(\boldsymbol{\hat{s}}\,|\,\boldsymbol{\hat{g}},\mathcal{M})-\frac{K}{2}\log\left(\frac{N}{2\pi}\right)-c_{\mathcal{M}},$
(3)

where $\boldsymbol{\hat{g}}$ are the maximum likelihood parameters, and
$c_{\mathcal{M}}$ is the geometric complexity of model $\mathcal{M}$. The
second term is a first order model complexity which is proportional to the
number $K$ of interactions in the model. The sum of these two complexity terms
is the MDL complexity.

Yet, computing the evidence in Eq. (2) or the complexity $c_{\mathcal{M}}$ for
a generic spin model $\mathcal{M}$ is computationally challenging [6]. In
addition, the number of potential models to compare is huge – in an $n$-spin
system, there exists $2^{n}-1$ interactions of arbitrary order, and therefore
$2^{2^{n}-1}$ potential models –, making the exhaustive search among all spin
models unfeasible, even for small systems (see Fig. 2). A possible solution to
this hurdle is to restrict model selection to a smaller family of models, such
as pairwise models. For this family, model selection turns into a problem of
graph reconstruction [10, 12], that aims at identifying the network of
interactions among spins. A plethora of methods and results have been derived
within this framework [11, 14, 13]. Here, we propose to focus on a different
family of models, the minimally complex models (MCMs).

Figure 3: (colors online) Examples of MCMs and gauge transformations. Models
are represented with the same notations as in Fig. 1. All the models are MCMs.
In the first column, the first model is an independent model (IM) composed of
two independent operators, the second is a model composed of a single
independent complete component (ICC), and the last one is a MCM composed of
two ICCs. In each column, models of each row are obtained by the same gauge
transformation $\mathcal{T}$ of these models, and therefore have the same
properties than the original model (i.e., they are respectively, IM, MCM with
one ICC, and MCM with two ICCs).

## Minimally Complex Models (MCMs)

### Definition

MCMs are a subset of the spin models defined in Eq. (1), for which the set
$\mathcal{M}$ of interactions is the union of independent complete components
(ICCs) $\mathcal{M}_{a}$:

$\mathcal{M}=\bigcup_{a\in\mathcal{A}}\mathcal{M}_{a}\,.$ (4)

A complete component is a set $\mathcal{M}_{a}$ of interactions verifying
that, for any two $\boldsymbol{\mu},\boldsymbol{\nu}\in\mathcal{M}_{a}$, the
interaction $\boldsymbol{\mu}\oplus\boldsymbol{\nu}$ associated with the spin
operator
$\phi^{\boldsymbol{\mu}\oplus\boldsymbol{\nu}}(\boldsymbol{s})\doteq\phi^{\boldsymbol{\mu}}(\boldsymbol{s})\,\phi^{\boldsymbol{\nu}}(\boldsymbol{s})$
also belongs to $\mathcal{M}_{a}$. Furthermore, complete components are said
independent if they do not overlap. Therefore, in Eq. (4), for any two
distinct $(a,a^{\prime})\in\mathcal{A}$,
$\mathcal{M}_{a}\bigcap\mathcal{M}_{a^{\prime}}=\emptyset$. Figure 3 shows
several examples of MCMs.

### Preferred basis of an MCM

A basis $\boldsymbol{b}$ of a model $\mathcal{M}$ is a minimal set of $r$
independent operators, $\boldsymbol{b}=(\phi_{1},\cdots,\phi_{r})$, that can
generate all the operators of $\mathcal{M}$, i.e., any operators
$\phi\in\mathcal{M}$ can be written as a unique product of elements of
$\boldsymbol{b}$. The cardinal $r$ of such basis is called the rank of
$\mathcal{M}$ (necessarily $r\leq n$). We recall that a set $\boldsymbol{b}$
of operators is said independent if none of its operators can be obtained as a
product of other operators of $\boldsymbol{b}$ [6]. For example,
$\boldsymbol{b}=\\{s_{1},s_{2},s_{3}\\}$ and
$\boldsymbol{b}^{\prime}=\\{s_{1},s_{1}s_{2},s_{1}s_{2}s_{3}\\}$ are two sets
of independent operators, and both sets are possible bases for the model
$\mathcal{M}=\\{s_{1},s_{2},s_{1}s_{2},s_{3}\\}$.

By definition, a complete component $\mathcal{M}_{a}$ is a special case of
model that contains all the $2^{r_{a}}-1$ operators generated by its basis
$\boldsymbol{b}_{a}=(\phi_{1},\cdots,\phi_{r_{a}})$. As a direct consequence,
the number of operators of an MCM is

$\displaystyle K=\sum_{a\in\mathcal{A}}(2^{r_{a}}-1)\,.$ (5)

Note that the choice of the basis of an ICC is not unique. For example,
$\boldsymbol{b}$ and $\boldsymbol{b}^{\prime}$ both generate the same ICC, and
are therefore two possible bases for that ICC.

Finally, as a result of the independence between the complete component
$\mathcal{M}_{a}$ of an MCM $\mathcal{M}$, the union of chosen bases
$\boldsymbol{b}_{a}$ for each $\mathcal{M}_{a}$ forms a basis $\boldsymbol{b}$
of $\mathcal{M}$,

$\displaystyle\boldsymbol{b}=\cup_{a\in\mathcal{A}}\boldsymbol{b}_{a}\,.$ (6)

In this basis, $\mathcal{M}$ exactly corresponds to a partition of the basis
variables into ICCs (see examples in Fig. 1, and in the first column of Fig.
3). We refer to such basis as a preferred basis of an MCM. A direct
implication of Eq. (6) is that the rank of an MCM is
$r=\sum_{a\in\mathcal{A}}r_{a}$.

### Properties

Let us now discuss properties of MCMs. We first observe that, as a consequence
of Eq. (6), the state probability distribution of an MCM $\mathcal{M}$
factorizes over the probability distribution of its ICCs, when expressed in
terms of the operators
$\boldsymbol{b}(\boldsymbol{s})=\cup_{a}\,\boldsymbol{b}_{a}(\boldsymbol{s})$
of a preferred basis of $\mathcal{M}$:

$p(\boldsymbol{s}\,|\,\boldsymbol{g},\mathcal{M})=\frac{1}{2^{n-r}}\,\prod_{a\in\mathcal{A}}p_{a}\left(\boldsymbol{b}_{a}(\boldsymbol{s})\,|\,\boldsymbol{g}_{a},\mathcal{M}_{a}\right).$
(7)

The prefactor corresponds to a $1/2$ probability for each of the $(n-r)$ spin
variables not modeled by $\mathcal{M}$ (see for instance the spins represented
by empty dots in the MCMs of the first column of Fig. 3). Note that each
$p_{a}$ is a probability distribution over $r_{a}$ spin variables.
Importantly, Eq. (7) implies that the evidence of an MCM factorizes over its
ICCs. In addition, the evidence of each ICC is remarkably simple to compute
(see Methods 1.1), leading to the following expression for the evidence of an
MCM $\mathcal{M}$ for any dataset $\boldsymbol{\hat{s}}$:

$P(\boldsymbol{\hat{s}}\,|\,\mathcal{M})=\frac{1}{2^{N(n-r)}}\prod_{a\in\mathcal{A}}\,\frac{\Gamma(2^{r_{a}-1})}{\Gamma(N+2^{r_{a}-1})}\;\prod_{\boldsymbol{b}_{a}}\frac{\Gamma(k_{\boldsymbol{b}_{a}}+\frac{1}{2})}{\sqrt{\pi}}\,,$
(8)

where $\Gamma$ is the gamma function, $N$ is the number of datapoints in
$\boldsymbol{\hat{s}}$, and $k_{\boldsymbol{b}_{a}}$ is the number of times
that the basis operators take the value ${\boldsymbol{b}_{a}}$ over the
dataset. The maximum likelihood distribution also takes a simple form (see
Methods 1.1):

$P(\boldsymbol{\hat{s}}\,|\,\boldsymbol{\hat{g}},\,\mathcal{M})=\frac{1}{2^{N(n-r)}}\prod_{a\in\mathcal{A}}\,\prod_{\boldsymbol{b}_{a}}\,\left(\frac{k_{\boldsymbol{b}_{a}}}{N}\right)^{k_{\boldsymbol{b}_{a}}}\,.$
(9)

This makes sampling from the maximum likelihood distribution a very easy task
(see Methods 1.3).

Another important property is that the family of MCMs is invariant under the
gauge transformations (GTs) defined in Ref. [6]. A GT $\mathcal{T}$ is a
bijection of the set of states into itself (automorphisms), where a new state
$\boldsymbol{s}^{\prime}$ is defined as a function of the original state
$\boldsymbol{s}$ by a set of $n$ independent operators:
$\boldsymbol{s}^{\prime}=\mathcal{T}(\boldsymbol{s})=\left(\phi_{1}(\boldsymbol{s}),\ldots,\phi_{n}(\boldsymbol{s})\right)$.
Such transformations are also bijections of the set of all operators into
itself333A transformed operator $\phi^{\prime}$ is obtained from the original
operator $\phi$ by
$\phi^{\prime}(\boldsymbol{s}^{\prime})=\phi(\mathcal{T}^{-1}(\boldsymbol{s}^{\prime}))$.,
and therefore of the set of models into itself444A transformed model
$\mathcal{M}^{\prime}=\mathcal{T}(\mathcal{M})$ is obtained by transforming
each operator of the original model $\mathcal{M}$.. Notice that the order of
an operator, i.e. the number of spins that occur in it, is not invariant under
GTs (see Fig. 3 for examples of GTs). However, the mutual relation between the
operators of a model is preserved under a GT555This mutual relation between
operators of a model is called the loop structure of the model in Ref. [6].
Models with he same loop structure have the same Fisher information matrix.
[6], i.e., if $\phi_{1}$ and $\phi_{2}$ are two operators and
$\phi_{1+2}=\phi_{1}\phi_{2}$, then the gauge transformed operator
$\phi_{1+2}^{\prime}$ is still the product of the transformed operators
$\phi_{1}^{\prime}$ and $\phi_{2}^{\prime}$. In particular, under a GT an ICC
maps into an ICC with same rank, and an IM maps into an IM with same rank. And
as a consequence, applying a GT to an MCM returns an MCM with the same number
ICCs and same rank sequence (see Fig. 3). On the contrary, the family of
pairwise models is not invariant under GTs. In BMS, invariance under GTs of
the considered family of models ensures that the outcome of the inference
process is independent of the gauge in which the data is expressed.

Finally, Ref. [6] shows that, among all models with the same number of
parameters, ICCs are the models with the lowest MDL complexity. As MCMs are
formed by the union of independent ICCs, their MDL complexity is the sum of
the complexity of each ICC, and is thus expected to be small. We conjecture
that MCMs are the models with the lowest MDL complexity among all spin models
with the same number of parameters and the same rank. We verified that this is
true for all models with up to $4$ spin variables (see Fig. S1). Although this
statement has yet to be proven in the general case, our belief in its likely
correctness motivated our choice of the name minimally complex models.

### Bayesian Model Selection

All together, these properties are good reasons for restricting BMS to the set
of MCMs. First, Eq. (7) shows that finding the best MCM provides sharp
predictions on independencies between variables of different ICCs. Second, the
fact that the computation of the evidence in Eq. (8) is straightforward from
the data greatly simplifies the task of BMS, and comparison between models can
be done without inferring any parameter and directly on the basis of their
evidence. Finally, invariance under GTs ensures that the selected models are
independent of the ways in which data is represented.

Yet the number of MCMs is still huge, even for moderate values of $n$ (growing
roughly as $2^{n^{2}}$, see Fig. 2). To find the MCM with maximal evidence, we
propose to divide the search in two steps. First, we look for the best model
among all Independent Models (IM), which are MCMs with $r_{a}=1$ for all
$a\in\mathcal{A}$ (see examples in Fig. 3). The maximization of the evidence
over this family can be done by finding the set $\boldsymbol{b}^{*}$ of the
$n$ most biased independent operators in the dataset (see Methods 1.2).
Although the number of IMs is also huge ($\sim 2^{n^{2}}$), finding the best
IM with this procedure only takes of the order of $n\,2^{n}$ operations. Once
the best IM $\boldsymbol{b}^{*}$ is singled out, we reduce our search for the
best MCM to the subset of models that admit $\boldsymbol{b}^{*}$ as a
preferred basis. This entails partitioning the vector $\boldsymbol{b}^{*}$
into ICCs, and finding which such partition maximizes the evidence. The
outcome of this algorithm can be represented as a factor graph with three
layers (see Fig. 4b), one connecting the original variables $\boldsymbol{s}$
with the basis vector $\boldsymbol{b}^{*}$, and the other connecting the
latter to the ICCs $\mathcal{M}_{a}$. A program implementing this algorithm
has been made available by the author at Ref. [20]. Note that the number of
IMs is not much smaller than the number of MCMs, whereas the number of MCMs
with the same basis (denoted as MCM* in Fig. 2) is much smaller. Indeed, the
number of such models is given by the Bell number of $n$, which grows slower
than $(n/\log(n))^{n}$ [21]. This makes exhaustive search possible for
moderate system sizes ($n\lesssim 15$). Larger systems require further
heuristics that will be discussed below.

Figure 4: Analysis of the US Supreme Court Data. Justices are represented by
circles labelled by their initials: Ruth Bader Ginsburg (RG), John P. Stevens
(JS), David Souter (DS), Stephen Breyer (SB), Sandra Day O’Connor (SO),
William Rehnquist (WR), Anthony Kennedy (AK), Clarence Thomas (CT), Antonin
Scalia (AS). Colours represent their political orientation taken from Ref.
[9]. a) The best MCM, represented in the original basis of variables (i.e.,
the justices’ votes). The best IM is composed of $8$ pairwise interactions,
represented by links between the nodes with width proportional to their
strength, and $1$ single body interaction represented by a square on CT. The
strongest interaction has $\langle s_{\rm CT}s_{\rm WR}\rangle\simeq 0.86$,
whereas the weakest has $\langle s_{\rm CT}\rangle\simeq-0.45$. The three
shaded circles indicate the overall best MCM, which was computed in the basis
of these interactions. The red and blue dotted squares indicate the best MCM
based on the original basis of spin variables (i.e., for which
$\boldsymbol{b}=\boldsymbol{s}$). b) Factor graph representation of the best
MCM. Spin variables $\boldsymbol{s}$ are represented by circles. The inference
procedure first identifies the best basis $\boldsymbol{b}^{*}$ of independent
operators, denoted by squares, and then the best clustering of these operators
into ICCs $\mathcal{M}_{a}$, denoted by triangles. Each square in the second
layer corresponds to one of the basis interactions represented in Fig. a.

## MCM selection for small systems

### Voting data from the US Supreme Court [9]

Let us first illustrate our method with the analysis of voting data from the
US Supreme Court database [22]. The dataset we consider is composed of the
votes of the $n=9$ justices on $N=895$ cases debated during the second
Rehnquist Court. As in Ref. [9], voting outcomes were re-labeled using the
political inclination (if known) of each debated case, such that, for each
case, each judge casts a vote $s_{i}=+1$ if the decision is conservative-
oriented, or $s_{i}=-1$ if it is liberal-oriented. Ref. [9] showed that a
fully connected pairwise model is able to explain higher order features of the
data. Ref. [23] further analysed the data within a more general scheme,
showing that pairwise interactions are indeed prevalent in this dataset. The
size of the system ($n=9$) is sufficiently small to perform the exhaustive
search described in the previous section666For $n=9$, the program in [20]
takes less than a second to find the best MCM on a simple laptop..

The best IM is displayed in Fig. 4.a and has a very large log-evidence ($\log
E=-3327$) compared to the usual IM composed of $9$ single spin operators
(i.e., the IM based in the original basis of the data
$\boldsymbol{b}=\boldsymbol{s}$, for which $\log E=-5258$). Interestingly, we
found that the $9$ most biased independent operators of the data are all
single- or two-body, although the selection procedure takes into account
interactions of all orders. This confirms the prevalence of low order
interactions in the system. Moreover, these $9$ interactions alone account for
more than $80\%$ of the information captured by the fully connected pairwise
model with $45$ interactions, and therefore they provide a meaningful
representation of the data. With the exception of the single spin interaction
on CT, all these interactions are pairwise. They connect justices with similar
political orientation [9] and define a network that spans the entire political
spectrum. Interestingly, the strength of the interactions increases towards
the extremes of the political spectrum and is weaker at the center. Without
the single spin interaction, which is much weaker than the other interactions,
the systems would be invariant under spin reversal ($s_{i}\longrightarrow-
s_{i}$). The single spin interaction on CT introduces a bias towards
conservative votes.

The best MCM can be found from the best IM by using algorithms for generating
all possible partitions of a set [24, 25, 26]. There are $115975$ different
MCMs that can be generated from a given IM with $n=9$. For the US Supreme
court data, the best MCM is composed of three ICCs (see Fig. 4). The first
consists solely of the interaction between AS and CT, suggesting that AS’s
voting patterns conditioned on that of CT is mostly independent of all the
others. The second ICC groups the three interactions between judges on the
liberal side of the political spectrum, suggesting that, conditioned on the
vote of SB, the voting patterns of RG, JS and DS is mostly independent of the
votes of the other judges. Finally, the third ICC groups interactions on the
conservative-center side of the political spectrum. For comparison, Fig. 4.a
also shows the best MCM whose preferred basis is the original basis of the
data (i.e., with $\boldsymbol{b}=\boldsymbol{s}$). This MCM divides the court
into two independent components that happen to match with the justices’
political orientations. Yet the evidence of this model ($\log
P(\hat{s}|\mathcal{M})=-3156.71$) is considerably smaller than that of the MCM
build from the best IM $\boldsymbol{b}^{*}$ ($\log
P(\hat{s}|\mathcal{M})=-3300.97$). The evidence of the best MCM in the
original basis is actually smaller than the evidence of most MCM with same
rank in the $\boldsymbol{b}^{*}$ basis, which highlights the advantage of
choosing the best IM as a basis for searching for the best MCM (see Fig. S2
for mode details). Finally, note that the best MCM identified by our algorithm
has $K_{\rm MCM}=(2^{5}-1)+(3^{3}-1)+1=39$ operators, which is smaller than
the number of parameters $K_{\rm pair}=9\times 10/2=45$ of a fully connected
pairwise model.

### Boolean description of a binary dataset

Next, we tested our approach on an artificial dataset of $n=9$ variables and
$N=10^{5}$ datapoints that we generated by sampling randomly three different
states, $S_{1}$, $S_{2}$ and $S_{3}$, with different probabilities,
respectively $p_{1}=0.2$, $p_{2}=0.3$ and $p_{3}=0.5$ (see Fig. 5.a). There
were no additional patterns or sources of noise in the data. The size of the
system ($n=9$) was chosen sufficiently small to easily perform the exhaustive
search described in the previous section. The factor graph of the best MCM
found for this dataset is shown in Fig. 5.b. It is composed of $8$ ICCs, $7$
of which are based on a single operator (green triangles).

Once it is fitted, this MCM can be used as a generative model. As discussed in
Methods 1.1, fitting an MCM does not require fitting the vector of parameters
$\boldsymbol{g}$, because one can directly fit the probability of the states
described by each ICC. Indeed, for an ICC $\mathcal{M}_{a}$, the likelihood is
maximized when the probability of occurrence of the states of
$\boldsymbol{b}_{a}$ in the model is equal to their probability of occurrence
in the data (see Eq. (9) and Methods 1.1). We thus obtain that all the single-
operator ICCs (green triangles in Fig. 5.b) are associated to statement that
are true with probability $1$ over the entire dataset, i.e. they are
associated to independent tautologies of the data. For instance, in this
dataset, the operator $\phi_{1}=s_{8}s_{9}$ is equal to $1$ with probability
$1$ and the operator $\phi_{5}=s_{3}s_{9}$ is equal to $-1$ with probability
$1$, i.e., the two statements $\phi_{1}=1$ and $\phi_{5}=-1$ are true over the
entire dataset. After taking into account the $7$ independent tautologies, we
are left with only four states that could be generated by the best MCM (see
Fig. 5.c), each of them corresponding to a different value of the basis
elements $\boldsymbol{b}_{a}=(\phi_{8},\phi_{9})$ of the last ICC (orange
triangle). Fitting the probabilities of occurrence of $\boldsymbol{b}_{a}$ on
the data, we obtained that $P(\phi_{8}=-1,\phi_{9}=1)=0.2$, $P(1,-1)=0.3$,
$P(1,1)=0.5$ and, by normalization, $P(-1,-1)=0$ (there are only three
independent probabilities). This finally provides a full description of the
rules used to generate the data. In this case, we observe that the best MCM
provides the shortest complete description of the data in terms of spin
models. Note that because each spin operator is associated to an elementary
boolean functions777A spin state can be re-written as $s_{i}=(-1)^{x_{i}}$
where $x_{i}\in\\{-1,+1\\}$. Thus the spin operator
$\phi^{\boldsymbol{\mu}}(\boldsymbol{s})=\prod_{i\in\boldsymbol{\mu}}s_{i}$
corresponds to the boolean function
$\phi^{\boldsymbol{\mu}}(\boldsymbol{x})=\bigoplus_{i\in\boldsymbol{\mu}}x_{i}$,
where $\oplus$ is the XOR operator, and
$\phi^{\boldsymbol{\mu}}(\boldsymbol{s})=(-1)^{\phi^{\boldsymbol{\mu}}(\boldsymbol{x})}$.,
the best MCM therefore also provides a complete decomposition of the data in
terms of such boolean functions. each associated with a parameter that is the
probability that this function takes value $1$ over the dataset ($7$ of which
being just probability $1$ or $0$).

Figure 5: Analysis of a simple artificial dataset. a) We generated a dataset
composed of three different states ($S_{1}$, $S_{2}$ and $S_{3}$) each sampled
with a different probability (respectively $p_{1}$, $p_{2}$ and $p_{3}$). Dark
pixels correspond to the spin state $s_{i}=-1$, and light pixels to
$s_{i}=+1$. b) Factor graph of the best MCM found for this dataset. c) Once it
is fitted, the best MCM can generate four different states. Three of them are
the original states and are sampled with the same probability than their
occurrence in the original dataset, the last one is sampled with probability
$0$.

## MC selection, heuristics for large systems

When $n$ is large and exhaustive search is unfeasible, we propose the
following heuristics. In order to find the best IM $\boldsymbol{b}^{\star}$,
we start from an initial basis $\boldsymbol{b}$ (e.g.,
$\boldsymbol{b}=\boldsymbol{s}$), and we build all interactions up to a chosen
order $k$. Among these interactions, we identify the set
$\boldsymbol{b}^{\prime}$ of $n$ independent operators that is maximally
biased, and replace $\boldsymbol{b}$ by $\boldsymbol{b}^{\prime}$. We repeat
this procedure until we find $\boldsymbol{b}^{\prime}=\boldsymbol{b}$. This
final basis is our best choice for the IM $\boldsymbol{b}^{\star}$. Although
the exploration of the space of IM is limited by the choice of $k$, the
iteration of this procedure is able, in principle, to explore the space of
operators to any order.

Next, we apply an hierarchical merging procedure to find the optimal MCM. We
start from the IM based on the basis operators $\boldsymbol{b}^{*}$ identified
above, which is an MCM with $n$ ICC of rank $r_{a}=1$. We merge two ICCs
$\mathcal{M}_{a}$ and $\mathcal{M}_{a^{\prime}}$ in all possible ways. Among
these combinations, we identify the pair that yields a maximal increase of the
evidence in Eq. (8) and merge the corresponding ICCs. This procedure generates
an approximation of the MCM that achieves a maximal value of the evidence
along the hierarchical merging process, as the number of ICCs varies from $n$
to $1$.

### Big Five Personality Test [27]

We applied this algorithm to several datasets, performing the iterative search
for the best basis at the $k=4^{\rm th}$ order. Figure 6 reports the resulting
MCM for the Big Five Personality Test [27]. The test consists of $n=50$
questions designed to probe the personality of individuals along five
different dimensions, that have been suggested as the main traits describing
personality. These traits are extraversion, neuroticism, agreeableness,
conscientiousness and openness to experience. Each trait is evaluated from the
answers to ten questions, on a scale of one (disagree) to five (agree), that
can be either positively or negatively associated with the trait888For
example, agreeableness is probed by questions such as “I sympathise with
others’ feelings” or “I insult people”.. Data on $N=1\,013\,558$ samples were
taken from Ref. [28], to which we refer for more details. We converted each
answer in binary format, depending on whether it was positively or negatively
associated with the trait with respect to the average score across the whole
sample. For this dataset, inference within the family of MCMs can reveal
whether or not the data confirms the hypothesis that these questions probe the
respondent’s personality along five dimensions, and how these dimensions are
associated with the questions of the test. Our result, obtained with the
heuristic algorithm, is shown in Fig. 6 top. Here as well, the best basis
$\boldsymbol{b}^{\star}$ contains exclusively single-body and two-body
operators. With the exception of one999We found a pairwise interaction between
the question “I’m full of ideas”, that should probe openness to experience,
and “I have little to say”, that probes (negatively) extraversion., all these
interactions are confined to questions relative to the same trait. The best
MCM we found has five ICCs, each mainly grouping together interactions related
to the same trait, apart from one ICC that significantly mixes
conscientiousness and openness. The log-evidence of this MCM ($\log
E\simeq-2.596\cdot 10^{7}$) is however still smaller than the log-evidence of
the MCM with 5 ICCs exactly grouping together all the interactions related to
the same trait ($\log E\simeq-2.518\cdot 10^{7}$). This result shows that this
latter MCM is indeed one of the best (see distribution of the log-evidence of
MCMs with 5 ICCs in Fig. S3), and overall, this analysis confirms that the
questions of the Big Five Personality Test accurately probe the desired five
independent traits [27]. Our analysis also found dependencies among the
questions within each trait, which could be used to further improve the design
or the analysis of the test. Finally, this analysis suggests that a measure
based on the log-evidence of MCMs can be used to assess the quality of a
personality test.

Figure 6:  Factor graph representation of the best MCM found with our
heuristic algorithm for the Big Five Personality Test dataset.

### Modeling the MNIST database [29]

Finally, we performed BMS with MCMs on the MNIST database [29], which is
composed of $N=60000$ images of hand-written digits. In order to reduce the
data to a manageable size, given our computational resources101010All
calculations were performed on a laptop computer., we coarse grained it in
cells of $2\times 2$ pixels, that were converted into binary values by
applying a threshold111111If the sum of the grey levels of the four pixels
exceeds 400, the coarse grained pixel is assigned the value one, otherwise it
is zero. Our code employs bitwise operations on 16 bytes integers, thus
limiting the size of the systems we could handle to $n\leq 128$.. We focus on
a central zone of $n=11\times 11=121$ pixels. Fig. 7.a shows the evidence as a
function of the number $A=|\mathcal{A}|$ of ICCs at different stages of the
hierarchical merging procedure. A maximum is achieved at $A=15$. The structure
of the inferred MCM at $A=15$ is rather complex, with several high order
interactions, which are grouped in relatively large ICCs (see Fig. 7.b and c).

The aim of this exercise is to test the efficiency of the inferred MCMs in
generalization. For this purpose, we sample digits from the maximum likelihood
distribution at different values of $A$. In order to generate a new digit
$\boldsymbol{s}$, we perform the following steps: i) we sample a random value
of $\boldsymbol{b}_{a}$ directly from the data, for each $a\in\mathcal{A}$;
ii) the sampled state
$\boldsymbol{b}^{*}=\\{\boldsymbol{b}_{a},\,a\in\mathcal{A}\\}$ is then
obtained by combining the states of each $\boldsymbol{b}_{a}$; finally iii) we
get the sampled digit $\boldsymbol{s}$ in the original basis, by inverting the
GT $\mathcal{T}$ that maps $\boldsymbol{s}$ to $\boldsymbol{b}^{*}$ and using
$\boldsymbol{s}=\mathcal{T}^{-1}(\boldsymbol{b}^{*})$. As shown in Methods
1.3, this last steps requires the inversion of an $n\times n$ binary matrix
(modulo 2), which needs only to be performed once. Step i) exploits the fact
the ICC $\mathcal{M}_{a}$ is complete on $\boldsymbol{b}_{a}$, and therefore,
computing $\boldsymbol{b}_{a}$ on a randomly drawn digit from the dataset
$\boldsymbol{\hat{s}}$ generates values of $\boldsymbol{b}_{a}$ that are
distributed according to Eq. (9). Hence, a new state of $\boldsymbol{b}^{*}$
can be generated in a straightforward manner from $A$ independent draws from
the dataset. And therefore, sampling from the likelihood of the inferred MCM
is remarkably simple. Notice that, for $A=1$ this procedure amounts to
sampling digits directly from the original dataset. For $A>1$, sampling
generates new patterns, as shown in the inset of Fig. 7 a). Although the
sampled images have some structure, their resemblance to digits fades away as
$A$ increases.

Figure 7:  a) Evidence of the selected MCM as a function of the number
$A=|\mathcal{A}|$ of ICCs during the merging procedure. The maximum is
attained at $A=15$. The inset shows sample digits drawn from the maximum
likelihood distribution at different values of $A$. The top row refers to the
original data ($A=1$). b) Histogram of the order of the interactions in the
best IM. There is a significant number of interactions that are of order
higher than pairwise. c) Histogram of the rank $r_{a}$ (“size”) of the $15$
ICCs of the best MCM found with the merging procedure.

## Discussion

Finding the model that best describes a dataset within Bayesian model
selection is a daunting task. We show that this task is considerably
simplified when the selection is restricted to the family of MCMs. These
models are simple in terms of their description length complexity, and they
are easy to infer. BMS within MCMs probes interaction of arbitrary order and
can reveal the presence of high order interactions or confirm that low order
interactions are prevalent in a dataset. In addition, MCMs disentangle the
structure of statistical dependence into independent components. As such, it
can be a useful preprocessing step to divide the inference problem into
smaller problems that can be analysed in more detail. Furthermore, this
structure provides a notion of locality that could be used in coarse-graining
procedures in the context of data renormalisation [30], going beyond pairwise
correlations. Finally, the invariance under GTs of the family of MCMs ensures
that the inference process is independent of how data is
represented121212Consider, for a purely illustrative purpose, an idealised
problem of inference in a gene regulatory network. Assume that whether gene
$i$ is expressed ($s_{i}=+1$) or not ($s_{i}=-1$) depends on whether its
$t_{i}$ transcription regulators are bound ($\sigma_{i,a}=+1$) or not
($\sigma_{i,a}=-1$) to the regulator binding region, i.e. that
$s_{i}=f_{i}(\sigma_{i,1},\ldots,\sigma_{i,t_{i}})$ is a boolean function of
the $\sigma_{i,a}$’s. If the relation between $\boldsymbol{s}$ and
$\boldmath{\sigma}$ is a GT, then our approach ensures that inference of the
gene regulatory network based on a dataset $\boldsymbol{\hat{s}}$ of gene
expression should give the same results as inference based on a dataset
$\boldmath{\hat{\sigma}}$ of binding on regulatory regions.. This is not true
when model selection is restricted to a family of models that is not invariant
under GTs, such as the family of pairwise models. In this case, it may be hard
to say whether the structure of the inferred model reflects statistical
dependencies of the data or the constraints imposed by the considered family
of models.

Our approach departs from the literature on graphical model reconstruction
[10, 31, 12, 13] in important ways. The latter aims at retrieving a model from
a dataset generated from it, by identifying which interactions among a pre-
assigned set (e.g. pairwise) are present. Here we do not make any assumptions
on the interactions present in the models. Most importantly, we are not
interested in reconstructing a specific model from which the data could have
been generated. Instead, we aim to find the most likely simple model for data,
which can inform us on the real structure of dependencies in the data. A
critical aspect in graphical model reconstruction is parameter estimation, for
which different approximate algorithms have been proposed [14]. We show that
there is no need to infer model parameters in order to perform model selection
among MCMs, which is a huge computational advantage. Moreover, MCMs can be
directly compared on the basis of their evidence, which allows to balance
goodness-of-fit and model complexity. In contrast, the model complexity of
graphical models is difficult to compute, which calls for adopting
regularisation methods [31, 14], which have no clear information theoretical
foundation. Finally, the set of models that can be compared in BMS within the
family of MCMs is much larger than the number of possible pairwise models
usually considered (see Fig. 2).

In summary, this paper offers a novel perspective on statistical inference of
high dimensional data, by focusing on simple representations of data. Further
extensions beyond binary variables, as well as the development of efficient
algorithms to find the MCM with maximal evidence for large systems, are
promising avenues of research.

## Materials and Methods

###  1.1 Maximum Likelihood and Evidence of an ICC

Let us consider an ICC $\mathcal{M}_{a}$ with rank $r_{a}$, and a basis
$\boldsymbol{b}_{a}(\boldsymbol{s})=(\phi_{1}(\boldsymbol{s}),\ldots,\phi_{r_{a}}(\boldsymbol{s}))$
of $\mathcal{M}_{a}$. The probability distribution of the states of the system
under $\mathcal{M}_{a}$ can be re-written in the basis of operators
$\boldsymbol{b}_{a}$ as:

$P(\boldsymbol{s}\,|\,\boldsymbol{g},\mathcal{M}_{a})=\frac{1}{2^{(n-r_{a})}}\,q\left(\boldsymbol{b}_{a}(\boldsymbol{s})\right),$
(10)

where
$q(\boldsymbol{b}_{a})=P(\boldsymbol{b}_{a}\,|\,\boldsymbol{g},\,\mathcal{M}_{a})$
is a normalised probability distribution over the states of
$\boldsymbol{b}_{a}\in\\{\pm 1\\}^{r_{a}}$. For an ICC $\mathcal{M}_{a}$,
there is a bijective map between the $2^{r_{a}}-1$ parameters
$g^{\boldsymbol{\mu}}$ of $\mathcal{M}_{a}$ and the $2^{r_{a}}-1$ states
probabilities $q(\boldsymbol{b}_{a})$ (the last probability being fixed by the
normalization). More precisely, the value of the parameters
$g^{\boldsymbol{\mu}}$ can be obtained from the states probabilities
$q(\boldsymbol{b}_{a})$ by taking the logarithm of Eq. (10), multiplying by
$\phi^{\boldsymbol{\mu}}(\boldsymbol{s})$ and summing over all the states of
$\boldsymbol{s}$. Finally using the orthogonality relation between operators,
one finds that for all interaction $\mu\in\mathcal{M}_{a}$:

$\displaystyle g^{\mu}$
$\displaystyle=\frac{1}{2^{n}}\,\sum_{\boldsymbol{s}}\phi^{\mu}(\boldsymbol{s})\,\log
q\left(\boldsymbol{b}_{a}(\boldsymbol{s})\right)\,,\qquad\forall\mu\in\mathcal{M}_{a}.$
(11)

Note that for $\mu\not\in\mathcal{M}_{a}$ this equation yields $g^{\mu}=0$,
because the function $q\left(\boldsymbol{b}_{a}(\boldsymbol{s})\right)$ is
expressed solely in terms of the operators $\mu\in\mathcal{M}_{a}$.
Statistical inference of ICCs, and therefore of MCMs, is more easily done in
the $q$-representation, i.e. in terms of the parameters
$q(\boldsymbol{b}_{a})$. The result can then be re-written in the original
spin model representation using Eq. (11).

Using Eq. (10), the likelihood function of an ICC $\mathcal{M}_{a}$ for a
dataset $\boldsymbol{\hat{s}}$ can be written as:

$P(\boldsymbol{\hat{s}}\,|\,\boldsymbol{g},\mathcal{M}_{a})=\frac{1}{2^{N(n-r_{a})}}\,\prod_{\boldsymbol{b}_{a}}\,q(\boldsymbol{b}_{a})^{k_{\boldsymbol{b}_{a}}}\,,$
(12)

where the product is over all possible $2^{r_{a}}$ values of
$\boldsymbol{b}_{a}$, and where $k_{\boldsymbol{b}_{a}}$ is the number of
times the basis operators take a given value ${\boldsymbol{b}_{a}}$ over the
entire dataset. Thanks to the bijection between the $2^{r_{a}}-1$ parameters
$\boldsymbol{g}$ and the $2^{r_{a}}-1$ parameters $q(\boldsymbol{b}_{a})$, the
maximum likelihood can be obtained by maximizing Eq. (12) directly over the
parameters $q(\boldsymbol{b}_{a})$. From this maximization, one finds that the
maximum likelihood estimate of $q$ is:

$\hat{q}({\boldsymbol{b}_{a}})=\frac{k_{\boldsymbol{b}_{a}}}{N}.$ (13)

Combined with Eq. (7) and Eq. (12), this leads to the result in Eq. (9) for
the maximum likelihood of an MCM.

The calculation of the evidence of an ICC is likewise straightforward. We
first exploit the invariance of Jeffreys prior under reparametrization [17],
which takes the following form in the $q$-representation:

$P_{0}(\boldsymbol{q}\,|\,\mathcal{M}_{a})=\frac{\Gamma(2^{r_{a}-1})}{\pi^{2^{r_{a}-1}}}\prod_{\boldsymbol{b}_{a}}\sqrt{q(\boldsymbol{b}_{a})}\;\delta\Big{(}\sum_{\boldsymbol{b}_{a}}q(\boldsymbol{b}_{a})-1\Big{)},$
(14)

where we denote by $\boldsymbol{q}$ the vector of $2^{r_{a}}$ probabilities
$q(\boldsymbol{b}_{a})$. Combined with the expression of the likelihood in Eq.
(12), this yields:

$\displaystyle P(\boldsymbol{\hat{s}}|\mathcal{M}_{a})$
$\displaystyle=\int\\!d\boldsymbol{g}\;P(\boldsymbol{\hat{s}}\,|\,\boldsymbol{g},\,\mathcal{M}_{a})\,P_{0}(\boldsymbol{g}\,|\,\mathcal{M}_{a})$
(15)
$\displaystyle=\frac{1}{2^{N(n-r_{a})}}\frac{\Gamma(2^{r_{a}-1})}{\pi^{2^{r_{a}-1}}}\int_{0}^{1}d\boldsymbol{q}\prod_{\boldsymbol{b}_{a}}q(\boldsymbol{b}_{a})^{k_{\boldsymbol{b}_{a}}+\frac{1}{2}}\delta\big{(}\sum_{\boldsymbol{b}_{a}}q(\boldsymbol{b}_{a})-1\big{)}$
(16)
$\displaystyle=\frac{1}{2^{N(n-r_{a})}}\frac{\Gamma(2^{r_{a}-1})}{\Gamma(N+2^{r_{a}-1})}\prod_{\boldsymbol{b}_{a}}\frac{\Gamma(k_{\boldsymbol{b}_{a}}+\frac{1}{2})}{\sqrt{\pi}}.$
(17)

This calculation easily generalizes to MCMs, in which the components
$\mathcal{M}_{a}$ are associated to independent basis elements
$\boldsymbol{b}_{a}$. We recall that the likelihood, the prior and the
evidence of an MCM factorizes over its independent components, using Eq. (7).
This finally leads to the expression of the evidence of an MCM given Eq. (8).

###  1.2 Search for the best IM

An independent model (IM) is an MCM in which all the components
$\mathcal{M}_{a}$ only contain a single operator ($r_{a}=1$ for all
$a\in\mathcal{A}$). Let us focus on the case where the number of components
$|\mathcal{A}|$ is equal to the number of spins $n$. Model selection among
these models only requires to compare their likelihood, because they all have
the same complexity terms. The maximum log-likelihood of an independent model
takes the simple form:

$\log
P(\boldsymbol{\hat{s}}\,|\,\boldsymbol{\hat{g}},\mathcal{M})=-N\sum_{a=1}^{n}H[m_{a}]$
(18)

where the sum is over the $n$ components, and where

$H[m]=-\frac{1+m}{2}\log\frac{1+m}{2}-\frac{1-m}{2}\log\frac{1-m}{2}$ (19)

and

$m_{a}=\frac{1}{N}\sum_{i=1}^{N}b_{a}(\boldsymbol{s}^{(i)})$ (20)

is the bias of operator $b_{a}$. The function $H[m]$ is maximal for $m=0$ and
it achieves its minimal value $H[m]=0$ when $m=1$ or $m=-1$. This implies that
the most likely independent model is given by the most biased set
$\boldsymbol{b}$ of independent operators.

The algorithm to find these operators is a recursive one: Let $b_{1}$ be the
most biased operator $\phi^{\mu}$, i.e. the one with the largest average value
of $|\phi^{\mu}|$ in the dataset. Add the next most biased operator $b_{2}$.
The third operator $b_{3}$ is the most biased one, excluding the operator
$b_{1}b_{2}$, which is not independent from the previous two. Proceeding in
this way, at step $r$, let $\boldsymbol{b}_{<r}=(b_{1},\ldots,b_{r-1})$ be the
set of independent operators already identified. This divides the set of all
operators in the subset $\mathcal{I}_{<r}(\boldsymbol{b}_{<r})$ that are
combinations of the operators $\boldsymbol{b}_{<r}$ and the set
$\mathcal{I}_{\geq r}(\boldsymbol{b}_{<r})$ of operators that are independent
of $\boldsymbol{b}_{<r}$. Choose the most biased operator
$b_{r}\in\mathcal{I}_{\geq r}(\boldsymbol{b}_{<r})$ and add it to the set
$\boldsymbol{b}_{<r}$. This gives
$\boldsymbol{b}_{<r+1}=(\boldsymbol{b}_{<r},b_{r})$. Iterate until $r=n$.

In order to show that this procedure generates the best independent model,
consider replacing $b_{r}$ with any other operator $b$. Since $\boldsymbol{b}$
is a complete basis for all operators, $b$ has to be expressed in terms of
them. With some abuse of notation, we write this as

$b=\prod_{a\in b}b_{a}=b_{r}b_{b-r}.$ (21)

The second relation expresses the fact that the operator $b_{r}$ necessarily
appears in this product, because the new basis $\boldsymbol{b}^{\prime}$ with
$b_{r}$ replaced by $b$, has to be complete. In Eq. (21) $b_{b-r}$ stands for
the product of the other operators in $\boldsymbol{b}$ excluding $b_{r}$, that
generate $b$.

It follows that the bias of the new operator

$m_{b}=\frac{1}{N}\sum_{i=1}^{N}b(\boldsymbol{s}^{(i)})=\frac{1}{N}\sum_{i=1}^{N}b_{r}(\boldsymbol{s}^{(i)})b_{b-r}(\boldsymbol{s}^{(i)})\leq
m_{r}.$ (22)

The last inequality derives from the fact that, by construction, $b_{r}$ is
the most biased operator among all those that contain $b_{r}$.

###  1.3 Sampling from an MCM

The probability distribution of a chosen MCM $\mathcal{M}$ at best fit reads:

$P(\boldsymbol{s}\,|\,\boldsymbol{\hat{g}},\mathcal{M})=\frac{1}{2^{n-r}}\prod_{a\in\mathcal{A}}\frac{k_{\boldsymbol{b}_{a}(\boldsymbol{s})}}{N},$
(23)

where $k_{\boldsymbol{b}_{a}}/N$ are the empirical probabilities of the states
of $\boldsymbol{b}_{a}$ in the data. To sample a configuration
$\boldsymbol{s}$ from this model, we draw a configuration of
$\boldsymbol{b}_{a}$ with probability
$q(\boldsymbol{b}_{a})=k_{\boldsymbol{b}_{a}}/N$ for each ICC
$\mathcal{M}_{a}$. This can be done by sampling directly from the empirical
distribution $q(\boldsymbol{b}_{a})$, or by drawing uniformly a configuration
$\boldsymbol{s}^{a}$ from the data and then computing the transformation
$\boldsymbol{b}_{a}(\boldsymbol{s}^{a})$ (Eq. (23) ensures that this procedure
produces a sample value of $\boldsymbol{b}_{a}$ with the correct
distribution). We then concatenate the $\boldsymbol{b}_{a}$ vectors to obtain
a value of the state vector
$\boldsymbol{b}=\cup_{a\in\mathcal{A}}\boldsymbol{b}_{a}$. The last $(n-r)$
spin states, which are not modeled by $\mathcal{M}$, are randomly sampled with
probability $1/2$ to be $+1$ or $-1$, and placed in the $(n-r)$-dimensional
vector $\bar{\boldsymbol{b}}$. Completing the basis
$\boldsymbol{b}(\boldsymbol{s})$ of $\mathcal{M}$ with $(n-r)$ independent
operators forms a basis of the whole $n$-dimensional spin space, which defines
a GT $\mathcal{T}$. Finally, we perform an inverse GT to obtain a sampled
value of
$\boldsymbol{s}=\mathcal{T}^{-1}(\boldsymbol{b}\cup\bar{\boldsymbol{b}})$.

This last part involves the inversion of an $n\times n$ binary matrix modulo
2. Indeed, spins can be represented as $s_{i}=(-1)^{x_{i}}$, with $x_{i}=0$ or
$1$. A spin configuration $\boldsymbol{s}=(s_{1},\cdots,s_{n})$ then
corresponds to a string of $n$ bits $\boldsymbol{x}=(x_{1},\cdots,x_{n})$. A
basis configuration $\boldsymbol{b}=(b_{1},\cdots,b_{n})$ can also be
represented in terms of a bit string $\boldsymbol{y}=(y_{1},\cdots,y_{n})$,
with $b_{j}=(-1)^{y_{j}}$. Then the GT $\boldsymbol{b}(\boldsymbol{s})$
corresponds to the matrix product:

$\boldsymbol{y}=\boldsymbol{x}\cdot\hat{\boldsymbol{T}}\,,$ (24)

where $\boldsymbol{x}$ and $\boldsymbol{y}$ are row vectors, and where
summations are performed modulo 2. Here, the elements of the matrix
$\hat{\boldsymbol{T}}$ are such that $T_{i,j}=1$ if the basis operator $b_{j}$
contains spin $s_{i}$, and $T_{i,j}=0$ otherwise. This way, one gets $y_{j}=1$
if and only if the number of negative spins that contribute to $b_{j}$ in
configuration $\boldsymbol{s}$ is odd. Therefore, to obtain the value of
$\boldsymbol{s}$ that corresponds to a value of the basis $\boldsymbol{b}$, we
need to invert the matrix $\hat{\boldsymbol{T}}$ and compute
$\boldsymbol{x}=\boldsymbol{y}\cdot\hat{\boldsymbol{T}}^{-1}$. This inversion
step can be done once, by Gaussian elimination. In summary, sampling from the
best MCM requires drawing $A$ independent configurations
$\boldsymbol{s}^{(i_{a})}$ from the dataset $\boldsymbol{\hat{s}}$ and a
matrix multiplication, modulo 2.

## Acknowledgements

We gratefully acknowledge Edward D. Lee and William Bialek for sharing the
data of [8]. We are grateful to the authors of Ref. [32] for sharing their
data. We are grateful to Iacopo Mastromatteo, Vijay Balasubramanian, and
Yasser Roudi for insightful discussions.

## Supplementary Materials

###  1.4 Background: useful results and definitions

This section provides a summary of definitions and results from Ref. [6] that
are used in the main text and in this document. We refer the reader to Ref.
[6] for comments and proofs of the results recalled in this section.

#### Spin operator and spin model

A spin model is a probabilistic model that describes the state of a system of
binary variables, called spins. The model assumes the existence of
interactions between the spins that constrain the states accessible to the
system. As no spatial organisation is assumed, interactions can be of
arbitrary order and of arbitrary range. To mathematically define a spin model,
each interaction of the model is associated with a spin operators.

More precisely, consider a system of $n$ spin variables,
$\boldsymbol{s}=(s_{1},\,\cdots,s_{n})$, that take random binary values
$s_{i}=\pm 1$. A spin operator $\phi^{\boldsymbol{\mu}}(\boldsymbol{s})$
associated with an interaction labeled by $\boldsymbol{\mu}$ is defined as the
product of the spin $s_{i}$ involved in the interaction:

$\phi^{\boldsymbol{\mu}}(\boldsymbol{s})=\prod_{i\in\boldsymbol{\mu}}\,s_{i}\,.$
(25)

In practice, $\boldsymbol{\mu}$ can be taken as an $n$-dimensional binary
vector encoding which spin is involved in the interaction (i.e., its $i$-th
element $\mu_{i}=1$ if spin $s_{i}$ is included in the interaction, and
$\mu_{i}=0$ otherwise). The product in Eq. (25) is then over all the indices
$i$ such that the $i$-th element of $\boldsymbol{\mu}$ is $1$. The number of
spin operators is therefore the number of non-null such binary vector
$\boldsymbol{\mu}$, which is $2^{n}-1$. With the addition of the identity
operator, $\phi^{\boldsymbol{0}}(\boldsymbol{s})=1$ for all $\boldsymbol{s}$,
the set of all operators $\\{\phi^{\boldsymbol{\mu}}\\}$ forms a group, in the
sense that the product
$\phi^{\boldsymbol{\mu}}\phi^{\boldsymbol{\nu}}=\phi^{\boldsymbol{\mu}\oplus\boldsymbol{\nu}}$
of any two operators is an operator131313The operator $\oplus$ is the bitwise
XOR Boolean operator applied to the two binary vectors $\boldsymbol{\mu}$ and
$\boldsymbol{\nu}$. . Notice that the square of any spin variable $s_{i}$ and
of any operator is equal to one. The set of operators forms also a complete
and orthogonal basis for all functions of $n$ binary variables, as

$\sum_{\boldsymbol{s}}\phi^{\boldsymbol{\mu}}(\boldsymbol{s})\,\phi^{\boldsymbol{\nu}}(\boldsymbol{s})=2^{n}\,\delta_{\boldsymbol{\mu},\boldsymbol{\nu}}\qquad{\rm
and}\qquad\sum_{\boldsymbol{\mu}=0}^{2^{n}-1}\phi^{\boldsymbol{\mu}}(\boldsymbol{s})\,\phi^{\boldsymbol{\mu}}(\boldsymbol{s}^{\prime})=2^{n}\,\delta_{\boldsymbol{s},\boldsymbol{s}^{\prime}}.$
(26)

Hence any function $F(\boldsymbol{s})=\sum_{\boldsymbol{\mu}\geq
0}f^{\boldsymbol{\mu}}\phi^{\boldsymbol{\mu}}(\boldsymbol{s})$ defined on the
spin configurations can be represented as a linear combination of spin
operators, with coefficients given by
$f^{\boldsymbol{\mu}}=2^{-(n-1)}\sum_{\boldsymbol{s}}\phi^{\boldsymbol{\mu}}(\boldsymbol{s})F(\boldsymbol{s})$.
These completeness relations hold also on the sub-space defined by a subset of
the spins, or by a set of independent operators.

A spin model is a maximum entropy model, whose Hamiltonian is a linear
combination of the elements of a set $\mathcal{M}$ of spin operators. The
probability distribution over the spin variables $\boldsymbol{s}$ under a spin
model $\mathcal{M}$ is therefore:

$\displaystyle
P(\boldsymbol{s}\,|\,\boldsymbol{g},\,\mathcal{M})=\frac{1}{Z_{\mathcal{M}}(\boldsymbol{g})}\,\exp\bigg{(}\sum_{\boldsymbol{\mu}\in\mathcal{M}}g_{\mu}\,\phi^{\boldsymbol{\mu}}(\boldsymbol{s})\bigg{)}\,,\qquad{\rm
with}\qquad
Z_{\mathcal{M}}(\boldsymbol{g})=\sum_{\boldsymbol{s}}\,\exp\left({\sum_{\boldsymbol{\mu}\in\mathcal{M}}g_{\boldsymbol{\mu}}\,\phi^{\boldsymbol{\mu}}(\boldsymbol{s})}\right)\,,$
(27)

where
$\boldsymbol{g}=\\{g_{\boldsymbol{\mu}},\,\boldsymbol{\mu}\in\mathcal{M}\\}$
is a vector of real parameters and where the partition function
$Z_{\mathcal{M}}(\boldsymbol{g})$ ensures normalisation. Each parameter
$g_{\boldsymbol{\mu}}$ modulates the strength of the interaction associated
with the operator $\phi^{\boldsymbol{\mu}}(\boldsymbol{s})$. In an $n$-spin
system, there are $2^{n}-1$ different spin operators, hence there are
$2^{2^{n}-1}$ different spin models that can be constructed. Each model is
identifies with the set $\mathcal{M}$ of spin operators that it contains.

#### Set of independent operators

The operators $\phi^{\mu}$ in a set $\mathcal{I}$ are said independent if none
of the operators of the set can be obtained as a product of any subset of
other operators of the set $\mathcal{I}$. For example, the sets
$\\{s_{1},\,s_{2},\,s_{3}\\}$ and $\\{s_{1},\,s_{1}s_{2},\,s_{1}s_{2}s_{3}\\}$
are two examples of such set, whereas $\\{s_{1},\,s_{2},\,s_{1}s_{2}\\}$ is a
counter-example, because any operator in this set equals the product of the
two other operators. In an $n$-spin system, the set of the $n$ spin variables
$\mathcal{I}_{s}=\\{s_{1},\ldots,s_{n}\\}$ is set of independent operators.
All the operators of the system are generated as products of the $n$ basis
spin elements $(s_{1},\cdots,s_{n})$. Since there is no other independent
operator, a set of independent operators can have at maximum $n$ elements.

#### Gauge transformation

If
$\mathcal{I}=\\{\phi_{1}(\boldsymbol{s}),\ldots,\phi_{n}(\boldsymbol{s})\\}$
is a set of $n$ independent operators, the transformation
$\boldsymbol{s}\to\boldsymbol{s}^{\prime}$ where
$s_{i}^{\prime}(\boldsymbol{s})=\phi_{i}(\boldsymbol{s})$ establishes a one to
one mapping between the set of states onto itself. This also establishes a
bijection of the set of operators onto itself and of the set of models.
Following Ref. [6], we call this a gauge transformations (GTs). A GT can be
thought of as a change of basis from the original spin representation to a new
one. Figure 3 of the main text shows some examples of gauge transformations.

Ref. [6] shows that a GT leaves the partition function
$Z_{\mathcal{M}}(\boldsymbol{g})$ of a model invariant up to permutation of
its parameters. Likewise, the Fisher Information matrix

$I_{\mu,\nu}(\boldsymbol{g})=\partial_{g_{\mu}}\partial_{g_{\nu}}\log
Z(\boldsymbol{g})$ (28)

enjoys the same invariance property. This shows that the model complexity

$c_{\mathcal{M}}=\int d\boldsymbol{g}\sqrt{{\rm det}I(\boldsymbol{g})}$ (29)

is invariant under GTs. This allows to classify all models into equivalence
classes, characterised by the same complexity $c_{\mathcal{M}}$ (see Ref.
[6]).

#### Independent Models

An Independent Model (IM) is a spin model defined by a set of $r$ independent
spin operators, $\mathcal{M}_{ind}=\\{\phi_{1},\ldots,\phi_{r}\\}$, where
necessarily $r\leq n$. A model with $r$ different single-body interactions
$\mathcal{M}_{ind}=\\{s_{i_{1}},\ldots,s_{i_{r}}\\}$ is a straightforward
example of an independent spin model. All the IM with $r$ interactions can be
obtained by gauge transformations of this model. Among all models with $r$
operators, independent models achieve the maximal value of the complexity,
which is given by $r\log\pi$ [6].

#### Independent Complete Component Models

An Independent Complete Component (ICC) model $\mathcal{M}$ is such that for
any two operators of $\mathcal{M}$,
$\boldsymbol{\mu},\boldsymbol{\nu}\in\mathcal{M}$, the product operator
$\boldsymbol{\mu}\oplus\boldsymbol{\nu}$ is also an operator of $\mathcal{M}$.
An ICC model has $2^{r}-1$ operators, where the rank $r$ is the maximal number
of independent operators in $\mathcal{M}$. In other words, an ICC is a model
with all the $2^{r}-1$ operators generated by a basis of $r$ independent
operators, where $r\leq n$. The ICC model with $r=n$ is called the complete
model. A GT transforms an ICC model into another ICC model (see Fig. 3 for
examples), which therefore has the same complexity.

#### Minimally Complex Models

A MCM is a model composed of independent complete components (ICC):

$\mathcal{M}=\bigcup_{a\in\mathcal{A}}\mathcal{M}_{a},$ (30)

where each $\mathcal{M}_{a}$ is an ICC with rank $r_{a}$, and
$\mathcal{M}_{a}\bigcap\mathcal{M}_{a^{\prime}}=\emptyset$ for all $a\neq
a^{\prime}\in\mathcal{A}$. A GT maps the family of MCMs into itself,
preserving the list of ranks $r_{a}$ of each ICC.

###  1.5 Enumeration

#### Number of independent models

All IM with $r$ interactions can be obtained by gauge transformations of the
model $\mathcal{I}=\\{s_{1},\ldots,s_{r}\\}$. Therefore their number equals
the number of GTs that transforms this model into different ones, and it is
given by:

$\mathcal{N}_{ind}(n,r)=\frac{1}{r!}\prod_{i=0}^{r-1}\left(2^{n}-2^{i}\right).$
(31)

This corresponds to the number of ways of choosing $r$ independent operators
in the set of all possible $2^{n}-1$ operators, divided by the number of
possible permutations of these $r$ elements. The first operator can be chosen
in $2^{n}-1$ ways. After choosing $i$ independent operators, the $i+1^{\rm
st}$ operator can be chosen among $2^{n}-2^{i}$ ones, because $2^{i}$
operators are dependent on the first $i$ operators. For large values of $n$,
the number of IM of rank $r$ grows as $2^{nr}/r!\,$. Summing over all possible
values of $r$ (from $r=0$ to $n$), we evaluate that the total number of
independent models grows roughly as $2^{n^{2}}$ (see Fig. 2 of the main text),
which is faster than the number of pairwise models in the original basis $\sim
2^{n(n+1)/2}$.

#### Number of Independent Complete Components of rank $r$

In an $n$-spin system, all the ICCs of a given rank $r$ can be obtained by
taking all possible GTs of any chosen ICC $\mathcal{M}$ of rank $r$. As the
GTs that involve only the operators in $\mathcal{M}$ leaves $\mathcal{M}$
invariant, the number of ICCs is therefore given by

$\displaystyle\mathcal{N}_{\rm
ICC}(n,r)=\prod_{i=0}^{r-1}\frac{2^{n}-2^{i}}{2^{r}-2^{i}}\,.$ (32)

Here the numerator $\prod_{i=0}^{r-1}(2^{n}-2^{i})$ counts the number of ways
of choosing the $r$ basis operators of the ICC. The denominator
$\prod_{i=0}^{r-1}(2^{r}-2^{i})$ counts the number of GT that transforms the
basis of the ICC while leaving $\mathcal{M}$ invariant. The total number of
models with a single ICC in an $n$-spin system is obtained by summing Eq. (32)
over all possible values of $r$ ($r=0$ to $n$) and is displayed in Fig. 2 of
the main text (see curve labelled “SCM”).

#### Number of Minimally Complex Models

For an MCM on $n$ spins with $m_{r_{a}}$ ICCs of same rank ${r_{a}}$ (for all
${r_{a}}$ from $1$ to $n$), one can combined the two previous arguments and
get:

$\displaystyle\mathcal{N}_{\rm
MCM}(n,\\{m_{r_{a}}\\})=\frac{\displaystyle\prod_{i=0}^{r-1}\left(2^{n}-2^{i}\right)}{\displaystyle\prod_{r_{a}=1}^{n}m_{r_{a}}!\,\left[\prod_{i=0}^{r_{a}-1}(2^{r_{a}}-2^{i})\right]^{m_{r_{a}}}}$
(33)

where the factor $1/m_{r_{a}}!$ accounts for permutations among ICC of the
same size. In order to count the total number of MCMs with $n$ spins (as shown
in Fig. 2 of the main text), it is necessary to sum over all the classes of
MCMs with different degeneracies $m_{r_{a}}$. The number of classes with
different degeneracy corresponds to the number of partitions of the rank $r$
of the MCM in the form $r=\sum_{r_{a}}{r_{a}}m_{r_{a}}$, and summing over all
$r\leq n$. For instance, $r=8$ admits $22$ partitions which are:

$\displaystyle 8,71,62,611,53,521,5111,44,431,422,4211,41111,332,3311,$
$\displaystyle 3221,32111,311111,2222,22211,221111,2111111,11111111\,.$ (34)

The partition “$8$” corresponds to the class of MCMs that contains only one
single ICC with rank $r=8$. The partition “$422$” is the class of MCMs formed
of three ICC models, with $m_{r_{a}}=0$ except for $m_{2}=2$ and $m_{4}=1$.

#### Number of Minimally Complex Models with the same preferred basis

For a given choice of $n$ basis operators, the number of MCMs of rank $r=n$
that admits that basis as a preferred basis is given by the number of possible
partitions of this set of $n$ basis elements, which is given by the Bell
number $B_{n}$. Ref. [21] gives an upper-bound for the Bell number,
$B_{n}<\big{(}\frac{0.792\,n}{\ln(n+1)}\big{)}^{n}$. The number of MCM∗
(indicated in Fig. 2 of the main text) corresponds the total number of MCMs
that share the same preferred basis, i.e. the MCMs of any rank from $r=0$ to
$r=n$. This number is given by:

$\displaystyle\mathcal{N}_{MCM^{*}}(n)=\sum_{r=0}^{n}\binom{n}{r}\,B_{r}\,.$
(35)

For example, for $n=9$ there are $115975$ MCMs in total that shares the same
preferred basis (i.e., MCM∗), $21147$ of them have a rank $r=9$.

Figure 8: (colors online) Complexity $c_{\mathcal{M}}$ of all spin models with
$n=4$ spins as a function of the number of parameters $K$ in the model.
Adaptation of Fig. 4 from Ref. [6], where we marked with a circle the classes
of MCM. Models of the same class have the same complexity and are represented
by either a triangle or a circle (see Ref. [6] for more details). Model
representants of some of the classes are displayed using the same diagrammatic
notation as in Fig. 1 of the main text. The complexity value of any model lies
between two limit lines: one corresponds to the complexity of IM (solid line),
the other to the complexity of SCM (dashed line). The colors indicated the
number $r$ of independent variables in the model (rank): red for 0, pink for
1, green for 3 and violet for 4. Among all models with same number $K$ of
parameters and same rank $r$, MCMs are the simplest.  Figure 9: (colors
online) Distribution of the log-evidence of MCMs for the US Supreme court
data: for all the MCMs that admits the best IM $\boldsymbol{b}^{*}$ as a
preferred basis (top graph, in light blue) and for all the MCMs that admits
the original basis of the data as a preferred basis (bottom graph, in light
blue). The best MCMs overall is based on $\boldsymbol{b}^{*}$ and is shown on
the top graph; it has three ICCs, each highlighted by a dotted red, blue or
orange square, and has a log-evidence of $\log E\simeq-3327$ (indicated by the
vertical red line on the top graph). The best MCM based on the original basis
of the justices is shown on the bottom graph; it has two ICCs, each
highlighted by a dotted blue or red square and a log-evidence of $\log
E\simeq-5258$ (indicated by the vertical blue line on the bottom graph). In
both graphs, the distribution in blue corresponds to the distribution of all
the MCMs with a chosen preferred basis: with $n=9$ there are $115975$ such
models (models of all rank from $r=0$ to $r=9$). In both graphs, a part of
this distribution has been highlighted in red; it corresponds to the MCMs with
rank $r=9$: with $n=9$ there are $21147$ such models. It is interesting to
observe that most MCMs based on the original basis (blue distribution in the
bottom graph) have a significantly lower log-evidence than the MCMs of rank
$r=9$ based on the best IM $\boldsymbol{b}^{*}$ (red distribution in the top
graph).  Figure 10: (colors online) Distribution of the log-evidence of MCMs
with 5 ICCs for the Big 5 data in its original basis: $10^{4}$ partitions of
the $n=50$ questions into $|\mathcal{A}|=5$ parts (ICCs) were sampled
uniformly at random. The log-evidence of the MCM partitioning the questions
into the big 5 personalities ($\log E\simeq-2.518\cdot 10^{7}$, indicated by a
vertical red line) is larger than the best log-evidence found by our heuristic
algorithm ($\log E\simeq-2.596\cdot 10^{7}$, vertical black line). However,
both values of the log-evidence are to the right of the distribution, deep in
its tail. There are $\binom{50}{5}=2\,118\,760$ possible partitions of the
$50$ questions in $5$ parts. Among the $10^{4}$ randomly sampled partitions,
the largest evidence found was $\log E\simeq-2.700\cdot 10^{7}$ and the lowest
evidence was $\log E\simeq-3.399\cdot 10^{7}$.

## References

  * [1] George EP Box.  Science and statistics.  Journal of the American Statistical Association, 71(356):791–799, 1976.
  * [2] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  nature, 521(7553):436–444, 2015.
  * [3] Eugene P Wigner.  The unreasonable effectiveness of mathematics in the natural sciences.  In Mathematics and Science, pages 291–306. World Scientific, 1990\.
  * [4] Peter D Grünwald and Abhijit Grunwald.  The minimum description length principle.  MIT press, 2007.
  * [5] In Jae Myung, Vijay Balasubramanian, and Mark A Pitt.  Counting probability distributions: Differential geometry and model selection.  Proceedings of the National Academy of Sciences, 97(21):11170–11175, 2000.
  * [6] Alberto Beretta, Claudia Battistin, Clélia De Mulatier, Iacopo Mastromatteo, and Matteo Marsili.  The stochastic complexity of spin models: Are pairwise models really simple?  Entropy, 20(10):739, 2018.
  * [7] Elad Schneidman, Michael J Berry II, Ronen Segev, and William Bialek.  Weak pairwise correlations imply strongly correlated network states in a neural population.  Nature, 440(7087):1007, 2006.
  * [8] Cristina Savin and Gašper Tkačik.  Maximum entropy models as a tool for building precise neural controls.  Current opinion in neurobiology, 46:120–126, 2017.
  * [9] Edward D Lee, Chase P Broedersz, and William Bialek.  Statistical mechanics of the us supreme court.  Journal of Statistical Physics, 160(2):275–301, 2015.
  * [10] Andrea Montanari and Jose A Pereira.  Which graphical models are difficult to learn?  In Advances in Neural Information Processing Systems, pages 1303–1311, 2009.
  * [11] Simona Cocco and Rémi Monasson.  Adaptive cluster expansion for the inverse ising problem: convergence, algorithm and tests.  Journal of Statistical Physics, 147(2):252–314, 2012.
  * [12] John P Barton, Eleonora De Leonardis, Alice Coucke, and Simona Cocco.  Ace: adaptive cluster expansion for maximum entropy graphical model inference.  Bioinformatics, 32(20):3089–3097, 2016.
  * [13] Marc Vuffray, Sidhant Misra, Andrey Lokhov, and Michael Chertkov.  Interaction screening: Efficient and sample-optimal learning of ising models.  In Advances in Neural Information Processing Systems, pages 2595–2603, 2016.
  * [14] H Chau Nguyen, Riccardo Zecchina, and Johannes Berg.  Inverse statistical problems: from the inverse ising problem to data science.  Advances in Physics, 66(3):197–261, 2017.
  * [15] Iacopo Mastromatteo.  On the typical properties of inverse problems in statistical mechanics.  arXiv preprint arXiv:1311.0190, 2013.
  * [16] John Kruschke.  Doing Bayesian data analysis: A tutorial with R, JAGS, and Stan.  Academic Press, 2014.
  * [17] Harold Jeffreys.  An invariant form for the prior probability in estimation problems.  Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 186(1007):453–461, 1946.
  * [18] Jorma Rissanen.  Stochastic complexity and modeling.  The annals of statistics, pages 1080–1100, 1986.
  * [19] Jorma J Rissanen.  Fisher information and stochastic complexity.  IEEE transactions on information theory, 42(1):40–47, 1996.
  * [20] Clélia de Mulatier.  Mincompspin.  https://github.com/clelidm/MinCompSpin, 2021.
  * [21] Daniel Berend and Tamir Tassa.  Improved bounds on bell numbers and on moments of sums of random variables.  Probability and Mathematical Statistics, 30(2):185–205, 2010.
  * [22] H.J. Spaeth, L. Epstein, T.W. Ruger, K. Whittington, J.A. Segal, and A.D. Martin.  Supreme court database – version 2011 release 3.  2011\.
  * [23] Luigi Gresele and Matteo Marsili.  On maximum entropy and inference.  Entropy, 19(12), 2017.
  * [24] Gideon Ehrlich.  Loopless algorithms for generating permutations, combinations, and other combinatorial configurations.  Journal of the ACM (JACM), 20(3):500–513, 1973.
  * [25] Donald E. Knuth.  The Art of Computer Programming, Volume 4, Combinatorial Algorithms: Part 1.  Addison-Wesley Professional, 2011.
  * [26] B. Djokić, M. Miyakawa, S. Sekiguchi, I. Semba, and I. Stojmenovi?  Short Note: A Fast Iterative Algorithm for Generating Set Partitions.  The Computer Journal, 32(3):281–282, 01 1989.
  * [27] Lewis R Goldberg.  The development of markers for the big-five factor structure.  Psychological assessment, 4(1):26, 1992.
  * [28] Bojan Tunguz.  Big five personality test.  https://openpsychometrics.org/_rawdata/IPIP-FFM-data-8Nov2018.zip.  Version 1. Created: 2020-02-17. Accessed: 2020-06-20.
  * [29] Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324, 1998.
  * [30] Leenoy Meshulam, Jeffrey L Gauthier, Carlos D Brody, David W Tank, and William Bialek.  Coarse graining, fixed points, and scaling in a large population of neurons.  Physical review letters, 123(17):178103, 2019.
  * [31] Pradeep Ravikumar, Martin J Wainwright, John D Lafferty, et al.  High-dimensional ising model selection using l1-regularized logistic regression.  The Annals of Statistics, 38(3):1287–1319, 2010.
  * [32] Hanne Stensola, Tor Stensola, Trygve Solstad, Kristian Frøland, May-Britt Moser, and Edvard I Moser.  The entorhinal grid map is discretized.  Nature, 492(7427):72–78, 2012.