A = B, B = C, but A ≠ C. Depending on N, of course!

Statistics are weird. But it all makes sense.

Let's say we have 3 variables and 600 measurements, in other words 600 data points with three dimensions or a 600-row matrix with three columns; or this chart:

These are simulated data, drawn from three Normal distributions with variance 1 and means $\mu_A = 0, \mu_B = 0.25,$ and $\mu_C = 0.5$. The distributions are:

There's considerable overlap between these distributions, so any single point is insufficient to determine any relationships between $A$, $B$, and $C$.

Let's say we want to have "99.9 percent confidence" in our assertions. What does that "99.9 percent confidence" mean? The statistical meaning is that there's at most 0.1 percent chance that if the data were generated by "the null hypothesis" (which in our case is that any two distributions are the same) we'd see a test statistic above the critical value.

Test statistic?! Critical value?!

Yes. A test statistic is something we'll compute from the data (a 'statistic') in order to test it. In our case it'll be the difference of the empirical, or sample, means divided by the standard error of those empirical means.

If that number, a summary of the difference between the samples, scaled to deal with the randomness, is greater than some value --- the critical value --- we trust that it's too big to be the result of the random variation. At least we trust it to be too big 99.9 percent of the time.

Okay, that statistics refresher aside, what can we tell from a sample of 25 points? Here are the distributions of those means:

Note how the distributions of the means are much narrower than the distributions of the data. (The means don't change.) That's the effect of averaging over 25 data points. The variance of the mean is $1/25$ and the standard deviation of the mean, called the standard error to avoid confusion with the standard deviation of the data, is $1/\sqrt{25}$.

Someone with a vague memory of a long-forgotten statistics class may recall seeing a $\sigma/\sqrt{n-1}$ in this context and try to argue that 25 should be a 24. And they'd be right if we were estimating the standard error of the population data from the standard error of the sample data; but we're not. Our data is simulated, therefore we know the standard error and we're using that to simplify things like this. Another one of which is the next one: the critical value.

Knowing the distributions lets us bypass one of the most overused jokes in statistics, the Student T ("how do you make it? Boil water, steep the student in it for 3 minutes"; student tea, get it?). More seriously, when the standard error is estimated from the sample data, the critical value is derived from a Student's T distribution; in our case we'll pick one derived from the Normal distribution, which has the advantage of not depending on the size of the sample (or, as it's called in statistics, the number of degrees of freedom in the estimation).

Now for the critical value. We're going to choose a single-sided test, so when we say $A \neq B$, we're really testing for $A < B$.

So how do we test whether some estimate of $(\mu_B - \mu_A)/(\sigma/\sqrt{n})$ is statistically greater than zero? We test the difference in empirical means, $M_A - M_B$, instead of $\mu_B - \mu_A$; since $M_A$ and $M_B$ are averages of Normal variables their difference is a Normal variable; and dividing the difference by the standard error makes it a random variable that is normally distributed with mean zero and variance one, a standard Normal variable, usually denoted by $Z$, hence this sometimes is called a z-test.

Observation: we're dividing the difference of the means by $\sigma/\sqrt{n}$; with $\sigma=1$, we're multiplying that difference by $\sqrt{n}$.

All we need to do now is determine whether $(M_B - M_A)\, \sqrt{25}$ is above or below the point $z^{*}$ in a standard Normal distribution where $F_{Z}(z^{*}) = .999$. That point is the critical value.

(If that scaled difference falls into the blue shaded area, then we can't reject the possibility that it was generated by randomness, instead of actual difference, with the probability that we selected; in the diagram it's 0.99, for our purposes in this post will be 0.999.)

Thanks to the miracle of computers we no longer need to look up critical values in books of statistical tables, like the peasants of old. Using, for example, the inverse normal distribution function of Apple Numbers, we learn that $F_Z(z^*) = 0.999$ implies $z^{*} = 3.09$.

So, with a sample of 25 data points, what can we conclude?

Between $A$ and $B$ our test statistic is $0.25 \times 5 = 1.25$, well below 3.09. So, $A = B$.
Between $B$ and $C$ our test statistic is $0.25 \times 5 = 1.25$, well below 3.09. So, $B = C$.
Between $A$ and $C$ our test statistic is $0.5 \times 5 = 2.5$, still below 3.09. So, $A = C$.

That was a lot of work to prove what we could see from the picture with the distributions for the average of 25 points: too much overlap and no way to tell the three variables apart from a sample of 25 points.

Ah, but we're leveling up! 100 points. Here are the distributions for the averages of 100-point samples:

Between $A$ and $B$ our test statistic is $0.25 \times 10 = 2.5$, well below 3.09. So, $A = B$.
Between $B$ and $C$ our test statistic is $0.25 \times 10 = 2.5$, well below 3.09. So, $B = C$.
Between $A$ and $C$ our test statistic is $0.5 \times 10 = 5$,  well above 3.09. So, $A \neq C$.

Now this is the weird case that gets people confused: $A = B$, $B = C$, but $A \neq C$! Equality is no longer transitive. And it does depend on $N$.

But wait, there's more. 300 more data points, and we get the 400 point case, with the following distributions:

Between $A$ and $B$ our test statistic is $0.25 \times 20 = 5$, well above 3.09. So, $A \neq B$.
Between $B$ and $C$ our test statistic is $0.25 \times 20 = 5$, well above 3.09. So, $B \neq C$.
Between $A$ and $C$ our test statistic is $0.5 \times 20 = 10$, well above 3.09. So, $A \neq C$.

Equality is again transitive. So there's only a small range of $N$ for which statistics are weird. (Not hard to figure out that range: consider it an entertaining puzzle.) This gets more complicated if those variables have different variances.

One has to be carefull with drawing inferences about equality (real equality) from statistical non-significant differences. Especially when there are small data sets and test values close to the critical values.

Yeah, about that exponential economy...

There's a lot of management and technology writing that refers to "exponential growth," but I think that most of it is a confusion between early life cycle convexity and true exponentials.

Here's a bunch of data points from what looks like exponential growth:

Looks nicely convex, and that red curve is an actual exponential fit to the data,
\[
y = 0.0057 \, \exp(0.0977 \, x)   \qquad  [R^2 = 0.971].
\]
Model explains 97.1% of variance. I mean, what more proof could one want? A board of directors filled with political apparatchiks? A book by [a ghostwriter for] a well-known management speaker? Fourteen years of negative earnings and a CEO that consumes recreational drugs during interviews?

Alas, those data points aren't proof of an exponential process, rather, they are the output of a logistic process with some minor stochastic disturbances thrown in:
\[
y = \frac{1}{1+\exp(-0.1 \, x+5)} + \epsilon_x \qquad \epsilon_x \sim \text{Normal}(0,0.005).
\]
The logistic process is a convenient way to capture growth behavior where there's a limited potential: early on, the limit isn't very important, so the growth appears to be exponential, but later on there's less and less opportunity for growth so the process converges to the potential. This can be seen by plotting the two together:

This difference is important because — and this has been a constant in the management and technology popular press — in the beginning of new industries, new segments in an industry, and new technologies, unit sales look like the data above: growth, growth, growth. So, the same people who declared the previous ten to twenty s-shaped curves "exponential economies" at their start come out of the woodwork once again to tell us how [insert technology name here] is going to revolutionize everything.

Ironically, knowledge is one of the few things that shows a rate of growth that's proportional to the size of the [knowledge] base. Which would make knowing stuff (like the difference between the convex part of an s-shaped curve and an exponential) a true exponential capability.

But that would require those who talk of "exponential economy" to understand what exponential means.

Fat tails and extremistan - not the same thing

Extremistan and mediocrestan

What, are we making up words, now? (All words are made up. Think about it.)

Extremistan and mediocrestan are characterizations of distributions; a simple way to think about them is that very large events either totally dominate (extremistan) or don't (mediocrestan):

Height is in mediocrestan: if the average height in a room with ten people is 200 cm, that's probably from ten people between 190 and 210 cm tall and not nine people 100 cm tall and one person 1100 cm tall.  

Wealth is in extremistan: if the average wealth in a room with 10 people is 2 billion dollars, that's more likely to be one billionaire with 20 billion and nine average income people than ten billionaires with 2 billion each.

This classification determines whether you can estimate relevant population parameters from samples (mediocrestan yes, extremistan no) and how well-behaved order statistics (maximum, second place, etc) are (mediocrestan nicely predictable, extremistan not so much).

There's a fairly common error that people make when they learn about extremistan: they think that because distributions in extremistan have fat tails and are dominated by extreme values, then — and this is the error — distributions that have fat tails, especially those with extreme values, are in extremistan.

Note the error: $a \Rightarrow b$ is being used to assert $b \Rightarrow a$.

As we'll see next, not all fat-tailed and extreme-valued distributions are in extremistan.

A tale of two tails

Let us compare (a) the probability that $n$ similar outcomes of large size $M$ add up to a combined event of size $nM$ (or, equivalently, average to $M$) with (b) the probability of an extreme event of size $nM$ and $n-1$ events of size 0 add up to that combined event $nM$. If the first is higher than the second, we're in mediocrestan, if the second is higher than the first, we're in extremistan.

For the Normal distribution, the probabilities (a) denoted $P(\text{Similar})$ and (b) denoted $P(\text{Extreme})$ are:
\begin{eqnarray*}
P(\text{Similar}) &=& \frac{1}{(2 \pi)^{n/2}} \exp(- n \, M^2/2) \\
P(\text{Extreme})  &=& \frac{1}{(2 \pi)^{n/2}} \exp( - n^2 \,  M^2/2)
\end{eqnarray*}
It's trivial to see that for the Normal we have
\[
P(\text{Similar}) > P(\text{Extreme}).
\]
Unsurprisingly enough, with its reference excess kurtosis of 0, the Normal distribution is well inside mediocrestan.

For our fat-tailed, extreme-valued distribution, we'll use the Gumbel distribution, which is also known as Extreme Value Type I. A simple form of this distribution has the following pdf:
\[
f_X(x) = \exp(- x - \exp(-x))
\]
As shown here, its variance is $\pi^2/6$, while the Normal above has variance 1, but since we're comparing within class (Normal with Normal and Gumbel with Gumbel), that makes no difference and saves a lot of unnecessary clutter if we just use that pdf as is.

For Gumbel we have the following probabilities:
\begin{eqnarray*}
P(\mathrm{Similar}) &=& \exp(-nM - n \, \exp(-M))
\\
P(\mathrm{Extreme}) &=&  \exp(-nM - \exp(-nM) -n+1)
\end{eqnarray*}
Since for large $M$ we have $\exp(-nM) \approx 0$ and $\exp(-M) \approx 0$, then $\exp(-nM) +n-1 > n \, \exp(-M)$, for Gumbel we also have
\[
P(\text{Similar}) > P(\text{Extreme}).
\]
The Gumbel distribution belongs in mediocrestan, despite its fat tails and extreme values.

Really makes us think about the specialness of scale independent distributions, where we can bet on a big event to overwhelm all the small events (i.e. an extremistan distribution). Those are the distributions for which a trading strategy of enduring many small losses to capture the one big win can beat a strategy of consistent small wins.

What about the maximum?

In many cases the maximum is more relevant than the mean or median. So, how do fat tails influence the maxima?

When you look at the maximum of something, say the fastest kid in a class, the larger the class, the higher the maximum will be, on average. So the fastest kid in a group of 100 is on average faster than the fastest kid in a group of 10, for example.

In mediocrestan this increase is concave on the number of kids (the difference between the fastest kids in classes of 100 and 200 kids is bigger than the difference between the fastest kids in classes of 1100 and 1200 kids, on average); in extremistan there are no guarantees.

But once again, fat tails and extreme value distributions (the Gumbel, here scaled to have variance 1) have well-behaved maxima:

This nice concavity (note the logarithmic horizontal scale) makes things predictable; since many real-world metrics are known to be fat-tailed, it's comforting to know that their maxima don't explode all of a sudden.

Note that there's an effect of the extreme value: the maxima are larger and they grow faster, with less concavity than for the Normal.

And the point is…?

There are a number of people who assert that all sorts of research and social metrics are unusable because their analysis is based on mediocrestan (either by using sample statistics to estimate population statistics or by assuming regular behavior from order statistics), but — so goes the argument — these real world metrics have fat tails, so they are in extremistan.

The point of the above was to show that this form of argument (usually punctuated with gratuitous insults, expletives, and Mathematica-based math or other forms of using pretend-math to bully one's audience) is wrong, tout court.

Only a small subset of fat-tailed, extreme-valued distributions is in extremistan. For all the rest, we can use our usual tools.

A family has two children. One is a boy. Now, do the math!

Problem

A family has two children. One is a boy. How likely is it that the other child is a boy?

Popular yet wrong solution

"There are four possible cases: two boys, a boy and a girl, a girl and a boy, and two girls. But because one child is a boy, it can't be the last case (two girls), so there are only three cases. Therefore the probability is one-third."

This solution is popular. Among others, Nassim Nicholas Taleb (on a since deleted tweet), vlogbrother Hank Green in an old SciShow episode (IIRC), probability instructors trying to show how interesting their class is to bored undergraduates, and people interviewing job candidates have used this solution.

This solution is fun because it's counter-intuitive; because of that it also looks like a smart solution.

This solution is wrong.

It's wrong because after we use "one is a boy" to eliminate the possibility of a family with two girls, we can no longer divide the probability equally among the remaining three possibilities. Equal division of probability can be used in a case of no information, but not in a case when information has already been used to change the set of possibilities.

The more attentive reader will notice that this is the same error most people make in the Monty Hall three-door problem. As a general rule, it's a bad idea to try to solve math problems by hand-waving.

If it's a math problem, do the math.*

Frequentist approach

Let's say we have a large number of cases, 4000 families for example. That's 1000 each for each combination of children: $(B,B), (B,G), (G,B)$, and $(G,G)$. Now we look at all the possibilities where we observe one of the children at random:

1000 $(B,B)$ families yield a total of 1000 boys;
1000 $(B,G)$ families yield a total of 500 boys;
1000 $(G,B)$ families yield a total of 500 boys;
1000 $(G,G)$ families yield a total of 0 boys.

We have a total of 2000 observed boys, and 1000 of these boys come from the case when the family has two boys, $(B,B)$. Half the time we observe a boy the underlying family has two boys; therefore the probability of a second boy is 1/2.

If instead of 4000 we had generic $N$ families, and called them "cases," this argument would be the frequentist derivation of the result. In frequentist parlance, the 2000 total boys are called the "possibles" and the 1000 boys from $(B,B)$ are called the "favorables." The probability is calculated as the ratio of favorables to possibles.

(The frequentist approach is how most people learn about probability and combinatorics.)

Bayesian approach

Frequentist arguments become unwieldy with more elaborate problems, so we can use this puzzle to illustrate a more elegant approach, Bayesian inference.†

First let's call things by their name: $(B,B), (B,G), (G,B)$, and $(G,G)$ are the unobserved states of the world. "One is a boy," which we'll represent by $B$, is an observed event.

Some events are uninformative, for example "one is blond," in that they don't help answer the question. Others like "one is a boy," $B$, are informative, because they help answer the question. But how can we tell?

Event $B$ is informative because it happens with different probabilities in different states of the world; therefore observing $B$ gives information about what states we're more likely to be in:

$\Pr(B|(B,B)) = 1$;
$\Pr(B|(B,G)) = 1/2$;
$\Pr(B|(G,B)) = 1/2$;
$\Pr(B|(G,G)) = 0$.

We don't know the unobserved state of the world (that is, in which of those four states the family in question falls), so in this situation we can assign equal probabilities to all four (we could look up demographics tables and confirm the numbers, but let's keep this simple):

$\Pr((B,B)) = \Pr((B,G)) = \Pr((G,B)) = \Pr((G,G)) = 1/4$.

What we want is the probability of the state $(B,B)$ having observed the event $B$; this is the conditional probability $\Pr((B,B)|B)$, which can be computed using the Bayes formula,

\[
\Pr((B,B)|B) = \frac{\Pr(B|(B,B)) \Pr((B,B))}{\Pr(B)}.
\]
Because the $\Pr(B)$ trips a lot of people, let's be clear about what it is: it's the probability that you will observe a boy in general, not in this particular case; sometimes called the a-priori probability or the unconditional probability. This is the probability that if we picked a two-child family at random and then picked one of the children at random, that child would be a boy. It's not "one, because we observe a boy," a common error.

To compute $\Pr(B)$ we must consider all four states of the world and add up ("integrate over the space of states" in expensive wording) the probability of observing a boy in each of these states weighed by the probability of the state itself:

$\begin{array}{rl}\Pr(B) =& \Pr(B|(B,B)) \Pr((B,B)) + \\
 & \Pr(B|(B,G)) \Pr((B,G)) +  \\
& \Pr(B|(G,B)) \Pr((G,B)) + \\
&\Pr(B|(G,G)) \Pr((G,G)) \\
=& 1/2
\end{array}$

(Unsurprisingly, it's 1/2, since half of the children are boys.)

Now we can compute our quantity of interest $\Pr((B,B)|B)$ by replacing the numbers in the Bayes formula. In fact, we can do that for all the states,

$\Pr((B,B)|B) = 1/2$;
$\Pr((B,G)|B) = 1/4$;
$\Pr((G,B)|B) = 1/4$;
$\Pr((G,G)|B) = 0$.

(As they used to say in the Soviet Union, trust but verify: check those numbers to be sure.)

If it's a math problem, do the math.

-- -- -- --
* "Do the math" means apply the rules of math, not just the notation and numbers.

† There's a bit of a schism in statistical modeling between frequentists and Bayesians. I'll let you figure out which side I'm on.

From value to share (in math)

There are a few places online (and in some books) where a kind of magic happens. At some point we go from something like this
\[
v(x) = w_1 x_1 + \ldots + w_K x_K
\]
to, in a single magic jump,
\[
\Pr(\text{Buy } x) = \frac{\exp(v(x))}{1 + \exp(v(x))}.
\]
The $v(x) = w_1 x_1 + \ldots + w_K x_K$ is easy enough to understand: the value of some object $x$ is a weighted sum of the dimensions of $x$; there are $K$ dimensions $x_1, \ldots, x_K$ and associated weights, $w_1, \ldots, w_K$.

But then there's something along the lines of "using appropriate assumptions, the probability of a customer purchasing $x$ is given by" that fraction.

This magic jump appears in so many places discussing analytics, data science, big data, machine learning, artificial intelligence, and other fashionable ways to refer to statistics (because that's what's there) that we could be forgiven for suspecting that the people making the jump don't know where it comes from.

It's fine not to know where something comes from, as long as one is aware of that. But it's instructive to work through the simple cases, like this one. After all, a personal trainer needn't be a super-athlete, but should at least be able to exercise.

So, how do we go from the weighted sum to those exponentials (called the logit formula)?

First, that formula is only correct if, among other things: (a) the choice is between buying $x$ or nothing at all; and (b) the $v(\cdot)$ for buying nothing is set to be zero, denote that as $v(0) =0$. Okay, the (a) makes this a simple choice and (b) is just a threshold; to move it around we can always have a constant $x_0$ added to the $v(x)$.

Now, why isn't the probability either one or zero, then? After all, either $v(x) > 0 = v(0)$ and the customer should always buy or $v(x) < 0 = v(0)$ and the customer should never buy.

That's the second step: there are things that can change the $v(x)$ that we don't observe. Maybe the customer's mood changes and the $w_i$ change with it, for example. We don't know, so we say that the customer optimizes "utility" $u(x)$ instead of value $v(x)$, and the difference is a stochastic disturbance (expensive wording for random stuff) $\epsilon$:
\[
u(x) = v(x) - \epsilon.
\]
(Let's assume that $u(0) = v(0) = 0$. It's not necessary but makes things simpler.) This second step changes the nature of the problem: instead of a decision (buy or not), we can only predict a probability:
\[
\Pr(\text{Buy } x) = \Pr(u(x)>0) = \Pr(v(x)> \epsilon).
\]
Now we can compute that probability: it's just the cumulative distribution for the $\epsilon$ variable evaluated at the point $v(x)$, denoted $F_{\epsilon}(v(x))$.

Unfortunately, by their own nature the $\epsilon$ aren't known, so we need to assume a probability distribution for $\epsilon$, one that's acceptable and is computationally convenient.

We could try the Normal distribution, which most people accept without questioning, so it gets us the first desirable characteristic. Unfortunately the $\exp(-z^2/2)$ makes integrating it difficult. What we need is a distribution that looks more or less like a Normal but can be integrated into a manageable form.

Enter the logistic distribution.

It looks enough like a Normal (it has fat tails but it's hard to tell by looking). And its probability density function
\[
f_{\epsilon}(z) = \frac{\exp(- z)}{(1+\exp(-z))^2}
\]
integrates easily to
\[
F_{\epsilon}(v) = \int_{-\infty}^{v} f_{\epsilon}(z) \, dz = \frac{1}{1 + \exp(-v)} = \frac{\exp(v)}{1+\exp(v)},
\]
which, when replaced in the formula for probability of buy gives the logit formula:
\[
\Pr(\text{Buy } x) = \Pr(v(x)> \epsilon) = F_{\epsilon}(v(x)) = \frac{\exp(v(x))}{1+\exp(v(x))}.
\]
When you know the math, you understand the magic.

Hidden factor correlation

Correlation is not causation; everyone learns to say that. But if there's a correlation, there's probably some sort of causal relationship hiding somewhere, unless it's a spurious correlation.

If two variables, $A$ and $B$ are correlated, the three simplest causal relationships are: $A$ causes $B$; $B$ causes $A$; or $A$ and $B$ are caused by an unseen factor $C$. There are many more complicated causation relationships, but these are the three basic ones.

The third case, where an unseen variable $C$ is the real source of the correlation, is what we're interested in this post. To illustrate the case let's say $C$ is a standard normal random variable, and $A$ and $B$ are noisy measures of $C$,

$ \qquad A = C + \epsilon_A$ and $ B = C + \epsilon_B$,

where the $\epsilon_i$ are drawn from a normal distribution with $\sigma_{\epsilon} = 0.05$.

To illustrate we generate 10,000 draws of $C$ and create the 10,000 $A$ and $B$ using R:

hidden_factor = rnorm(10000)
var_A_visible = hidden_factor + 0.05 * rnorm(10000)
var_B_visible = hidden_factor + 0.05 * rnorm(10000)

Now we can plot $A$ and $B$, and the correlation is obvious

And we can regress $A$ on $B$ to get the correlation and test statistics for the estimates using a linear model,

model_no_control = lm(var_A_visible~var_B_visible)
summary(model_no_control)

With the result:

Call:
lm(formula = var_A_visible ~ var_B_visible)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.271466 -0.047214 -0.000861  0.047400  0.302517 

Coefficients:
               Estimate Std. Error  t value Pr(>|t|)    
(Intercept)   0.0005294  0.0007025    0.754    0.451    
var_B_visible 0.9975142  0.0006913 1442.852  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.07025 on 9998 degrees of freedom
Multiple R-squared:  0.9952, Adjusted R-squared:  0.9952 
F-statistic: 2.082e+06 on 1 and 9998 DF,  p-value: < 2.2e-16

So, both the model and the graph confirm a strong correlation ($p < 0.0001$) between $A$ and $B$. And in many real-life cases, this is used to support the idea that either $A$ causes $B$ or $B$ causes $A$.

Now we proceed to show how the hidden factor is relevant. First, let us plot the residuals, $A-C$ against $B-C$:

The apparent correlation has now disappeared. And a linear model including the hidden factor confirms this:

model_with_control = lm(var_A_visible~var_B_visible+hidden_factor)
summary(model_with_control)

With the result

Call:
lm(formula = var_A_visible ~ var_B_visible + hidden_factor)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.18347 -0.03382 -0.00021  0.03410  0.17780 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.0007378  0.0004986   1.480    0.139    
var_B_visible 0.0004082  0.0100560   0.041    0.968    
hidden_factor 0.9997573  0.0100707  99.274  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04985 on 9997 degrees of freedom
Multiple R-squared:  0.9976, Adjusted R-squared:  0.9976 
F-statistic: 2.072e+06 on 2 and 9997 DF,  p-value: < 2.2e-16

Hidden factors are easy to test for, as seen here, but they are not always apparent. For example, in nutrition papers there's often an hidden factor relating to how health-conscious an individual is that is more often than not causing both observables (say exercising regularly and eating salads; high correlation, but exercising doesn't cause eating salads and eating salads doesn't cause exercise).

Correlation is not causation, but generally one can find a causal relationship behind a correlation, possibly one that involves hidden factors or more complex relationships.

Calculating God?

I don't believe that the God of any earthly religion is the creator of the Universe.

But I really dislike a lazy and innumerate argument commonly used to "prove" the non-existence of God, which can be summarized in the following false dichotomy:

Either there is no God and the universe just 'poofed' into existence, or there's an infinite number of Gods, because the plane of existence for each God has to be created by a higher-level God.

This is a false dichotomy: it could well be that our universe was created by a powerful being from a higher-order universe, but that universe poofed into existence without a creator. Or maybe it did have a creator, whose universe poofed into existence; or that third universe may have had a creator...

Hey, this looks like dynamic programming. I know dynamic programming.

Let's say that universes are recursively nested until one of them just poofs into existence. Of course we can't see outside our universe, but we can build simple models.

So, our universe either poofed into existence (say with probability $p$) or it was created by some higher being (with probability $1-p$). Now we iterate the process: 'level 2' universe either poofed into existence (with some probability $q$) or was created by a 'level 3' universe being (with probability $1-q$); and so on.

Time for a simplifying assumption, or as non-mathematicians call it, making things up. Let's assume that all these universes share the poofed/created probabilities, so that for any 'level $k$' universe, it poofed into existence with probability $p$ and was created by a being from a 'level $k+1$' universe with probability $1-p$.

Note that it's still possible to have an infinite number of universes, but with this formulation, the probability of a 'level $k$' universe (with us being 'level 1') being the last level is

$p (1-p)^{k-1}$.

This probability gets small pretty quickly, which suggests the 'infinite regress of universes' argument gets thin very fast.

Now we can compute the expected number of universes as a function of $p$:

$\mathbb{E}(n) = p + 2(1-p)p + 3 (1-p)^2 p + \ldots N (1-p)^{N-1} p + \ldots$

or

$\mathbb{E}(n) = p/(1-p) \times ( \text{ sum of series } N (1-p)^N )$

The sum of series $N (1-p)^N$ is $(1-p)/p^2$, so

$E(n) = 1/p$

Therefore, if we believe that the probability of a universe poofing into existence is 0.1, there are an expected ten universes; for 0.2, five universes; for 0.5, two universes.

Very far from 'turtles all the way down.'

Of course, these calculations were unnecessary, because as we know from the revelations of the prophet Terry Pratchett, it's four elephants on the back of the Great A'Tuin swimming in the Sea of Stars.