A Tetration Distribution

MSC Classification: 60E05, 26A18, 37E05

Introduction

In considering the natural operators that we typically work with in mathematical settings, we often overlook their iterative extensions—the so-called hyper-operators. They extend the original operators in an iterative way. Consider the following:

In Equation (1) a few symbols are taken for granted, particularly the binary function of addition represented by the ‘+’ symbol, and the value of the Arabic numeral ‘1’. The function $S(x)$ is known as the successor function, and it simply returns the argument of the function incremented by one.

In Equation (2) the ‘+’ symbol is defined, and we can see how addition is simply the iterative application of succession. Multiplication also extends naturally if we follow the sequence:

Next we come to exponentiation, which is not commutative, and therefore the second equality from Equations (2) and (3) does not hold. Exponentiation is right-associative:

Finally we come to the operation of tetration, upon which the distribution introduced in this paper is based:

The name tetration comes from the root prefix ‘tetra-’ and the suffix ‘-ation’ and refers to the nature of the operator as arising from iterating previous operators. Tetration is the first of the hyper-operators, a hierarchy first introduced by Goodstein [1]; the double up-arrow notation is due to Knuth [2]. The convergence theory of infinite power towers goes back to Euler [3], who showed that the iteration $y \mapsto x^y$ converges precisely for $x \in[e^{-e}, e^{1/e}]$; the constant $e^{-e}$ will play a central role in what follows. A detailed survey of iterated exponentials and infinite power towers was given by Knoebel [4]. More recent probabilistic investigations of random power towers appear in Dalthorp [5]. Ackermann’s foundational 1928 work [6] provided an early functional form that captures iterated operators. Independently, Dijoux [7] constructed a family of lifetime distributions based on continuous iteration of the exponential-minus-one function; the odd-height tetration CDF formulation of the present paper is related to but distinct from that construction.

Here we introduce a distribution that arises when we consider the map

on the unit interval. For a random variable $X$, the above map defines a valid cumulative distribution function, and its derivative therefore defines a valid probability density function:

This is the key insight of the paper. We derive a closed recursive form for the derivative in § 2 and prove that it is always positive on $(0,1)$ and integrates to one. The resulting density exhibits a characteristic shape—a steep initial rise to a peak near the origin, followed by a dip and a gradual recovery on the unit interval—which is investigated in § 3. § 4 identifies the limiting distribution as the height tends to infinity and shows that the location of the spike is governed by Euler’s constant $e^{-e}$. § 5 develops estimation and simulation for the family, reports a Monte Carlo study, fits two real datasets, and delimits the range of problems for which the family is—and is not—appropriate. § 6 notes a few qualitative morphological analogies. § 7 summarizes the results and lists open questions.

Construction and Definition

Consider the form of the cumulative distribution function. Several members of the family are plotted in Figure 1. These curves begin at zero and increase monotonically to one over their domain, so they are candidates for CDFs; their derivatives will be valid probability density functions provided they remain positive on $(0,1)$. We establish this positivity below.

Throughout we write

so that $y_m(x)=x \mathbin{\uparrow\uparrow} m$ for every positive integer $m$, whether or not $m$ is odd. The symbols $F_n$ and $f_n$ are reserved exclusively for odd $n$, where

denote the CDF and density of the height-$n$ tetration distribution. The distinction matters: when $m$ is even the function $y_m$ is a perfectly well-defined real-analytic function on $(0,1)$, but it is not a CDF, since it is not non-decreasing on the whole of $[0,1]$ (indeed $x^x\to1$ as $x\to0^+$). We return to this point in Remark 2.4. To prove validity for odd heights it is enough to verify the Kolmogorov axioms on the unit interval. We first derive the density recursion.

Proof. The relation $y_m=x^{y_{m-1}}$ implies

Differentiating both sides with respect to $x$ and multiplying through by $y_m$ immediately yields the claimed recursion. The base case $y_1=x$ gives $\frac{d y_1}{d x}=1$.

Throughout the remainder of this section we fix $x \in(0,1)$ and write $t=\ln(1/x) \in(0,\infty)$, so that the tower recursion becomes $y_1=x$ and $y_m=e^{-t y_{m-1}}$ for $m \ge2$. The following elementary bounds on the towers are the key to the positivity proof.

Case 2: $0<g_{n-2}<1/t$. Then $0<1-tg_{n-2}<1$, so $0<g_{n-1}<y_{n-2}$. Since $n-2\ge3$ is odd, Lemma 2.2 gives $ty_{n-2}<1$, whence $tg_{n-1}<ty_{n-2}<1$.

In either case $tg_{n-1}<1$, and therefore $g_n=y_{n-1}(1-tg_{n-1})>0$. By induction $g_n>0$, and hence $y'_n>0$, for every odd $n$.

Remark 2.4 (Status of the even-index towers) The recursion of Theorem 2.1, and equivalently the density recursion (11) below, evaluates the height-$n$ density for odd $n$ in terms of the quantities $y_{n-1}$ and $y'_{n-1}$ of even index. It is worth stating explicitly that this involves no circularity and no appeal to any object of doubtful standing.

For even $m$, the function $y_m$ defined by (8) is a composition of finitely many real-analytic maps and is therefore itself real-analytic, finite and strictly positive on $(0,1)$; likewise $y'_m$ is real-analytic and finite there. What fails for even $m$ is only the monotonicity required of a CDF, and consequently $y'_m$ is not a probability density: by (10), $g_2=x(1-t)$ is negative exactly when $t>1$, that is on $x\in(0,e^{-1})$, so $y'_2<0$ there. The even-index towers are used here purely as algebraic intermediates in a recursion among analytic functions, and carry no probabilistic interpretation whatsoever. This is why the notation of (9) reserves $F$ and $f$ for odd indices.

Crucially, the proof of Theorem 2.3 is arranged so that no positivity is ever required of an even-index quantity. Case 1 of the inductive step explicitly admits $g_{n-1}\le0$, and the conclusion $g_n=y_{n-1}(1-tg_{n-1})>0$ is reached in that case precisely because $g_{n-1}$ is non-positive, since then $1-tg_{n-1}\ge1>0$ and $y_{n-1}>0$. The sign changes of the even-index derivatives are thus harmless: the induction passes through them rather than around them, and the positivity of the odd-index densities is established independently of the sign of any intermediate term. The same observation applies verbatim to the numerical evaluation of (11), which involves only finite analytic quantities at every step.

Properties of the Density

Probability density functions for several of these distributions are shown in Figure 2. The density $f_n$ is obtained from the recursive formula derived in § 2:

where $F_n(x)=x \uparrow\uparrow n$ denotes the CDF and $y_{n-1}$ is the even-index auxiliary tower of (8); see Remark 2.4. From the graphs one observes a characteristic shape: an initial spike near zero, followed by a dip below height one, and then a gradual rise. The height of the power tower primarily affects the sharpness of the initial spike and how quickly the density falls below one.

Mode

Numerical maximization of the density yields the approximate modes reported in Table 1.

$n$	3	5	7	9	11	15	31	101	501	2001
Mode	0.06399	0.06344	0.06316	0.06303	0.06298	0.06300	0.06346	0.06464	0.06560	0.06588
$10^3\delta_n$	2.00	2.55	2.83	2.96	3.01	2.98	2.53	1.35	0.38	0.11

Although these values are obtained numerically and therefore carry a small discretization error, every one of them lies just below $e^{-e}\approx0.065988$. The mode decreases until roughly height 11–15 and then rises again, approaching $e^{-e}$ from the left rather than drifting toward zero as the height increases: the deficit $e^{-e}-\text{mode}$ peaks near $n=11$ and thereafter falls by more than an order of magnitude between $n=101$ and $n=2001$. This is not a numerical coincidence. § 4 shows that $e^{-e}$ is exactly Euler’s convergence boundary for infinite power towers and that the limiting density has a one-sided singularity there, which forces the modes to converge to $e^{-e}$ from below.

Mean and Variance

Approximate means and variances, obtained by numerical integration, are reported in Table 2.

$n$	3	5	7	9	11	15	31	101	501	2001	$\infty$
Mean	0.4269	0.4024	0.3912	0.3850	0.3812	0.3769	0.3717	0.3692	0.3687	0.3687	0.3687
Variance	0.0911	0.0970	0.1007	0.1031	0.1048	0.1068	0.1096	0.1110	0.1113	0.1113	0.1113

As the height increases the mean declines slowly while the variance grows, both quantities converging to the corresponding moments of the limit distribution identified in § 4.

At first sight the behaviour of the variance is counter-intuitive: the mode sharpens dramatically as $n$ grows, and a sharper mode is normally associated with less dispersion, yet the variance increases from 0.0911 to 0.1113. The resolution is that the sharpening peak is not drawing mass in from the whole interval; it is separating the distribution into two groups at very different scales. As Theorem 4.3 below shows, in the limit a proportion exactly $1/e \approx0.368$ of the total mass is compressed into the short interval $[0,e^{-e}]$, a region of length under 0.066, while the remaining $1-1/e \approx0.632$ is spread across $[e^{-e},1]$. It is this growing separation between a tightly concentrated low-value group and a broadly dispersed high-value group—rather than any spreading of the peak itself—that drives the variance upward. The naive intuition “sharper mode implies smaller variance” is reliable only for unimodal densities whose mass remains concentrated about a single scale, which is precisely what fails here.

Endpoint behaviour

We now examine the density $f_n(x)$ near the endpoints of the unit interval.

As x0+.

Although the power tower itself vanishes at the origin for every odd height $n \ge3$, the density does not: evaluation of the recursion (11) shows that $f_n(x) \to1$ as $x \to0^+$ for every odd $n$. From this common boundary value the density rises extremely rapidly to its mode just below $e^{-e}$ and then falls. The taller the tower, the more violent this initial spike: the peak becomes both taller and narrower as the height increases.

As $x \to1^{-}$

At the right endpoint the situation is far more regular. Direct evaluation of the recursive formula shows

for every height. Moreover the density approaches this value smoothly from below after the mid-interval dip. There is no singularity or rapid oscillation near 1; the behaviour is essentially the same for all odd heights once $x$ is bounded away from zero.

The Limit Distribution and Euler's Constant $e^{-e}$

We now identify the limiting behaviour of the family as the height tends to infinity through odd integers. This resolves the shape of the spike observed in § 3.1 and explains the role of the constant $e^{-e}$.

Euler's convergence boundary

Fix $x \in(0,1)$ and set $t = \ln(1/x) > 0$. The tower recursion (8) is the orbit of $y_1 = x$ under the map

so that $y_m = \varphi_x^{m-1}(y_1)$. The map $\varphi_x$ is a decreasing contraction-or-expansion whose unique fixed point $y^*(x)$ solves $y = x^y$, equivalently $ty e^{ty} = t$, and is therefore given in closed form by the Lambert $W$ function [8] as

Since $\varphi_x'(y) = -t e^{-ty}$, the multiplier at the fixed point is $\varphi_x'(y^*) = -t y^* = -W(t)$. The fixed point is therefore attracting precisely when $W(t) < 1$, that is when $t < e$, that is when $x > e^{-e}$; it loses stability exactly at

using $W(e) = 1$. This is Euler’s classical convergence boundary for infinite power towers [3, 4]: the multiplier passes through $-1$, which is a period-doubling bifurcation. For $x < e^{-e}$ the fixed point is repelling and the orbit is attracted instead to a two-cycle $\{y_-(x), y_+(x)\}$ with $y_-(x) < y^*(x) < y_+(x)$, the two points being the solutions of $y = \varphi_x(\varphi_x(y))$ other than $y^*$.

Because $y_1 = x$ is small for small $x$ while $y_2 = x^x$ is close to 1, the odd-index iterates track the lower point of this two-cycle and the even-index iterates track the upper one. This dichotomy is the exact analytic counterpart of the even/odd distinction imposed in § 2.

The following observation is what makes the analysis elementary, and in particular rules out any behaviour more complicated than a two-cycle.

Proof. By the chain rule,

the product of two negative factors. A strictly increasing continuous self-map of an interval has only monotone orbits: if $\psi_x(y)>y$ then applying the increasing map $\psi_x$ repeatedly preserves the inequality, so the orbit is increasing, and similarly in the decreasing case. Being monotone and bounded ($\varphi_x$ maps into $(0,1)$), each orbit converges, necessarily to a fixed point of $\psi_x$ by continuity. This excludes cycles of period greater than two for $\varphi_x$, and in particular excludes any further period-doubling cascade.

For the final claim, $y_3>y_1$ reads $x^{x^x}>x$. Taking logarithms gives $x^x\ln x>\ln x$, and dividing by $\ln x<0$ reverses the inequality to $x^x<1$, which holds for every $x\in(0,1)$. □

Proof. Pointwise convergence. Fix $x\in(0,1)$. By Lemma 4.1 the odd-index subsequence $y_1,y_3,y_5,\ldots$ is a strictly increasing orbit of $\psi_x$ and therefore converges to the smallest fixed point of $\psi_x$ lying at or above $y_1=x$.

We first record that

Indeed, writing $w=W(t)>0$ so that $t=we^w$, the inequality $x=e^{-t}<W(t)/t=y^*$ is equivalent to $te^{-t}<w$, that is to $e^{w-e^w}<1$, that is to $w(1-e^w)<0$, which holds because $w>0$ and $e^w>1$.

By Euler’s classical convergence theorem [3, 4], the full tower sequence $(y_m)_{m\ge1}$ converges exactly when $x\ge e^{-e}$, and for $x<e^{-e}$ it splits into two subsequential limits forming a two-cycle of $\varphi_x$. We use this only to know whether a two-cycle exists; the identification of the branch is then forced by monotonicity, as follows.

If $x\ge e^{-e}$, the whole sequence converges, necessarily to the unique fixed point, so the odd subsequence has limit $y^*(x)=W(t)/t$. (At $x=e^{-e}$ the fixed point is neutral and the convergence is algebraic rather than geometric, but convergence still holds.)

If $x<e^{-e}$, then $\varphi_x$ has the three fixed points $y_-<y^*<y_+$. Since $\psi_x$ is increasing, the sign of $\psi_x(y)-y$ is constant between consecutive fixed points and alternates across them; as $\varphi_x(0^+)=e^{-t}>0$ we get $\psi_x(y)>y$ on $(0,y_-)$ and $\psi_x(y)<y$ on $(y_-,y^*)$. Orbits starting in $(y_-,y^*)$ are therefore decreasing. But by Lemma 4.1 our orbit is strictly increasing, and by (15) it starts at $x<y^*$; hence $x\in(0,y_-)$ and the orbit increases to $y_-(x)$.

$F_{\infty}$ is a continuous CDF. Each $F_n$ is non-decreasing by Theorem 2.3, so the pointwise limit $F_{\infty}$ is non-decreasing. As $x \to0^{+}$ we have $t \to\infty$ and $0 < F_{\infty}(x) \le y^{*}(x) = W(t)/t \to0$, so $F_{\infty}(0+) = 0$; and $W(t)/t \to1$ as $t \to0^{+}$, so $F_{\infty}(1-) = 1$. On $(e^{-e}, 1)$ the function $W(t)/t$ is real-analytic, and on $(0, e^{-e})$ the branch $y_{-}$ is real-analytic by the implicit function theorem, the relevant derivative $\partial_y(\psi_x(y)-y)$ being non-zero at $y_{-}$ away from the bifurcation point. Continuity at $x = e^{-e}$ follows from the expansion of Proposition 4.4 below, in which $y_{-}(x) \to1/e = y^{*}(e^{-e})$ as $x \uparrow e^{-e}$. Hence $F_{\infty}$ is a continuous CDF on $[0,1]$.

Uniformity. Pointwise convergence of distribution functions to a continuous limit is exactly convergence in distribution, and by Pólya’s theorem such convergence is automatically uniform on $[0,1]$. □

Figure 3 displays the convergence. By height $n = 2001$ the tower is numerically indistinguishable from $F_{\infty}$ at the resolution of the plot. Away from the critical point the convergence is geometric and very fast: in double precision $|F_{1001} - F_{\infty}|$ is already at the level of rounding error (below $10^{-15}$) for $x \le0.06$ and for $x \ge0.1$. Inside that window the rate degrades as $e^{-e}$ is approached, the multiplier $|-W(t)|$ tending to $1$; at $x = e^{-e}$ itself convergence is algebraic, and $|F_{1001} - F_{\infty}| \approx2.8 \times10^{-2}$ there. This slow approach precisely at the critical point is the finite-height signature of the singularity described in Proposition 4.4.

The critical point

Proof. At $x=e^{-e}$ we have $t=e$, and $W(e)=1$ since $1\cdot e^1=e$. Substituting into (13) gives $y^*=1/e$, and by Theorem 4.2 this value is $F_\infty(e^{-e})$. $\square$

Theorem 4.3 is the exact statement underlying the appearance of $1/e$ in the numerical results of § 3.2. It is worth emphasising that $1/e$ arises here as a probability, not as the mean; see Remark 4.6.

Proof sketch. At the critical parameter the multiplier of $\varphi_x$ at $y^*$ equals $-1$, so the bifurcation is a generic period doubling: writing $u=y-y^*$ and $t=e+\delta$, the second-iterate map $\varphi_x\circ\varphi_x$ has the normal form $u\mapsto u+a\delta u+b u^3+O(\delta u^2,u^4)$ with $a,b\ne0$, whence the two-cycle points separate from $y^*$ at rate $|u|\asymp\sqrt{\delta}$. Since $t=\ln(1/x)$ is analytic with $dt/dx=-1/x=-e^e\ne0$ at $x=e^{-e}$, we have $\delta\asymp h$ and therefore $y^*-y_-\asymp\sqrt{h}$. Differentiating the expansion gives the stated behaviour of $f_\infty$. Analyticity away from the critical point is the implicit function theorem, as in Theorem 4.2. $\square$

The constant $C$ was obtained numerically: computed in 60-digit arithmetic, the ratio $(1/e-F_\infty(e^{-e}-h))/\sqrt{h}$ increases monotonically through the values 2.127450, 2.127591, 2.127635, 2.127649 and 2.127653 at $h=10^{-8},10^{-9},\ldots,10^{-12}$. The $O(h)$ term in the expansion predicts that this ratio should approach $C$ from below at rate $O(\sqrt{h})$, that is, by a factor of $\sqrt{10}\approx3.16228$ per decade in $h$. The observed successive increments are $1.4026\times10^{-4}$, $4.4349\times10^{-5}$, $1.4024\times10^{-5}$ and $4.4347\times10^{-6}$, shrinking by factors 3.16272, 3.16242 and 3.16232—converging to $\sqrt{10}$, which both confirms the square-root branching of Proposition 4.4 and justifies Richardson extrapolation of the sequence. That extrapolation is itself stable, returning 2.127655467, 2.127655458, 2.127655457 and 2.127655457 from the four consecutive pairs, giving $C=2.1276555\ldots$. We have not identified a closed form for $C$; extracting one from the period-doubling normal form of Proposition 4.4 is a natural next step.

Proposition 4.4 accounts for the entire qualitative picture of § 3.1. The limiting density is unbounded as $x\uparrow e^{-e}$ and bounded on the other side, so the maximum of $f_n$ must occur to the left of $e^{-e}$ for every finite $n$ and must migrate toward $e^{-e}$ as $n$ grows. This is exactly the pattern of Table 1, in which every mode lies below $e^{-e}$ and the deficit decreases by more than an order of magnitude from $n=101$ to $n=2001$. The finite-height densities are smooth approximations to a function with an integrable singularity: the peak grows without bound while the area it encloses tends to zero, which is why the limit law has no atom at $e^{-e}$ despite the diverging spike.

Limiting moments

Proof. The integral representations are the standard tail formulas for a non-negative random variable bounded by $1$, valid since $F_{\infty}$ is a continuous CDF supported on $[0,1]$ by Theorem 4.2. The numerical values were obtained by splitting the range at $e^{-e}$ and applying the substitution $x=e^{-e}-u^{2}$ on the lower piece, which removes the square-root singularity of Proposition 4.4 and renders both integrands smooth; the quoted digits are stable under refinement and agree with direct quadrature against the height-6401 tower.

Statistical Inference

The preceding sections treat the family analytically. We now turn to inference: estimation of the height from data, a simulation study of the resulting estimators, and an application to two real datasets bounded on the unit interval. The outcome of the application is negative—the family is decisively outperformed by the Beta and Kumaraswamy distributions—and Section 5.5 explains structurally why, which we regard as the more useful contribution of this section.

Maximum likelihood

Given a sample $x_1,\ldots,x_m$ from the unit interval, the log-likelihood of the height-$n$ tetration distribution is

with $f_n$ given by the recursion (11). The parameter is a discrete one—an odd positive integer—so maximum likelihood is an exact search over a grid of odd heights rather than a continuous optimisation, and no regularity conditions or derivative computations are required.

A practical remark makes this cheap. Evaluating (11) separately for each candidate height costs $O(N^2)$ passes over the data when heights up to $N$ are considered, but the recursion builds every height on the way to the tallest, so recording $(F_k,f_k)$ at each odd $k$ during a single sweep delivers the entire profile in $O(N)$ passes. The whole profile over $n\le101$ is therefore obtained at the cost of evaluating the tallest tower once.

A two-parameter extension.

The one-parameter family is very rigid (§ 5.5), so we also consider the exponentiated version, in the sense of Gupta and Kundu [9],

which is a valid CDF for every odd $n$ and every $\alpha>0$, and reduces to the $\operatorname{Beta}(\alpha,1)$ law at $n=1$. Its density is $f_{n,\alpha}=\alpha F_n^{\alpha-1}f_n$, so

For fixed $n$ the score equation in $\alpha$ solves in closed form,

the denominator being strictly negative because $0<F_n<1$ on $(0,1)$. The two-parameter MLE is thus again an exact profile search over odd $n$, with $\alpha$ eliminated analytically at each step; no numerical optimiser is needed. Since $\partial^2\ell/\partial\alpha^2=-m/\alpha^2$, the observed information gives the conditional standard error $\operatorname{se}(\hat{\alpha})=\hat{\alpha}/\sqrt{m}$, valid given $\hat{n}$.

Method of moments.

Matching the sample mean to $\mathbb{E}[X_n]$ is possible in principle but of very limited use, because the attainable means form a short interval: by Table 2 the mean decreases monotonically from $0.5$ at $n=1$ to the limiting value $0.3686906\ldots$ of Proposition 4.5, so

for every odd $n.

A moment estimator therefore exists only for samples whose mean lies in that interval, which excludes both datasets analysed below. This is the first symptom of the rigidity discussed in § 5.5. The exponentiated family (17) does not suffer from it: its mean ranges over all of $(0,1)$ as $\alpha$ varies.

Simulation

Because $F_n$ is continuous and strictly increasing on $(0,1)$ by Theorem 2.3, inverse-transform sampling applies without qualification: if $U \sim\operatorname{Unif}(0,1)$ then $F_n^{-1}(U)$ has the height-$n$ tetration distribution, and $F_n^{-1}$ is evaluated by bisection, which is unconditionally convergent here. Bisecting in $\log x$ rather than $x$ is advisable, since for tall towers the lower quantiles are extremely small. For the exponentiated family one simply draws $F_n^{-1}(U^{1/\alpha})$.

Monte Carlo study

Table 3 reports the behaviour of the maximum likelihood estimator of the height, over 2000 replicates at each combination of true height and sample size, with the search restricted to $n \in\{1,3,\ldots,41\}$. Bias and mean squared error both decrease steadily with the sample size and the probability of exact recovery increases, as consistency requires. Taller towers are harder to identify, which is expected: the CDFs $F_n$ converge to $F_\infty$, so neighbouring tall heights are increasingly similar. Even so the profile likelihood is sharply peaked rather than flat—for a sample of 5000 drawn at $n=11$, the log-likelihood at $n=15$ is already 24 units below the maximum—so the height is genuinely identified, not merely bounded.

$n$	$m$	mean($\hat{n}$)	bias	MSE	$\Pr(\hat{n}=n)$
3	20	4.993	1.993	47.562	0.391
3	50	3.740	0.740	11.012	0.568
3	100	3.332	0.332	1.452	0.718
3	500	3.056	0.056	0.116	0.971
5	20	7.528	2.528	75.384	0.254
5	50	6.009	1.009	18.410	0.351
5	100	5.356	0.356	4.136	0.499
5	500	5.160	0.160	0.620	0.848
11	20	13.730	2.730	116.668	0.115
11	50	12.574	1.574	48.432	0.193
11	100	11.738	0.738	20.464	0.222
11	500	11.236	0.236	2.748	0.509

Table 4 reports the corresponding study for the exponentiated family, focusing on $\hat{\alpha}$ from (18). The bias is small already at $m=50$ and the mean squared error falls at the expected rate.

$n$	$\alpha$	$m$	mean($\hat{\alpha}$)	bias	MSE	$\Pr(\hat{n}=n)$
3	1.5	20	1.6287	0.1287	0.2390	0.323
3	1.5	50	1.5477	0.0477	0.1060	0.446
3	1.5	100	1.5105	0.0105	0.0607	0.541
3	1.5	500	1.5131	0.0131	0.0114	0.871
5	0.8	20	0.8543	0.0543	0.0575	0.208
5	0.8	50	0.8180	0.0180	0.0206	0.298
5	0.8	100	0.8126	0.0126	0.0097	0.429
5	0.8	500	0.8037	0.0037	0.0018	0.806

Application to data

We fit the family to two real datasets supported on the unit interval, both distributed with the datasets package of R [10]. The first is the pore-shape index of 48 rock samples from a petroleum reservoir (rock$shape; image analysis by R. Katz, University of Oxford, from BP Research data), a dimensionless measure of pore cross-sectional shape recorded in $(0,1)$. ¹ The second is the proportion of army draftees with education beyond primary school in each of 47 French-speaking Swiss provinces around 1888 (swiss$Education, divided by 100), from Mosteller and Tukey [11]. Both were chosen because they are right-skewed with mass concentrated at small values, which is the regime in which the tetration density is least unlike real data; no dataset on which the family performs better was found.

We compare the one-parameter tetration family and its exponentiated extension against the Beta distribution and the Kumaraswamy distribution [12], the two standard two-parameter competitors on the unit interval. Table 5 reports maximised log-likelihoods, AIC, BIC, and the Kolmogorov–Smirnov statistic with its $p$-value; Figure 4 shows the fitted densities.

Data	Model	log $L$	$k$	AIC	BIC	KS	$p$
rock$shape<br>$m=48$ | Beta | 55.603 | 2 | −107.205 | −103.463 | 0.1428 | 0.256
rock$shape<br>$m=48$ | Kumaraswamy | 52.494 | 2 | −100.988 | −97.246 | 0.1534 | 0.188
rock$shape<br>$m=48$ | Tetration | 12.891 | 1 | −23.783 | −21.912 | 0.4130 | $<10^{-4}$
rock$shape<br>$m=48$ | Exp-tetration | 13.174 | 2 | −22.347 | −18.605 | 0.4099 | $<10^{-4}$
swiss$Education<br>$m=47$ | Beta | 59.607 | 2 | −115.214 | −111.514 | 0.1667 | 0.131
swiss$Education<br>$m=47$ | Kumaraswamy | 58.672 | 2 | −113.345 | −109.644 | 0.1734 | 0.105
swiss$Education<br>$m=47$ | Exp-tetration | 45.780 | 2 | −87.559 | −83.859 | 0.3006 | 0.0003
swiss$Education<br>$m=47$ | Tetration | 40.181 | 1 | −78.363 | −76.513 | 0.4385 | $<10^{-4}$

The conclusion is unambiguous and negative. On both datasets the Beta and Kumaraswamy fits are superior by 40–80 units of AIC, and the Kolmogorov–Smirnov test rejects both tetration models at any conventional level while rejecting neither competitor. The fitted heights are $\hat{n}=5$ and $\hat{n}=19$; the exponentiated extension recovers only a small part of the gap. We also fitted the remaining unit-interval datasets in the same package, with the same outcome; on four of them the maximised likelihood occurs at $n=1$, which is the uniform distribution, so that the family has nothing to offer beyond the trivial member. The one-parameter tetration law is also beaten by one-parameter competitors such as the power law $F(x)=x^{a}$ and the truncated exponential, so the deficit is not an artefact of comparing one parameter against two.

Scope of applicability

The failure is systematic rather than accidental, and its cause is visible in § 3: the endpoint values of the density are the same for every member of the family,

and the interior mode is confined to a neighbourhood of $e^{-e}$, converging to it as the height grows (Proposition 4.4). The height $n$ controls only the sharpness of the spike. The family therefore has no freedom to adjust either tail: it cannot produce a density vanishing at an endpoint, nor one diverging there, whereas the Beta, Kumaraswamy and power families all do so freely through their shape parameters. The right-hand panel of Figure 4 shows the consequence directly—the fitted tetration density is still at height $\approx0.6$ at $x=0.5$, where the data have essentially no mass, purely because it must return to 1 at $x=1$.

Two honest conclusions follow. First, the tetration family should not be recommended as a general-purpose model for bounded data; on the evidence here the Beta and Kumaraswamy families dominate it comprehensively. Second, the family is plausibly useful only in the narrow regime its geometry actually describes: data that are close to uniform over most of the interval but carry a pronounced spike at small values near $e^{-e}$, with no systematic decay toward either endpoint. We have not identified a dataset of that kind. The contribution of the present paper is accordingly the analytic structure of § 2–§ 4—in particular the limit law and the role of Euler’s constant—rather than a competitive modelling tool, and the estimation theory above should be read as establishing what the family can and cannot do.

Qualitative Analogies

There are a few natural phenomena that exhibit a qualitatively similar trend: an initial response that overshoots a baseline level, drops rapidly, and then gradually recovers. Two examples that first suggested themselves when examining the form of the density are the classical illustration of the Dunning–Kruger effect and the technology hype cycle.

In the usual schematic of the Dunning–Kruger effect the horizontal axis represents an individual’s actual knowledge of a subject while the vertical axis represents expressed confidence. The resulting curve shows a high initial peak at low knowledge, followed by a trough and a later rise. The tetration density produces a qualitatively similar profile. Although more recent work has questioned the empirical accuracy of that particular curve, the originally postulated shape remains a useful morphological analogy for the density studied here.

Likewise, the technology hype cycle traces a rapid rise in attention or notoriety, a subsequent trough of disillusionment, and a gradual recovery. Time series of public interest in novel technologies (for instance as measured by Google Ngram frequencies) often display this same general shape. The tetration density offers a simple one-parameter family that reproduces the characteristic peak–dip–rise pattern.

Two caveats should be stated plainly. First, these analogies are morphological only; no claim is made that either phenomenon is generated by the mechanism studied here. Second, and more concretely, the curves just described are drawn on arbitrary axes—time, competence, confidence—whereas the tetration distribution lives strictly on $(0,1)$. Any use of the family on such a domain therefore requires an affine transformation

which has CDF $F_Y(y)=F_n((y-a)/(b-a))$ and density $f_Y(y)=(b-a)^{-1}f_n((y-a)/(b-a))$ on $[a,b]$. Under this transformation the mode sits at $a+(b-a)\operatorname{mode}(f_n)$, so the characteristic peak occurs at roughly 6.6% of the way along the range, and the limiting mass below it is $1/e$ by Theorem 4.3. The endpoints $a$ and $b$ are genuine additional parameters and would have to be estimated alongside the height $n$ in any applied use of the family. It should be said plainly that the affine map rescales the interval but does not alter the shape of the density, so it does not relieve the rigidity documented in § 5.5; these analogies remain illustrations of a shape, not evidence that the family models the phenomena.

Conclusions

We have introduced a new one-parameter family of distributions on the unit interval generated by odd-height tetration. The cumulative distribution function is given explicitly by the power tower of height $n$, and the density satisfies a simple recursive relation. Numerical investigation reveals a characteristic shape consisting of a sharp spike near the origin, a subsequent dip, and a gradual rise toward the right endpoint.

The main structural result is the identification of the limit as the height tends to infinity through odd integers. The limiting CDF is $W(t)/t$ with $t = \ln(1/x)$ above Euler’s constant $e^{-e}$ and the lower branch of the attracting two-cycle of $y \mapsto x^y$ below it, and the convergence is uniform. The constant $e^{-e}$ is precisely the parameter value at which $y \mapsto x^y$ undergoes a period-doubling bifurcation, which produces a one-sided inverse-square-root singularity in the limiting density; this explains the spike observed at every finite height and shows that the modes converge to $e^{-e}$ from the left. At the critical point the limit law satisfies $F_\infty(e^{-e}) = 1/e$ exactly.

The limiting mean and variance are $0.3686906\ldots$ and $0.1112922\ldots$. The mean is close to $1/e$ but is provably different from it; the constant $1/e$ enters this problem as a probability rather than as a mean. No closed form has been identified for either limiting moment, nor for the branching constant $C$ of Proposition 4.4, and closed-form expressions for the moment-generating function and for moments of order greater than two remain unknown.

On the inferential side, § 5 shows that estimation itself poses no difficulty. The height is discrete, so maximum likelihood is an exact profile search, computable for all heights up to $N$ in $O(N)$ passes over the data; in the exponentiated extension the second parameter is removed in closed form; and a Monte Carlo study confirms that bias and mean squared error decay as consistency requires. What fails is the fit. On two real unit-interval datasets the Beta and Kumaraswamy distributions outperform the family by 40–80 units of AIC, and a Kolmogorov–Smirnov test rejects the tetration models while rejecting neither competitor. The reason is structural: $f_n(0^+) = f_n(1) = 1$ for every height, so no member of the family can adapt either tail, and the height governs only the sharpness of the spike near $e^{-e}$. We therefore do not recommend the family as a general-purpose model for bounded data, and we have stated in § 5.5 the narrow geometry for which it could be suitable.

Several analytic directions remain open. Closed forms for the two limiting moments and for the branching constant $C$ of Proposition 4.4, and the rate at which the finite-height modes approach $e^{-e}$, are natural questions suggested by the period-doubling picture. A further possibility, suggested by the rigidity just described, is to enlarge the family by composing it with a flexible transformation of the unit interval—the affine map (20) being the simplest—so as to free the endpoint behaviour while retaining the tetration structure; whether anything of the characteristic geometry survives such an enlargement is not clear to us.

The distinctive peak–dip–rise profile of the density bears a qualitative resemblance to certain temporal or competence-related curves that exhibit an early overshoot followed by recovery. These morphological parallels are noted only as illustrations of the shape; formal statistical modelling of such phenomena is left for future work.

Declarations

Funding

The author declares that no funds, grants, or other support were received during the preparation of this manuscript.

Conflict of interest

The author has no relevant financial or non-financial interests to disclose.

Data availability

The two datasets analysed in § 5.4 are publicly available and require no special access: both are distributed with the datasets package of R [10] as rock$shape and swiss$Education, the latter originating with Mosteller and Tukey [11]. All other numerical results are reproducible directly from the formulas given herein; the scripts used to generate the figures, tables and simulation study are available from the author upon reasonable request.