Why the simplest explanation is always the best

Preamble

The simplest explanation is always the best among the explanations that are representative.  It is called the principle of parsimony,  Occam's razor. It is the bedrock of scientific enlightenment. There is a recent development that people start to do a category error to discard this principle based on the following setting that lead to misunderstanding: if we say that the simplest model that can explain the data in one representation but another representation requires a more complex model, then we are choosing more complex model. This is obviously wrong. The simplest explanation is to be chosen from the models that captures the complexity on the given representations, not the simplest over both representations. Filter explanations first based on representations then selection follows. Here we shown the core idea via an illustrative example. 

Figure: Circle has a zero Pearson
correlation. (Wikipedia)

Revisiting Occam’s Razor: A case of correlation and geometry

In order to understand this category error, we will work on a concrete example. Let's say we have a data, which is a circular shape $\mathscr{D}(x, y)$, as in figure. We have 3 models: 

$$\mathscr{M}_{1} : y = a x + b $$

$$\mathscr{M}_{2} : y = \sqrt{1-x^{2}} $$

$$\mathscr{M}_{3} : y ~ NN(x)  $$

And $\mathscr{U}$ a utility function, Pearson correlation $C(x, y)$.  We consider performance measure in ranking for the utility and a measure of representation. $\mathscr{M}_{3}$ is a neural network with a lot of parameters.

There is no error in choosing $\mathscr{M}_{1}$ based on the similar correlations these models produce. Principle of parsimony is not violated at all.  This is correct if we are only considering numerical representation. Pearson correlation as a utility works well for purely numerical representation. What about geometric representation? Then we need to change our utility function (representation measure or performance function). 

Let's say if we use curvature as a utility, $\kappa(x,y)$, in this case $\mathscr{M}_{1}$ fails to capture curvature and its is filtered out before Occam's razor can be applied. Then we left with  $\mathscr{M}_{2}$ and $\mathscr{M}_{3}$ . 

Correlation and geometric explanations are two different things. Two vastly different geometries can produce the same correlations.  A model can be quite good in explaining correlation but fails to capture geometric complexity. In this setting, it does not mean that Occam's razor is wrong. We need to apply Occam's razor across representations: numeric, geometric, algebraic, or symbolic, depending on the purpose. Always keep in mind the purpose or utility of model when invoking Occam's razor. Simplest explanation over the required representations that are relevant are the best. 

On the utility, performance and representations measure

Performance function, representation measure and utility functions can be different in real life, here for illustration purposes we consider them interchangeably. 

Conclusion

Nature minimises cost over complexity but under utility constraints. Minimal cost without satisfying utility or representation measure won't be chosen despite being the simplest. We need to filter first based on utility or representation measure before applying Occam's razor. The simplest explanation is always the best among the explanations that are representative. 

 Cite as 

 @misc{suzen25occam, 

     title = { Why the simplest explanation is always the best}, 

     howpublished = {\url{https://science-memo.blogspot.com/2025/11/simplest-explanation-always-best.html}}, 

     author = {Mehmet Süzen},

     year = {2025}

}  

Overfitting is about complexity ranking of inductive biases : Algorithmic recipe

Preamble

    Figure: Moon patterns
human brain
 invents. (Wikipedia)

Detecting overfitting is inherently a comparison problem of the complexity of multiple objects, i.e., models or an algorithm capable of making predictions. A model is overfitted (underfitted) if we only compare it to another model. Model selection involves comparing multiple models with different complexities. The summary of this approach with basic mathematical definitions is given here.

Misconceptions: Poor generalisation is not synonymous with overfitting. 

None of these techniques would prevent us from overfitting: Cross-validation, having more data, early stopping, and comparing test-train learning curves are all about generalisation. Their purpose is not to detect overfitting.

We need at least two different models, i.e., two different inductive biases, to judge which model is overfitted. One distinct approach in deep learning, called dropout, prevents overfitting while it alternates between multiple models, i.e., multiple inductive bias. For judgment, dropout implementation has to compare those alternating model test performances during training to judge overfitting. 

What is an inductive bias? 

There are multiple inceptions of inductive bias. Here, we concentrate on a parametrised model, $\mathscr{M}(\theta)$ on a dataset $\mathscr{D}$, the selection of a model type, or modelling approach, usually manifest as a functional form $\mathscr{M}=f(x)$ or as a function approximation, i.e., for example neural network, are all manifestation of inductive biases. Different parameterisation of model learned on the subsets of the dataset are still the same inductive bias.

Complexity ranking of inductive biases: An Algorithmic recipe 

We are sketching out an algorithmic recipe for complexity ranking of inductive biases via informal steps:

1. Define a complexity measure $\mathscr{C}$($\mathscr{M}$) over an inductive bias.
2. Define a generalisation measure  $\mathscr{G}$($\mathscr{M}$, $\mathscr{D}$) over and inductive bias and dataset.
3. Select a set of inductive biases, at least-two, $\mathscr{M}_{1}$ and $\mathscr{M}_{2}$.
4. Produce complexity and generalisation measures on ($\mathscr{M}$, $\mathscr{D}$): Here for two inductive biases: $\mathscr{C}_{1}$, $\mathscr{C}_{2}$,   $\mathscr{G}_{1}$, $\mathscr{G}_{2}$.
5. Ranking of  $\mathscr{M}_{1}$ and $\mathscr{M}_{2}$:  $argmax \{ \mathscr{G}_{1}, \mathscr{G}_{2}\}$ and $argmin \{ \mathscr{C}_{1}, \mathscr{C}_{2}\}$

The core concept appears as when generalisations are close enough we pick out the inductive bias that is less complex. 

Conclusion & Outlook

In practice,  probably due to hectic delivery constraints, or mere laziness, we still rely on simple holdout method to build models, only single test and train split, not even learning curves, specially in deep learning models without practicing Occam's razor. A major insight in this direction appears to be that, holdout approach can only help us to detect generalisation, not overfitting. We clarify this via the concept of inductive bias distinguishing that different parametrisation of the same model doesn't change the inductive bias introduced by the modelling choice. 

In fact, due to resource constraints of model life-cycle, i.e., energy consumption and cognitive load of introducing a complex model, practicing proper Occam's razor: complexity ranking of inductive biases, is much more important than ever for sustainable environment and human capital.

Further reading

Some of the posts, reverse chronological order, that this blog have tried to convey what overfitting entails and its general implications. 

• Empirical risk minimization is not learning :A mathematical definition of learning and re-understanding of overfitting and Occam's razor in machine learning
• Critical look on why deployed machine learning model performance degrade quickly.
• Bringing back Occam's razor to modern connectionist machine learning.
• Understanding overfitting: an inaccurate meme in Machine Learning

Critical look on why deployed machine learning model performance degrade quickly

Illustration of William of Ockham 
(Wikipedia)

One of the major problems in using so called machine learning model, usually a supervised model, in so called deployment, meaning it will serve new data points which were not in the training or test set,  with great astonishment, modellers or data scientist observe that model's performance degrade quickly or it doesn't perform as good as test set performance. We earlier ruled out that underspecification would not be the main cause. Here we proposed that the primary reason of such performance degradation lies on the usage of hold out method in judging generalised performance solely.

Why model test performance does not reflect in deployment? Understanding overfitting

Major contributing factor is due to inaccurate meme of overfitting which actually meant overtraining and connecting overtraining erroneously to generalisation solely.  This was discussed earlier here as understanding overfitting. Overfitting is not about how good  is the function approximation compared to other subsets of the dataset of the same “model” works. Hence, the hold-out method (test/train) of measuring performances  does not  provide sufficient and necessary conditions to judge model’s generalisation ability: with this approach we can not detect overfitting (in Occam’s razor sense) and as well the deployment performance. 

How to mimic deployment performance?

This depends on the use case but the most promising approaches lies in adaptive analysis and detected distribution shifts and build models accordingly. However, the answer to this question is still an open research.

Understanding the empirical law of large numbers and the gambler's fallacy

    Platonic dices (Wikipedia).

One of the misconceptions in our understanding of 

statistics, or a counter-intuitive guess, fallacy, appears in the assumption of the existence of the law of averages. Imagine we toss a fair coin many times, most people would think that the number of heads and tails would be balanced over the increasing number of trails, which is wrong. If you don't, then you might have a very good statistical intuition. Briefly,  we will illustrate this, a kind of gambler's fallacy with a simple simulation approach and discuss the empirical law of large numbers.

Empirical law of large numbers

If we repeat an experiment long enough, we would approach to expected outcome. The simplest example is a coin-toss experiment, that an expected 
fair coin toss would lead to equal likelihood for head and tail, 1 or 0.  This implies that, the ratio of head and tails will approach to one with increasing number of repeats. Let's say, we toss the coin $N$ times. The number of  heads and tails would be $n_{1}$ and $n_{0}$. The emprical law of large numbers states

$ \frac{n_{1}}{n_{0}} \to 1.0 $ if $N$ is very large, $N >>1$


But note that, the absolute difference, $|n_{1}-n_{0}|$ does not approach to any constant, on the contrary, it will increase with increasing number of repeats. This is classic example of gambler's fallacy that an outcome would balance out as there are more repeats. 

Fair coin-toss: Evolution of Bernoulli process

The Bernoulli process expresses binary outcomes, 1 or 0, i.e., success or failure, true or false. Bernoulli distribution reads,

$Ber(p, k) = p^{k} (1-p)^{1-k}$ for $k \in \{0,1\}$.

Figure 2: No law of averages. Absolute difference
of occurances
increases over repeats.

Figure 1: Empirical law of large numbers, ratio of
occurances approach to a constant.

$p$ is the probability of success. We draw 50K samples from this distribution to get a Bernoulli process with $p=0.5$ and repeat the experiment 50 times, in order to obtain a "generalised" behaviour with uncertainty. This situation corresponds to a fair coin-toss experiment.  

Results


Empirical results of ratio of two outcomes and their absolute difference over repeats are reported in Figure 1 and 2 respectively..



Appendix: Source codes

R and Rcpp functions are shown in this section to reproduce the plots in this post. Source files are also available on github (here).

#'
#'  Misconception of the law of averages: Understanding the empirical law of large numbers
#' 
#'  (c) 2016 Mehmet Suzen 
#'  

rm(list=ls())
set.seed(987231)
library(Rcpp)          # law level implementations
library(LaplacesDemon) # for drawing samples from Bernoulli distribution
library(matlab)        # for tic/toc
# Set working directory to script directory 
Rcpp:::sourceCpp("empirical_law_of_large_numbers.cpp") # diff_binary_vec2
#'
#' Calculate the running difference in a binary vector 
#' 
diff_binary_vec  <- function(bvec) {
 ll <- length(bvec) 
 diff_vec  <- vector("double", ll)
 diff_vec  <- sapply(1:ll, function(i) { abs(i-2.0*sum(bvec[1:i])) } )
 diff_vec
}
#'
#' Calculate the running ratio in a binary vector 
#'
ratio_binary_vec  <- function(bvec) {
 ll         <- length(bvec) 
 ratio_vec  <- vector("double", ll)
 ratio_vec  <- sapply(1:ll, function(i) { abs(i/sum(bvec[1:i])-1.0) } )
 ratio_vec  <- sapply(ratio_vec, function(rv) {if(is.infinite(rv)) { rv <- 0; }; rv })
 ratio_vec
}

#' Wall-call timing difference 
tb <- rbern(20000, 0.5) # a fair-coin
tic()
t1 <- diff_binary_vec(tb)
toc()
tic()
t2 <- diff_binary_vec2(tb)
toc()
tic()
r1 <- ratio_binary_vec(tb)
toc()
tic()
r2 <- ratio_binary_vec2(tb)
toc()
#' 
#' Generate Bernoulli Process
#' 
nr        <- 50    # repeats
nt        <- 20000 # Bernoulli trails
tic()
bern_df <- data.frame(trail=c(), diff=c(), ratio=c())
for(i in 1:nr) {
  cat("repeat:",i, "\n")
  trail   <- rbern(nt, 0.5) # a fair-coin
  diff    <- diff_binary_vec2(trail)
  ratio   <- ratio_binary_vec2(trail)
  bern_df <- rbind(bern_df, cbind(1:nt, diff, ratio))
}

#' Now plot ratio and diff evolution with local regression in ggplot2
names(bern_df) <- c("trail", "diff", "ratio")
library(ggplot2)
p_diff <- ggplot(data=bern_df, aes(x=trail, y=diff)) + geom_smooth(formula="y~x") + 
          theme(
                panel.background = element_blank(), 
                axis.text.x      = element_text(face="bold", color="#000000", size=11),
                axis.text.y      = element_text(face="bold", color="#000000", size=11),
                axis.title.x     = element_text(face="bold", color="#000000", size=11),
                axis.title.y     = element_text(face="bold", color="#000000", size=11))  +
           xlab("Bernoulli Trails") + ylab("Difference between occurance of two outcomes") +
           ggtitle("No Law of Averages: Tail/Heads do not balanced out!")

png(file="no_law_of_averages.png")
p_diff
dev.off()

p_ratio <- ggplot(data=bern_df, aes(x=trail, y=ratio)) + geom_smooth(formula="y~x") + 
          theme(
                panel.background = element_blank(), 
                axis.text.x      = element_text(face="bold", color="#000000", size=11),
                axis.text.y      = element_text(face="bold", color="#000000", size=11),
                axis.title.x     = element_text(face="bold", color="#000000", size=11),
                axis.title.y     = element_text(face="bold", color="#000000", size=11))  +
           xlab("Bernoulli Trails") + ylab("Difference between occurance of two outcomes") +
           ggtitle("Empirical Law of Large Numbers: Ratio of Tails/Heads approach to one")

png(file="law_of_large_numbers.png")
p_ratio
dev.off()

#include <Rcpp.h>
#include <stdlib.h>
using namespace Rcpp;

//
// Sum a numeric vector up to an index
//
double sum_nv(NumericVector x, int s) {
  int i=0;
  double ss=0.0;
  for(i=0;i<(s+1);i++) {
   ss=ss+x[i];
  }
  return(ss);
}

//
// Calculate the running difference in a Bernoulli process (binary_vec)
//
// [[Rcpp::export]]
//
NumericVector diff_binary_vec2(NumericVector binary_vec) {
  int ll = binary_vec.size();
  int i  = 0;
  NumericVector diff_vec(ll, 0.0);
  for(i=0;i<ll;i++) {
    diff_vec[i] = std::abs(i+1.0-2.0*sum_nv(binary_vec, i));
  }
  return(diff_vec);
}

//
// Calculate the running ratio in a Bernoulli process (binary_vec)
//
// [[Rcpp::export]]
NumericVector ratio_binary_vec2(NumericVector binary_vec) {
  int ll = binary_vec.size();
  int i  = 0;
  NumericVector ratio_vec(ll, 0.0);
  for(i=0;i<ll;i++) {
    ratio_vec[i] = std::abs(((i+1.0)/sum_nv(binary_vec, i))-1.0);
  }
  return(ratio_vec);
}

Football Goal Distributions and Statistical Mechanics

With the advent of complex network science and its allied approaches in the last decade or so with the data driven research, using statistical mechanical techniques out side of materials or atomic physics becomes a standard and quite a popular practice. One of the most interesting of this usage was in football (or soccer for Americans). This is probably because of mass football-mania in the UK and rest of the Europe and the new-world of course. However, using statistics is not a new thing at all but finding similarities with the atomic systems. Specially quantitative approach to sports is very well known, a recent film featuring Brad Pitt, Moneyball, shows this. Every major sports club (merchandise in the US) is now running a statistics division, considering the sky high salaries of players. There are interesting works in the goal statistics, being non-Gaussian [link] [link] [link] and passing network for football strategies [link] among other works. This kind of research is classified as econophysics.

Imputation of missing data: Recursive 1D discrete KNN algorithm

Any generated data is often have missing components or values. Probably, the most common occurrence manifest in time series data where there is no value available on the given time point, hence a NaN is placed in general (or NA in R). There is a large literature on how a statistical analysis must be performed in such data sets. For example, a seminal book by Little & Rubin called Statistical Analysis With Missing Data provides very detailed exposure to the field. Probably the simplest of all methods is called imputation. For example there are high quality R packages like imputation or mi that does the job for you. Similarly knnimpute from MATLAB bioinformatics toolbox provides similar solution.
  
k-nearest neighbour algorithm (KNN) is the most common approach to discover the closest available value in the data vector. However, often implementations of KNN contains a lot of options that are not needed for simple imputations in 1D and Euclidean metric. Here I propose 1D discrete KNN recursive algorithm that scans a given vector and determines the closest available value to given index). The main idea is assuming periodic boundary conditions for the vector index boundaries (Two-way linear search on the ring, see simple sketch). This is an $O(\bf{N})$ algorithm, we scan the vector twice in the worst case scenario. I have implemented this idea in MATLAB, for a given matrix. Files are available on the matworks file exchange [link]. This tool could be useful in imputing data in regression design matrices.
  
   

Hands on Crash Tutorial for Supervised Learning : Multinomial Logistic Regression (MNL)

If you are a computer scientist, probably you would call a task supervised learning what others might call classification based on data. Data being anything that has a well defined structure and associated classes. We will work on hands on example for multi-class prediction with logistic regression i.e. multinomial logistic regression (MNL). What I meant by multi-class is here that we have $k \in \bf{Z}$ distinct classes for each observation, $k>1$. Actually if you think in terms of link functions from Generalized Linear Models (GLMs), the support of the distribution will tell you the distinction of the nature of the class. In statistics literature classes can be manifest as factors (or categories).

The following pointers would be helpful before you read further. From R perspective, I'd suggest German Rodrigez's notes for background reading. Data mining blog by Will Dwinnell is one of the clear descriptions in this direction for MATLAB. For sure, an authoritative resource on this is the book by Hastie et. al. called elements of statistical learning, you can obtain it from their Stanford page. An other excellent resource from Professor Ripley of Oxford, his books and R packages are known to be de facto standard in the field, particularly here I refer to nnet and applied statistics book.

Recall that a design matrix $\bf{X}$ is noting more then set of observations ($n$) in the rows. Observables ($p$) are placed along columns. The idea of supervised learning is that we train our model, here we choose to be multinomial GLM, against a data set and obtain coefficients of the model (parameters).  Let's do this with the simplest possible example. (R codes start with '>' and MATLAB '>>')

Statistical toolbox brings set of functionality to do multinomial logistic regression. mnrfit is the main function we would like to use. Lets use a simple data set, it is trivial and the statistics we would get from this data might be very poor, but we aimed at the concept in this post.

>>  X  =  [0.0 0.1 0.7 1.0 1.1 1.3 1.4 1.7 2.1 2.2]'; % design matrix
>>  Y  =  [1 2 1 3 1 2 1 3 1 1]'; % associated classes 

We can split this set to training and validation data set by choosing

>> trainIndex = [2 8 10 4 5];
>> validIndex = [1 3  6 7 9]; 

Now we can obtain the coefficients via mnrfit.

>>  betaHat    = mnrfit(X(trainIndex), Y(trainIndex), 'model', 'ordinal', 'interactions', 'off', 'link', 'logit');
Note that, model is ordinal, meaning that it takes discrete values. Switched off interactions generates only one coefficient per class (betaHat vector).  Then we can get the probabilities for the validation set.

>> predictProbs=mnrval(betaHat, X(validIndex), 'model', 'ordinal', 'interactions', 'off', 'link', 'logit');

ans =

    0.2980    0.1958    0.5061
    0.3636    0.2041    0.4324
    0.4242    0.2045    0.3713
    0.4346    0.2039    0.3615
    0.5084    0.1954    0.2961

So the highest probabilites form this matrix will give us [3, 3, 1, 1, 1] predictive set. With this approach we have predicted last two classes correctly.

Let's do this same example with R, even though the following approach may not be exactly the same procedure described above, however this is multinomial as well:
> require(nnet)
> # Design Matrix: 10 observations
> X            = matrix(c(0.0, 0.1, 0.7, 1.0, 1.1, 1.3, 1.4, 1.7, 2.1, 2.2));
> # Corresponding Classes: 3 ordinal classes
> Y            =  matrix(c(0, 1, 0, 0,
                                     0, 0, 1, 0,
                                     0, 1, 0, 0,
                                     0, 0, 0, 1,
                                     0, 1, 0, 0,
                                     0, 0, 1, 0,
                                     0, 1, 0, 0,
                                     0, 0, 1, 0,
                                     0, 1, 0, 0,
                                     0, 1, 0, 0  
                                    ), nrow = 10, ncol=4, byrow=TRUE)
 
Similarly we can obtain training and validation sets

># Generate training and validation data sets for X, Y
> trainIndex   = c(2, 8, 10, 4, 5);
> validIndex   = c(1, 3, 6, 7, 9);
Now, we can use nnet package fitting function for multinomial log-linear.

> mfit = multinom(formula = Y[trainIndex, ] ~ X[trainIndex] + 0)

Note that we put 0 for the intercepts while first column of Y is dummy. 

> predict(mfit, X[validIndex])
[1] 2 2 2 2 2

Resulting classification is not that great. The above example shown a basic approach to supervised learning in MATLAB and R. However, one must check statistical values, like residuals, deviance, p-values etc. for the quality of results. 

Gaussians in n-dimensions

New algorithms in Gaussian Mixture Model (GMM) may sound quite oxymoron to analyze n-dimensional data sets but considering principle of parsimony it may be
the chosen one among mixture models.