Jaynes-Gibbs Entropic Convex Duals and Orthogonal Polynomials

Abstract

The univariate noncentral distributions can be derived by multiplying their central distributions with translation factors. When constructed in terms of translated uniform distributions on unit radius hyperspheres, these translation factors become generating functions for classical families of orthogonal polynomials. The ultraspherical noncentral t, normal N, F, and ${\chi }^2$ distributions are thus found to be associated with the Gegenbauer, Hermite, Jacobi, and Laguerre polynomial families, respectively, with the corresponding central distributions standing for the polynomial family-defining weights. Obtained through an unconstrained minimization of the Gibbs potential, Jaynes’ maximal entropy priors are formally expressed in terms of the empirical densities’ entropic convex duals. Expanding these duals on orthogonal polynomial bases allows for the expedient determination of the Jaynes–Gibbs priors. Invoking the moment problem and the duality principle, modelization can be reduced to the direct determination of the prior moments in parametric space in terms of the Bayes factor’s orthogonal polynomial expansion coefficients in random variable space. Genomics and geophysics examples are provided.

1. Introduction

We shall argue that the four noncentral univariates t, F, the normal N, and ${\chi }^2$ distributions $\rho \left(r\right|r_o),$ with r being the one-dimensional random space variable, and $r_o$ the one-dimensional noncentrality parameter of the respective distributions, can all be constructed in a modular fashion by multiplying their central counterparts $\rho \left(r\right|r_o=0)$ with a factor $T\left(r\right|r_o)$ effecting a central distribution translation, that is,

$$\rho \left(r\right|r_o)=T\left(r\right|r_o)\rho \left(r\right|0).$$

In statistical parlance, the noncentral distributions are needed to estimate or modelize effect sizes [1]. With the exception of the normal distribution, these translations are non-shape-preserving. The derivation of the translation factors $T\left(r\right|r_o)$ can be carried out in two manners, depending on whether primacy is put upon translated uniform density distributions on unit radius hyperspheres, as is done in this manuscript, or on translated normal distributions, as done classically. We shall review both derivations herein, with emphasis on the hyperspherical distributions.

We choose to place primacy upon the simple uniform density distribution on the unit radius hypersphere ${\mathcal{S}}^{\nu }$, where $\nu$ stands for either the dimension of the hypersphere surface ${\mathcal{S}}^{\nu }\in {\Re }^{\nu +1}$, or, equivalently, the degrees of freedom (dof) of the specified distribution. It is known that the projection of a unit radius uniform density hyperspherical distribution on ${\mathcal{S}}^{\nu }$ on any given polar axis readily provides us with the central t distribution [2], and that such a projection converges with the central normal distribution $N(0,1)$ with null mean when $\nu$ tends to infinity [3] (p. 59). This observation provides us with the needed building principle used throughout this manuscript: use the uniform density distribution on the unit radius hypersphere ${\mathcal{S}}^{\nu }\in {\Re }^{\nu +1}$ to derive modular expressions for the central and noncentral t, F, N, and ${\chi }^2$ distributions. In order to proceed, one needs to master some very simple notions concerning the hypersphere geometry. The projection of a random unit vector $\mathbf{x}$ on the unit radius hypersphere ${\mathcal{S}}^{{\nu }_2}$ on any given unit polar axis $\mathbf{p}$ naturally defines a polar angle $\theta$ through the scalar product $cos\theta =\mathbf{x}·\mathbf{p}$. As such, the central t distribution with ${\nu }_2$ degrees of freedom can be drawn on the compact support $-1<cos\theta <1$, on which it acquires a simple expression in trigonometric terms: it is simply proportional to ${sin}^{{\nu }_2-1}\theta$ after integration of the azimuthal coordinates. See Saville and Wood [2] for a very extensive digression on the subject. On the familiar sphere ${\mathcal{S}}^2\in {\Re }^3$, the latter provides us with the well-known spherical surface element $sin\theta \mathrm{d}\theta$ after the integration of the azimuthal coordinate. Similarly, the central $F^{({\nu }_1,{\nu }_2)}$ distribution becomes proportional to ${cos}^{{\nu }_1-1}\theta {sin}^{{\nu }_2-1}\theta$ on the compact domain $0\le cos\theta <1$, where $\theta$ is the angle between a random vector $\mathbf{x}$ on the hypersphere ${\mathcal{S}}^{{\nu }_1+{\nu }_2}$ and a secant hyperplane defining the subspace ${\mathcal{S}}^{{\nu }_1}$. In statistical parlance, the ${\mathcal{S}}^{{\nu }_1}$ and ${\mathcal{S}}^{{\nu }_2}$ subspaces refer to the between-class and within-class variance spaces in an analysis variance (ANOVA), as reflected by the F-statistic redefinition in trigonometric terms, found in Section 2. As it turns out, these two central uniform hyperspherical distributions are all that is needed to proceed with the derivation of all the univariate noncentral t, F, N, and ${\chi }^2$ distributions where, for uniformity of designation, a normal distribution $N(\delta ,1)$ with non-vanishing mean $\delta$ and unit variance will be simply referred to as the noncentral normal distribution. Finally, in the theory of orthogonal polynomials, the designation ultraspherical polynomials—also known as Gegenbauer polynomials [4]—has prevailed over that of hyperspherical polynomials, and we shall abide by this nomenclature. In order to distinguish the ultraspherical t and F distributions from the ones derived classically from the normal N and ${\chi }^2$ distributions, we shall designate the former densities by the Greek letter $\upsilon$ (upsilon), as in υπερσφαίρα (ypersfaíra or hypersphere in English), and the latter densities by the Greek letter $\rho$.

In Bayesian inference, the determination of the Bayes prior is referred to as an inverse problem, and Jaynes’ data-constrained maximal entropy priors provide a principled solution to this inverse problem [5,6,7]. The maximal entropy is reached by minimizing the Gibbs potential, and the solution to this optimization problem requires the determination of the empirical density entropic convex dual. See Le Blanc [8] for an extensive review on the subject, with references. As derived from translated normal distributions, the classical noncentral distributions all have a submodular decomposition for their translation factor of the form

$$T\left(r\right|r_o)=E\left(r\right|r_o)e^{-r_0^2/2},$$

with $E\left(r\right|r_o)$ being a generalized hypergeometric function, a decomposition which does not readily provide regrouping of terms of similar order in the noncentrality parameter $r_o$. Conversely, the translation factors for the noncentral t, the normal N, F, and ${\chi }^2$ distributions, as derived from translated uniform distributions on unit radius hyperspheres, are found to be generating functions for the Gegenbauer, Hermite, Laguerre, and Jacobi orthogonal polynomial families ${\left\{P_n\right\}}_{n=0}^{\infty }$, respectively, and intrinsic properties of the orthogonal polynomials [9] allow for regrouping all terms of similar order in the noncentrality parameter $r_o,$ that is,

$$T\left(r\right|r_o)=\sum _{n=0}^{\infty }{c}_n{P}_n\left(r\right)r_o^n,$$

with the constants $c_n$ provided hereafter. To the best of our knowledge, the derivation of these translation factors and their identification as orthogonal polynomial family-generating functions has not been carried out before. As a consequence, one can expand the entropic convex duals on a generally small number of low-order orthogonal polynomials, an approach which greatly curtails the computational cost of determining the duals and, thus, the Jaynes–Gibbs priors. To the best of our knowledge, the expansion and discretization of the convex duals over orthogonal polynomial bases has not been proposed before. We adopt in this manuscript the convention that the polynomial family-defining weight functions should be provided by the corresponding normalized central distributions, a convention which results in much simplified expressions for the norm of the orthogonal polynomials.

In parametric Bayesian modeling, the translation factor $T\left(r\right|r_o)$ can be weighted by a Bayesian prior $\pi \left(r_o\right)$ to obtain the Bayes factor

$$\mathrm{BF}\left(r\right)=\int T\left(r\right|r_o)\pi \left(r_o\right)\mathrm{d}{r}_o$$

for the generic superposition density

$$\rho \left(r\right)=\int \rho \left(r\right|r_o)\pi \left(r_o\right)\mathrm{d}{r}_o=\int T\left(r\right|r_o)\rho \left(r\right|0)\pi \left(r_o\right)\mathrm{d}{r}_o=BF\left(r\right)\rho \left(r\right|0).$$

We will, in this manuscript, identify the normalized central distribution $\rho \left(r\right|0)$ with the orthogonal polynomial family-defining weight $w\left(r\right)\equiv \rho \left(r\right|0)$ [9], with, as a result, the rewriting of any generic density $\rho \left(r\right)$ as

$$\rho \left(r\right)=BF\left(r\right)w\left(r\right),$$

that is, the Bayes factor $BF\left(r\right)$ can stand as a substitute to the generic density $\rho \left(r\right).$ With $p\left(r\right)={\int }^r\rho \left(r^{\prime }\right|0)\mathrm{d}{r}^{\prime }$ standing for the cumulative density function of the respective central distributions $\rho \left(r\right|0)$, or, equivalently, the p-value of the null hypothesis statistical testing (NHST) procedure, we have that $\rho \left(p\right)=BF\left(p\right)$. That is, the Bayes factor $BF\left(p\right)$ stands for the generally nonuniform p-value distribution $\rho \left(p\right)$ of the above generic superposition density (see Le Blanc [10] for details of a proof), to be contrasted with the NHST framework which only considers the central distribution with its uninformative uniform p-value distribution. Now, if one’s goal is only to model the density $BF\left(p\right)$ or to compute the associated local false discovery rate $\mathrm{fdr}\left(p\right)=1/(1+BF(p\left)\right)$ [11], one then only needs prior moments to carry on with the modelization. Invoking the moment problem [12] and the duality principle [13], prior moments in parametric space will be shown to be readily provided by Bayes factor $BF\left(r\right)$ orthogonal polynomial expansion coefficients in random variable space—that is,

$${\int }_{{\mathcal{R}}_o}\pi \left(r_o\right)r_o^n\mathrm{d}{r}_o=\frac{1}{c_n{\parallel P_n\parallel }^2}{\int }_{\mathcal{R}}{P}_n\left(r^{\prime }\right)BF\left(r^{\prime }\right)w\left(r^{\prime }\right)\mathrm{d}{r}^{\prime }.$$

The paper is organized as follows. In Section 2, we review, for the sake of completeness, the derivation of the classical univariate noncentral distributions as derived from translated normal distributions. In Section 3, the univariate ultraspherical noncentral t, the normal N, F, and ${\chi }^2$ distributions are derived from translated uniform density distributions on unit radius hyperspheres, and are shown to be expressible as products of their central distribution times specific generating functions for the Gegenbauer, Hermite, Jacobi, and Laguerre orthogonal polynomial families, respectively. We argue, in Section 4, that the determination of the Gibbs priors in terms of empirical densities’ entropic convex duals is much simplified when these duals are expanded on a small number of low-order orthogonal polynomials. We also discuss how prior moments in parametric space are directly provided by the Bayes factor orthogonal polynomial expansion coefficients in random variable space. Section 5 and Section 6 are devoted to applications in genomics and geophysics, respectively.

2. The Classical Noncentral Distributions

The four central F, ${\chi }^2$, t, and normal N distributions can be given by

$$\begin{array}{}\mathrm{(1a)}& {\upsilon }_F^{({\nu }_1,{\nu }_2)}\left(\theta \right|0)& =2\frac{\mathsf{\Gamma }(\frac{{\nu }_1+{\nu }_2}{2})}{\mathsf{\Gamma }\left(\frac{{\nu }_1}{2}\right)\mathsf{\Gamma }\left(\frac{{\nu }_2}{2}\right)}{cos}^{{\nu }_1-1}\theta {sin}^{{\nu }_2-1}\theta ,0\le cos\theta <1,\hfill \mathrm{(1b)}& {\rho }_{{\chi }^2}^{\left({\nu }_1\right)}\left(r\right|0)& =\frac{1}{2^{{\nu }_1/2}\mathsf{\Gamma }\left(\frac{{\nu }_1}{2}\right)}{r}^{({\nu }_1-2)/2}{e}^{-r/2},0\le r<\infty ,\hfill \mathrm{(1c)}& {\rho }_t^{\left({\nu }_2\right)}\left(\theta \right|0)& =\frac{\mathsf{\Gamma }\left(\frac{{\nu }_2+1}{2}\right)}{\mathsf{\Gamma }\left(\frac{1}{2}\right)\mathsf{\Gamma }\left(\frac{{\nu }_2}{2}\right)}{sin}^{{\nu }_2-1}\theta ,-1<cos\theta <1,\hfill \mathrm{(1d)}& {\rho }_N\left(r\right|0)& =\frac{1}{\sqrt{2\pi }}{e}^{-\frac{1}{2}{r}^2},-\infty <r<\infty ,\hfill \end{array}$$

respectively, where we have redefined the conventional t and F statistics so that

$$t=\sqrt{{\nu }_2}\frac{cos\theta }{sin\theta },-1<cos\theta <1,$$

(2)

and

$$F=\frac{{\nu }_2}{{\nu }_1}\frac{{cos}^2\theta }{{sin}^2\theta },0\le cos\theta <1,$$

(3)

so as to re-scale the corresponding t and F distributions on the specified finite compact domains. As discussed in the Introduction, the above central t and F distributions are projections of uniform distributions on the unit radius ${\mathcal{S}}^{{\nu }_2}$ and ${\mathcal{S}}^{{\nu }_1+{\nu }_2}$ hyperspheres, respectively. The reasons for such rescaling will become apparent when we link in Section 3 the ultraspherical noncentral t and F distributions with the Gegenbauer and Jacobi orthogonal polynomial families, respectively, families which both have similar finite compact support domains. It is straightforward to verify that the above central t and F distributions converge with the central normal N and ${\chi }^2$ distributions, respectively, in the limit ${\nu }_2\to \infty$: see the discussions leading to Equations (19) and (36) below. In geometric terms, the limit ${\nu }_2\to \infty$ allows one to restrict the study of the central t and F distributions on a restricted domain around the angle $\theta =\pi /2$, around which these distributions will concentrate their respective distribution weights in concordance with the notion that a high-dimensional hypersphere concentrates its surface on a narrow equatorial band at its equator [3]. In this limit, with the respective variable substitutions $r=\sqrt{{\nu }_2-1}cos\theta$ and $r=({\nu }_2-2){cos}^2\theta$, the normal N and ${\chi }^2$ respective central distributions’ definition domains consequently extend to infinity as $cos\theta \to 1$, as reflected by the respective distribution domains provided above.

The noncentral t distribution has been classically derived by considering the ratio of a random variable which distributes according to a noncentral normal distribution $N(\delta ,1)$ with non-vanishing mean $\delta$, over that of a random variable distributed according to a central ${\chi }_{{\nu }_2}$ distribution with ${\nu }_2$ degrees of freedom. Similarly, the noncentral F distribution is classically derived by considering the ratio of a noncentral ${\chi }^2$ random variable with noncentrality parameter $\mathsf{\Lambda }$ and ${\nu }_1$ degrees of freedom over that of a central ${\chi }^2$ random variable with ${\nu }_2$ degrees of freedom. See, for example, Walck et al. [14] for explicit steps for such derivations. As a result, as discussed in the Introduction, the classical distributions have a submodular decomposition for their translation factor of the form $T\left(r\right|r_o)=E\left(r\right|r_o)e^{-r_0^2/2}$, with $E\left(r\right|r_o)$ being a generalized hypergeometric function [15]. Upon recalling that the F and ${\chi }^2$ noncentral parameter $\mathsf{\Lambda }$ stands for the square of the noncentral parameter $\delta$ for the noncentral t and normal distributions, we have that

$$\begin{array}{}\mathrm{(4a)}& {\rho }_F^{({\nu }_1,{\nu }_2)}\left(\theta \right|\mathsf{\Lambda })& =E_F^{({\nu }_1,{\nu }_2)}\left(\theta \right|\mathsf{\Lambda })e^{-\mathsf{\Lambda }/2}{\rho }_F^{({\nu }_1,{\nu }_2)}\left(\theta \right|0),0\le cos\theta <1,0\le \mathsf{\Lambda }<\infty ,\hfill \mathrm{(4b)}& {\rho }_{{\chi }^2}^{\left({\nu }_1\right)}\left(r\right|\mathsf{\Lambda })& =E_{{\chi }^2}^{\left({\nu }_1\right)}\left(r\right|\mathsf{\Lambda })e^{-\mathsf{\Lambda }/2}{\rho }_{{\chi }^2}^{\left({\nu }_1\right)}\left(r\right|0),0\le r<\infty ,0\le \mathsf{\Lambda }<\infty ,\hfill \mathrm{(4c)}& {\rho }_t^{\left({\nu }_2\right)}\left(\theta \right|\delta )& =E_t^{\left({\nu }_2\right)}\left(\theta \right|\delta )e^{-{\delta }^2/2}{\rho }_t^{\left({\nu }_2\right)}\left(\theta \right|0),-1<cos\theta <1,-\delta <\delta <\infty ,\hfill \mathrm{(4d)}& {\rho }_N\left(r\right|\delta )& =E_N\left(r\right|\delta )e^{-{\delta }^2/2}{\rho }_N\left(r\right|0),-\infty <r<\infty ,-\delta <\delta <\infty ,\hfill \end{array}$$

with the respective factors’ $E\left(r\right|r_o)$ McLaurin expansions — computed following the steps elaborated by Walck et al. [14]—given by

$$\begin{array}{}\mathrm{(5a)}& E_F^{({\nu }_1,{\nu }_2)}\left(\theta \right|\mathsf{\Lambda })& =\frac{\mathsf{\Gamma }({\nu }_1/2)}{\mathsf{\Gamma }(1/2)}\sum _{j\mathrm{even}\ge 0}^{\infty }\frac{\mathsf{\Gamma }\left(\right(j+1)/2)}{\mathsf{\Gamma }((j+{\nu }_1)/2)}\frac{\mathsf{\Gamma }((j+{\nu }_1+{\nu }_2)/2)}{\mathsf{\Gamma }(({\nu }_1+{\nu }_2)/2)}\frac{{(\sqrt{2\mathsf{\Lambda }}cos\theta )}^j}{j!}\hfill \mathrm{(5b)}& & ={}_{1}F_1(\frac{{\nu }_1+{\nu }_2}{2};\frac{{\nu }_1}{2};\frac{\mathsf{\Lambda }{cos}^2\theta }{2}),\mathrm{(5c)}& E_{{\chi }^2}^{\left({\nu }_1\right)}\left(r\right|\mathsf{\Lambda })& =\frac{\mathsf{\Gamma }({\nu }_1/2)}{\mathsf{\Gamma }(1/2)}\sum _{j\mathrm{even}\ge 0}^{\infty }\frac{\mathsf{\Gamma }\left(\right(j+1)/2)}{\mathsf{\Gamma }((j+{\nu }_1)/2)}\frac{{\left(\mathsf{\Lambda }r\right)}^{j/2}}{j!}={\mathcal{I}}^{({\nu }_1-2)/2}\left(\mathsf{\Lambda }r\right)\hfill \mathrm{(5d)}& & ={}_{0}F_1(;\frac{{\nu }_1}{2};\frac{\mathsf{\Lambda }r}{4}),\mathrm{(5e)}& E_t^{\left({\nu }_2\right)}\left(\theta \right|\delta )& =\sum _{j=0}^{\infty }\frac{\mathsf{\Gamma }((j+1+{\nu }_2)/2)}{\mathsf{\Gamma }((1+{\nu }_2)/2)}\frac{{(\sqrt{2}\delta cos\theta )}^j}{j!}\hfill & & ={}_{1}F_1(\frac{{\nu }_2+1}{2};\frac{1}{2};\frac{{\delta }^2{cos}^2\theta }{2})+\sqrt{2}\delta cos\theta \frac{\mathsf{\Gamma }\left(\frac{{\nu }_2+2}{2}\right)}{\mathsf{\Gamma }\left(\frac{{\nu }_2+1}{2}\right)}{}_{1}F_1(\frac{{\nu }_2+2}{2};\frac{3}{2};\frac{{\delta }^2{cos}^2\theta }{2}),\mathrm{(5f)}& E_N\left(r\right|\delta )& =\sum _{j=0}^{\infty }\frac{{\left(\delta r\right)}^j}{j!}={}_{0}F_0(;;r)=e^{\delta r}\hfill \mathrm{(5g)}& & ={}_{0}F_1(;\frac{1}{2};\frac{{\delta }^2{r}^2}{4})+\delta r{}_{0}F_1(;\frac{3}{2};\frac{{\delta }^2{r}^2}{4})=cosh\delta r+sinh\delta r=e^{\delta r},\end{array}$$

where

$${}_{p}F_q(a_1,\dots ,a_p;b_1,\dots ,b_q;x)\sum _{n=0}^{\infty }\frac{{\left(a_1\right)}_n\dots {\left(a_p\right)}_n}{{\left(b_1\right)}_n\dots {\left(b_q\right)}_n}\frac{x^n}{n!}.$$

(6)

stands for the generalized hypergeometric function, and where ${\mathcal{I}}^{({\nu }_1-2)/2}\left(\mathsf{\Lambda }r\right)$ stands for the normalized modified Bessel function [16]. The respective first expansions for the multiplicative factors E are successively simpler versions of the Maclaurin expansion for the noncentral F distribution multiplicative factor $E_F^{({\nu }_1,{\nu }_2)}$, such that

$$E_{{\chi }^2}^{\left({\nu }_1\right)}\underset{{\nu }_2\to \infty }{\leftarrow }{E}_F^{({\nu }_1,{\nu }_2)}\underset{{\nu }_1=1}{\to }{E}_t^{\left({\nu }_2\right)}\underset{{\nu }_2\to \infty }{\to }{E}_N,$$

(7)

with the understanding that the respective Maclaurin summations involve either only the even integers for the F and ${\chi }^2$ cases, or both odd and even integers for the t and normal cases. It is straightforward to verify that the classical noncentral t and F distributions converge with the noncentral normal and ${\chi }^2$ distributions, respectively, in the limit ${\nu }_2\to \infty$. The expressions for the classical noncentral distributions do not readily regroup terms of similar order in their noncentrality parameter, contrarily to the ultraspherical noncentral distributions, which do so, as is discussed next.

3. The Ultraspherical Noncentral Distributions

Using geometrical arguments, the ultraspherical noncentral t-distribution for the t-statistic (2) on the hypersphere ${\mathcal{S}}^{{\nu }_2}$ was shown by Le Blanc [10] to be given by

$${\upsilon }_t^{\left({\nu }_2\right)}\left(\theta \right|\delta )=T_t^{\left({\nu }_2\right)}\left(\theta \right|\delta ){\upsilon }_t^{\left({\nu }_2\right)}\left(\theta \right|0),$$

(8)

where ${\upsilon }_t^{\left({\nu }_2\right)}\left(\theta \right|0)$ is identical to the central t-distribution ${\rho }_t^{\left({\nu }_2\right)}\left(\theta \right|0)$ given in Equation (1), where the multiplicative distribution translation factor $T_t^{\left({\nu }_2\right)}\left(\theta \right|\delta )$ is given by

$$T_t^{\left({\nu }_2\right)}\left(\theta \right|\delta )=\frac{1-cos\theta cos{\theta }_{\delta }}{{[1-2cos\theta cos{\theta }_{\delta }+{cos}^2{\theta }_{\delta }]}^{({\nu }_2+1)/2}},$$

(9)

and where

$$cos{\theta }_{\delta }=\frac{\delta /\sqrt{{\nu }_2}}{\sqrt{{(\delta /\sqrt{{\nu }_2})}^2+1}}=\frac{\delta }{\sqrt{{\delta }^2+{\nu }_2}},0<{\theta }_{\delta }<\pi ,$$

(10)

in terms of the noncentrality parameter $\delta ,$$-\infty <\delta <\infty .$ As it turns out, the translation term $T_t^{\left({\nu }_2\right)}\left(\theta \right|\delta )$ is a generating function for the ultraspherical or, equivalently, Gegenbauer polynomials: redefining the variables so that $x=cos\theta$, $z=cos{\theta }_{\delta }$, and $b=({\nu }_2-1)/2$, we have that

$$\begin{array}{ccc}\hfill T_t^{\left(b\right)}\left(x\right|z)& =& \frac{1-xz}{{[1-2xz+z^2]}^{b+1}}=\sum _{n=0}^{\infty }\frac{2b+n}{2b}{C}_n^{\left(b\right)}\left(x\right)z^n\hfill \\ & =& {(1-xz)}^{-(2b+1)}{}_2{F}_1\left(\frac{2b+1}{2},\frac{2b+2}{2};b+\frac{1}{2};-\frac{(1-x^2)z^2}{{(1-xz)}^2}\right),\hfill \end{array}$$

(11)

where $C_n^{\left(b\right)}\left(x\right)$ are Gegenbauer polynomials with the explicit representation [4]

$$C_n^{\left(b\right)}\left(x\right)=\sum _{k=0}^{⌊n/2⌋}\frac{{(-1)}^k{\left(b\right)}_{n-k}}{k!(n-2k)!}{\left(2x\right)}^{n-2k},$$

(12)

where the Pochhammer symbol ${\left(x\right)}_n$ is defined by the equality

$${\left(x\right)}_n=\frac{\mathsf{\Gamma }(x+n)}{\mathsf{\Gamma }\left(x\right)}=x(x+1)(x+2)\cdots (x+n-1),$$

(13)

where $⌊y⌋$, the floor of y, is given by the lowest integer such that $y-1<⌊y⌋\leqq y$, and where the last equality given in terms of the hypergeometric function ${}_{2}F_1$ can be deduced from Rainville [17] (Equation (144.8)). The Gegenbauer polynomials are orthogonal:

$${\int }_{x=-1}^1{C}_m^{\left(b\right)}\left(x\right)C_n^{\left(b\right)}\left(x\right)w_t^{\left(b\right)}\left(x\right)\mathrm{d}x={\delta }_{m,n}{\parallel C_n^{\left(b\right)}\parallel }^2$$

(14)

with respect to the weight function

$$w_t^{\left(b\right)}\left(x\right)=\frac{\mathsf{\Gamma }(b+1)}{\mathsf{\Gamma }\left(\frac{1}{2}\right)\mathsf{\Gamma }(b+\frac{1}{2})}{(1-x^2)}^{b-1/2},$$

(15)

which is identical to the central t distribution (1), except for the change of the variable $x=cos\theta$. With this weight function normalization—which, note, differs from the usual weight function for the Gegenbauer polynomials [4]—the norm of the Gegenbauer polynomials simplifies to

$$\parallel C_n^{\left(b\right)}{\parallel }^2=\frac{b}{b+n}\frac{{\left(2b\right)}_n}{n!},$$

(16)

with, in particular, $\parallel C_0^{\left(b\right)}{\parallel }^2=1$. Note also that Equation (11) could be used to define a generalization $T_n^{\left(b\right)}\left(x\right)$ for the Chebyshev polynomials of the first kind, $T_n^{(b=0)}\left(x\right)$, with

$$\frac{1-xz}{{[1-2xz+z^2]}^{b+1}}=\sum _{n=0}^{\infty }{T}_n^{\left(b\right)}\left(z\right)z^n=\sum _{n=0}^{\infty }\frac{2b+n}{2b}{C}_n^{\left(b\right)}\left(x\right)z^n,$$

(17)

which encompasses the defining equation

$$\frac{1-xz}{[1-2xz+z^2]}=\sum _{n=0}^{\infty }{T}_n^{\left(0\right)}\left(x\right)z^n=\underset{b\to 0}{lim}\sum _{n=0}^{\infty }\frac{n}{2b}{C}_n^{(\Delta )}\left(x\right)z^n$$

(18)

for the Chebyshev polynomials of the first kind. Setting $x=r/{\left(2b\right)}^{1/2}$ and $z=\delta /{\left(2b\right)}^{1/2}$, the central t distribution converges with the normal central distribution

$${\upsilon }_N\left(r\right|0)=\frac{1}{\sqrt{2\pi }}{e}^{-\frac{1}{2}{r}^2}$$

(19)

in the limit $b\to \infty$. Applying the same limit to the Gegenbauer polynomials and their generating function (11), we have that

$$\begin{array}{ccc}\hfill \underset{b\to \infty }{lim}{T}_t^{\left(b\right)}\left(x\right|z)& =& \underset{b\to \infty }{lim}\frac{1-xz}{{[1-2xz+z^2]}^{b+1}}\hfill \\ & =& \underset{b\to \infty }{lim}\sum _{n=0}^{\infty }\frac{2b+n}{2b}{C}_n^{\left(b\right)}\left(\frac{r}{{\left(2b\right)}^{1/2}}\right){\left(\frac{\delta }{{\left(2b\right)}^{1/2}}\right)}^n\hfill \\ & =& \sum _{n=0}^{\infty }{\mathit{He}}_n\left(r\right)\frac{{\delta }^n}{n!}=e^{r\delta -{\delta }^2/2}=T_N\left(r\right|\delta ),\hfill \end{array}$$

(20)

where we have used the limit result [4]

$$\underset{b\to \infty }{lim}\frac{1}{{\left(2b\right)}^{n/2}}{C}_n^{\left(b\right)}\left(\frac{r}{{\left(2b\right)}^{1/2}}\right)=\frac{{\mathit{He}}_n\left(r\right)}{n!}.$$

(21)

The multiplicative factor $T_N\left(r\right|\delta )$ imparts a non-vanishing mean $\delta$ to the zero mean normal distribution, since

$$T_N\left(r\right|\delta ){\upsilon }_N\left(r\right|0)=e^{\delta r-{\delta }^2/2}\frac{1}{\sqrt{2\pi }}{e}^{-\frac{1}{2}{r}^2}=\frac{1}{\sqrt{2\pi }}{e}^{-\frac{1}{2}{(r-\delta )}^2}={\upsilon }_N\left(r\right|\delta ),$$

(22)

where one recognizes the generating function $e^{\delta r-{\delta }^2/2}$ for the Hermite polynomials ${\mathit{He}}_n\left(r\right)$, with the explicit representation [4]

$${\mathit{He}}_n\left(r\right)=n!\sum _{\ell =0}^{⌊n/2⌋}\frac{{(-1)}^{\ell }{r}^{n-2\ell }}{2^{\ell }\ell !(n-2\ell )!}.$$

(23)

The latter are orthogonal with respect to their defining weight function

$$w_N\left(r\right)={\upsilon }_N\left(r\right|0)=\frac{1}{\sqrt{2\pi }}{e}^{-\frac{1}{2}{r}^2},$$

(24)

with norm

$${\int }_{r=-\infty }^{\infty }{\mathit{He}}_m\left(r\right){\mathit{He}}_n\left(r\right)w_N\left(r\right)\mathrm{d}r={\delta }_{m,n}n!,$$

(25)

and with $\parallel {\mathit{He}}_0{\parallel }^2=1$ in particular.

Using geometrical arguments, the ultraspherical noncentral F-distribution for the F-statistic (3) was shown by Le Blanc [10] to be given by the integral representation

$${\upsilon }_F^{({\nu }_1,{\nu }_2)}\left(\theta \right|\mathsf{\Lambda })=T_F^{({\nu }_1,{\nu }_2)}\left(\theta \right|\mathsf{\Lambda }){\upsilon }_F^{({\nu }_1,{\nu }_2)}\left(\theta \right|0),$$

(26)

where ${\upsilon }_F^{({\nu }_1,{\nu }_2)}\left(\theta \right|0)$ is identical to the central F distribution ${\rho }_F^{({\nu }_1,{\nu }_2)}\left(\theta \right|0)$ given in Equation (1), where the multiplicative distribution translation factor $T_F^{({\nu }_1,{\nu }_2)}\left(\theta \right|\mathsf{\Lambda })$ is given by the integral

$$T_F^{({\nu }_1,{\nu }_2)}\left(\theta \right|\mathsf{\Lambda })={\int }_{\psi =0}^{\pi }{T}_F^{({\nu }_1,{\nu }_2)}(\theta ,\psi |\mathsf{\Lambda })\mathrm{d}\psi$$

(27)

with integrand

$$T_F^{({\nu }_1,{\nu }_2)}(\theta ,\psi |\mathsf{\Lambda })=\frac{(1-cos\theta cos\psi cos{\theta }_{\mathsf{\Lambda }})}{{[1-2cos\theta cos\psi cos{\theta }_{\mathsf{\Lambda }}+{cos}^2{\theta }_{\mathsf{\Lambda }}]}^{({\nu }_1+{\nu }_2)/2}}{\upsilon }_t^{({\nu }_1-1)}\left(\psi \right|0),$$

(28)

and where

$$cos{\theta }_{\mathsf{\Lambda }}=\sqrt{\mathsf{\Lambda }/(\mathsf{\Lambda }+{\nu }_2)},0\le \mathsf{\Lambda }<\infty ,0\le cos{\theta }_{\mathsf{\Lambda }}<1,$$

(29)

in terms of the noncentral parameter $\mathsf{\Lambda }$. The special case ${\nu }_1=1$ is given by

$${\upsilon }_F^{({\nu }_1=1,{\nu }_2)}\left(\theta \right|\mathsf{\Lambda })=\sum _{cos{\theta }^{\prime }=[cos\theta ,-cos\theta ]}\frac{1-cos{\theta }^{\prime }cos{\theta }_{\mathsf{\Lambda }}}{{[1-2cos{\theta }^{\prime }cos{\theta }_{\mathsf{\Lambda }}+{cos}^2{\theta }_{\mathsf{\Lambda }}]}^{({\nu }_2+1)/2}}{\upsilon }_t^{\left({\nu }_2\right)}\left(\theta \right|0).$$

(30)

Setting $a=({\nu }_1-2)/2,b=({\nu }_2-2)/2,x=cos\theta ,\xi =cos\psi ,z=cos{\theta }_{\mathsf{\Lambda }}$, and making use of the Gegenbauer polynomial-generating function expansion (11) together with the corresponding polynomials’ explicit expression (12), the integration in (27) can be carried out: we find

$$\begin{array}{ccc}\hfill T_F^{(a,b)}\left(x\right|z)& =& {\int }_{\xi =-1}^1\frac{1-x\xi z}{{[1-2x\xi z+z^2]}^{a+b+2}}{w}_t^{\left(a\right)}\left(\xi \right)\mathrm{d}\xi \hfill \\ & =& \sum _{n=0}^{\infty }\frac{a+b+1+n}{a+b+1}{\mathcal{F}}_n^{(a,b)}\left(x^2\right)z^{2n},\hfill \\ & =& \sum _{n=0}^{\infty }{(-1)}^n\frac{{(a+b+2)}_n}{{(a+1)}_n}{P}_n^{(a,b)}(1-2x^2)z^{2n}\hfill \\ & =& {(1+z^2)}^{-(a+b+2)}{}_{2}F_1(\frac{1}{2}(a+b+2),\frac{1}{2}(a+b+3);a+1;\frac{4x^2{z}^2}{{(1+z^2)}^2}),\hfill \end{array}$$

(31)

where the polynomials ${\mathcal{F}}_n^{(a,b)}\left(x^2\right)$ which one related to the Jacobi polynomials [4], through

$${\mathcal{F}}_n^{(a,b)}\left(x^2\right)={(-1)}^n\frac{{(a+b+1)}_n}{{(a+1)}_n}{P}_n^{(a,b)}(1-2x^2)$$

(32)

are provided with the explicit representation

$${\mathcal{F}}_n^{(a,b)}\left(x^2\right)=\sum _{k=0}^n\frac{{(-1)}^k}{k!}\frac{{(a+b+1)}_{2n-k}}{{(a+1)}_{n-k}}\frac{{\left(x^2\right)}^{n-k}}{(n-k)!},$$

(33)

and where the last equality—deduced from Rainville [17] (Equation (132.10))—provides us with a generating function for the polynomials in terms of the hypergeometric function ${}_{2}F_1$. The ${\mathcal{F}}_n^{(a,b)}\left(x^2\right)$ polynomials are orthogonal with respect to the weight function

$$w_F^{(a,b)}\left(x\right)=2\frac{\mathsf{\Gamma }(a+b+1)}{\mathsf{\Gamma }(a+1)\mathsf{\Gamma }\left(b\right)}{\left(x^2\right)}^{a+\frac{1}{2}}{(1-x^2)}^b,0\le x\le 1,$$

(34)

which is identical to the central F distribution (1), except for the change of the variable $x=cos\theta ,$ with norm

$$\parallel {\mathcal{F}}_n^{(a,b)}{\parallel }^2=\frac{a+b+1}{a+b+1+2n}\frac{{(a+b+1)}_n{(b+1)}_n}{{(a+1)}_{n}n!},$$

(35)

and with $\parallel {\mathcal{F}}_0^{(a,b)}{\parallel }^2=1$ in particular. Equation (31) can be verified to be valid for the special case (30) with $a=-1/2$ (${\nu }_1=1$). In the limit ${\nu }_2\to \infty$, the central F distribution converges with the central ${\chi }^2$ distribution

$${\upsilon }_{{\chi }^2}^{\left({\nu }_1\right)}\left(r\right|0)=\frac{1}{2^{{\nu }_1/2}\mathsf{\Gamma }\left(\frac{{\nu }_1}{2}\right)}{r}^{({\nu }_1-2)/2}{e}^{-r/2},0\le r<\infty ,$$

(36)

which, with $a=({\nu }_1-2)/2$, can be used as a normalized defining weight function,

$$w_{{\chi }^2}^{\left(a\right)}\left(r\right)=\frac{1}{2^{a+1}}\frac{1}{\mathsf{\Gamma }(a+1)}{r}^a{e}^{-r/2},a>-1,$$

(37)

for the Laguerre orthogonal polynomial family. The Laguerre polynomials can be given the explicit representation [4]

$$L_n^{\left(a\right)}\left(y\right)=\sum _{\ell =0}^n{(-1)}^{\ell }\frac{{(a+\ell +1)}_{n-\ell }}{(n-\ell )!}\frac{y^{\ell }}{\ell !},0\le y<\infty .$$

(38)

Since Equation (32) provides the noncentral F distribution translation factor $T_F^{(a,b)}\left(x\right|z)$ with an expansion in terms of the Jacobi polynomials $P_n^{(a,b)}(1-2x^2)$, and since the limit result

$$\underset{b\to \infty }{lim}{P}_n^{(a,b)}\left(1-\frac{2y}{b}\right)=L_n^{\left(a\right)}\left(y\right)$$

(39)

holds [4], one verifies with $x^2=r/2b$ and $z^2=\mathsf{\Lambda }/2b$ that

$$\begin{array}{ccc}\hfill \underset{b\to \infty }{lim}{T}_F^{(a,b)}\left(x\right|z)& =& T_{{\chi }^2}^{\left(a\right)}\left(r\right|\mathsf{\Lambda })={\mathcal{I}}^{\left(a\right)}\left(\mathsf{\Lambda }r\right)e^{-\mathsf{\Lambda }/2}{}_{0}F_1(;a+1;\frac{\mathsf{\Lambda }r}{4})e^{-\mathsf{\Lambda }/2}\hfill \\ & =& \sum _{n=0}^{\infty }\frac{{(-1)}^n}{{(a+1)}_n}{L}_n^{\left(a\right)}(r/2){(\mathsf{\Lambda }/2)}^n,\hfill \end{array}$$

(40)

which is a lesser-known expression for the noncentral ${\chi }^2$ distribution translation factor first derived by Tiku [18]. Under normalization (37) for their defining weight factor, the norm of the Laguerre polynomials $L_n^{\left(a\right)}(r/2)$ is given by

$${\int }_{r=0}^{\infty }{L}_m^{\left(a\right)}(r/2)L_n^{\left(a\right)}(r/2)w_{{\chi }^2}^{\left(a\right)}\left(r\right)\mathrm{d}r={\delta }_{m,n}{\parallel L_n^{\left(a\right)}\parallel }^2{\delta }_{m,n}=\frac{{(a+1)}_n}{n!}$$

(41)

with, in particular, $\parallel L_0^{\left(a\right)}{\parallel }^2=1.$

To summarize, the four ultraspherical noncentral F, ${\chi }^2$, t, and the normal N distributions are given by the product of translation factors T, which take the form of generating functions for specific orthogonal polynomial families, times the corresponding normalized central distributions which also stand for the corresponding polynomial family-defining weight functions. Thus, we have that

$$\begin{array}{}\mathrm{(42a)}& {\upsilon }_F^{(a,b)}\left(x\right|z)& =T_F^{(a,b)}\left(x\right|z)w_F^{(a,b)}\left(x\right),0\le x<1,0\le z<1,\hfill \mathrm{(42b)}& {\upsilon }_{{\chi }^2}^{\left(a\right)}\left(r\right|\mathsf{\Lambda })& =T_{{\chi }^2}^{\left(a\right)}\left(r\right|\mathsf{\Lambda })w_{{\chi }^2}^{\left(a\right)}\left(r\right),0\le r<\infty ,0\le \mathsf{\Lambda }<\infty ,\hfill \mathrm{(42c)}& {\upsilon }_t^{\left(b\right)}\left(x\right|z)& =T_t^{\left(b\right)}\left(x\right|z)w_t^{\left(b\right)}\left(x\right),-1<x<1,-1<z<1,\hfill \mathrm{(42d)}& {\upsilon }_N\left(r\right|\delta )& =T_N\left(r\right|\delta )w_N\left(r\right),-\infty <r<\infty ,-\infty <\delta <\infty ,\hfill \end{array}$$

with the normalized central distributions

$$\begin{array}{}\mathrm{(43a)}& w_F^{(a,b)}\left(x\right)& =2\frac{\mathsf{\Gamma }(a+b+1)}{\mathsf{\Gamma }(a+1)\mathsf{\Gamma }\left(b\right)}{\left(x^2\right)}^{a+\frac{1}{2}}{(1-x^2)}^b,\hfill \mathrm{(43b)}& w_{{\chi }^2}^{\left(a\right)}\left(r\right)& =\frac{1}{2^{a+1}}\frac{1}{\mathsf{\Gamma }(a+1)}{r}^a{e}^{-r/2},\hfill \mathrm{(43c)}& w_t^{\left(b\right)}\left(x\right)& =\frac{\mathsf{\Gamma }(b+1)}{\mathsf{\Gamma }\left(\frac{1}{2}\right)\mathsf{\Gamma }(b+\frac{1}{2})}{(1-x^2)}^{b-1/2}\hfill \mathrm{(43d)}& w_N\left(r\right)& ={\upsilon }_N\left(r\right|0)=\frac{1}{\sqrt{2\pi }}{e}^{-\frac{1}{2}{r}^2},\hfill \end{array}$$

providing the polynomial family-defining weights, and with the corresponding generating functions and orthogonal polynomial expansions given by

$$\begin{array}{}\mathrm{(44a)}& T_F^{(a,b)}\left(x\right|z)& ={(1+z^2)}^{-(a+b+2)}{}_{2}F_1(\frac{1}{2}(a+b+2),\frac{1}{2}(a+b+3);a+1;\frac{4x^2{z}^2}{{(1+z^2)}^2})\hfill & & =\sum _{n=0}^{\infty }{(-1)}^n\frac{{(a+b+2)}_n}{{(a+1)}_n}{P}_n^{(a,b)}(1-2x^2)z^{2n}=\sum _{n=0}^{\infty }\frac{a+b+1+n}{a+b+1}{\mathcal{F}}_n^{(a,b)}\left(x^2\right)z^{2n},\mathrm{(44b)}& T_{{\chi }^2}^{\left(a\right)}\left(r\right|\mathsf{\Lambda })& {}_{0}F_1(;a+1;\frac{\mathsf{\Lambda }r}{4})e^{-\mathsf{\Lambda }/2}={\mathcal{I}}^{\left(a\right)}\left(\mathsf{\Lambda }r\right)e^{-\mathsf{\Lambda }/2}=\sum _{n=0}^{\infty }\frac{{(-1)}^n}{{(a+1)}_n}{L}_n^{\left(a\right)}(r/2){(\mathsf{\Lambda }/2)}^n,\hfill \mathrm{(44c)}& T_t^{\left(b\right)}\left(x\right|z)& ={(1-xz)}^{-(2b+1)}{}_{2}F_1(\frac{2b+1}{2},\frac{2b+2}{2};b+\frac{1}{2};-\frac{(1-x^2)z^2}{{(1-xz)}^2})\hfill & & =\frac{1-xz}{{[1-2xz+z^2]}^{b+1}}=\sum _{n=0}^{\infty }\frac{2b+n}{2b}{C}_n^{\left(b\right)}\left(x\right)z^n,\hfill & \mathrm{(44d)}& T_N\left(r\right|\delta )& {}_{0}F_1(;;\delta r)e^{-{\delta }^2/2}=e^{r\delta -{\delta }^2/2}=\sum _{n=0}^{\infty }{\mathit{He}}_n\left(r\right)\frac{{\delta }^n}{n!}.\hfill \end{array}$$

We conclude this section by stressing the fact that the random variable space for the ultraspherical noncentral distributions are translated hyperspheres rather than translated normal distributions, as is assumed for the classical noncentral distributions. The tide model developed in Section 6 provides a concrete example of such a distribution on a translated sphere. As more extensively argued in [15], the ultraspherical and classical noncentral F and t distributions correspond to projections of translated hyperspheres and translated normal distributions, respectively; are identical in their central distributions when their noncentrality parameters are zero; and converge in high-dimensional spaces, but diverge in low-dimension spaces and for large noncentrality parameters. See Figure 1. These properties ultimately stem from the counterintuitive properties of the solid hypersphere which concentrates its volume on a thin ultraspherical shell in high-dimensional spaces [3], allowing one to use the ultraspherical distributions as surrogates for the classical noncentral t and F distributions in high-dimensional spaces.

4. Entropic Convex Duals Expansion in Terms of Orthogonal Polynomials

Bayes–Jaynes–Gibbs data-constrained maximal entropy priors—simply designated as Gibbs priors in the following—can be objectively computed for dense datasets. See Le Blanc [8] for an extensive review on the subject, from which we recall that the Bayes–Laplace prior and posterior update rules are rooted in the convex geometry of Shannon’s entropy function, with the Kullback–Leibler relative entropy being a Bregman divergence defined in terms of the former. The Gibbs priors can be formally expressed as

$$\pi \left(r_o\right)=\frac{1}{Z\left(\lambda \right)}exp\left({\int }_{\mathcal{R}}\lambda \left(r\right)\upsilon \left(r\right|r_o)\mathrm{d}r\right),$$

(45)

where the partition function $Z\left(\lambda \right)$ is given by

$$Z\left(\lambda \right)={\int }_{{\mathcal{R}}_o}exp\left({\int }_{\mathcal{R}}\lambda \left(r\right)\upsilon \left(r\right|r_o)\mathrm{d}r\right)\mathrm{d}{r}_o,$$

(46)

and where $\lambda \left(r\right)$ is the entropic convex dual of the empirical density $\upsilon \left(r\right)$. The latter is obtained through the unconstrained minimization of the Gibbs potential

$$\underset{\lambda }{inf}{G}_{\rho }\left(\lambda \right)=\underset{\lambda }{inf}\left[logZ\left(\lambda \right)-{\int }_{\mathcal{R}}\lambda \left(r\right)\upsilon \left(r\right)\mathrm{d}r\right],$$

(47)

which, by convex duality, corresponds to Jaynes data-constrained maximal entropy. The corresponding Gibbs–Jaynes model for a generic empirical density $\upsilon \left(r\right)$ is given by

$$\upsilon \left(r\right)={\int }_{{\mathcal{R}}_o}\upsilon \left(r\right|r_o)\pi \left(r_o\right)\mathrm{d}{r}_o=\upsilon \left(r\right|0){\int }_{{\mathcal{R}}_o}T\left(r\right|r_o)\pi \left(r_o\right)\mathrm{d}{r}_o=\upsilon \left(r\right|0)BF\left(r\right).$$

(48)

As posed, solving for the entropic convex dual function $\lambda \left(r\right)$ in (47) requires one to compute its value across the entire support domain of the empirical distribution $\upsilon \left(r\right)$, a task which can be expensive in computing terms. Now, since the ultraspherical noncentral distribution translation factors $T\left(r\right|r_o)$ are generating functions for orthogonal polynomial families, and since the corresponding central distributions $\upsilon \left(r\right|0)$ are the family-defining weight functions $w\left(r\right)\equiv \upsilon \left(r\right|0)$, one can rewrite the unconstrained minimization problem in simpler terms. Indeed, the exponentiated term in the partition function (46) can be rewritten as

$$\begin{array}{ccc}\hfill {\int }_{\mathcal{R}}\lambda \left(r\right)\upsilon \left(r\right|r_o)\mathrm{d}r& =& {\int }_{\mathcal{R}}\lambda \left(r\right)T\left(r\right|r_o)w\left(r\right)\mathrm{d}r={\int }_{\mathcal{R}}\lambda \left(r\right)\left[\sum _{n=0}^{\infty }{c}_n{P}_n\left(r\right)r_o^n\right]w\left(r\right)\mathrm{d}r\hfill \\ & =& \sum _{n=0}^{\infty }{c}_n\left[{\int }_{\mathcal{R}}\lambda \left(r\right)P_n\left(r\right)w\left(r\right)\mathrm{d}r\right]r_o^n=\sum _{n=0}^{\infty }{c}_n{\lambda }_n{\parallel P_n\parallel }^2{r}_o^n\hfill \\ & =& \sum _{n=0}^{\infty }{\tilde{\lambda }}_n{r}_o^n\mathrm{with}{\tilde{\lambda }}_n=c_n{\lambda }_n{\parallel P_n\parallel }^2,\hfill \end{array}$$

(49)

where the polynomials $P_n$, with their respective multiplicative coefficients $c_n$, are listed in (44), and where the ${\lambda }_n$ are the entropic convex dual expansion coefficients on family-wise orthogonal polynomials that remain to be determined. Similarly, the additive constraint term in (47) can be rewritten as

$${\int }_{\mathcal{R}}\lambda \left(r\right)\upsilon \left(r\right)\mathrm{d}r=\sum _{n=0}^{\infty }{\lambda }_n{\int }_{\mathcal{R}}{P}_n\left(r\right)\upsilon \left(r\right)\mathrm{d}r=\sum _{n=0}^{\infty }{\lambda }_n{\int }_{\mathcal{R}}{P}_n\left(r\right)BF\left(r\right)w\left(r\right)\mathrm{d}r=\sum _{n=0}^{\infty }{\tilde{\lambda }}_n{\tilde{\beta }}_n,$$

(50)

where $\upsilon \left(r\right)$ is the empirical density to be modeled, and where

$${\tilde{\beta }}_n=\frac{1}{c_n{\parallel P_n\parallel }^2}{\int }_{\mathcal{R}}{P}_n\left(r\right)BF\left(r\right)w\left(r\right)\mathrm{d}r,$$

(51)

with, in particular, ${\tilde{\beta }}_0=1$. Within the division by the factor $c_n$, the latter equality is the Bayes factor $BF\left(r\right)$ expansion on the orthogonal polynomials in random sample space. The unconstrained optimization problem can thus be reformulated as

$$\underset{\left\{{\tilde{\lambda }}_n\right\}}{inf}\left[log\left({\int }_{{\mathcal{R}}_o}exp(\sum _{n=0}^{\infty }{\tilde{\lambda }}_n{r}_o^n)\mathrm{d}{r}_o\right)-\sum _{n=0}^{\infty }{\tilde{\lambda }}_n{\tilde{\beta }}_n\right]$$

(52)

in terms of orthogonal polynomial expansion coefficient sets ${\left\{{\tilde{\lambda }}_n\right\}}_{n=0}^{\infty }$ for the continuous entropic convex dual function $\lambda \left(r\right)$, sets which can be restricted to a small finite number of coefficients, as can be assessed by the Kullback–Leibler divergence between the empirical density and its model in terms of orthogonal polynomials. At the minimum of the Gibbs potential (52), one has the simple condition

$$\frac{{\int }_{{\mathcal{R}}_o}exp\left({\sum }_n{\tilde{\lambda }}_n{r}_o^n\right)r_o^n\mathrm{d}{r}_o}{{\int }_{{\mathcal{R}}_o}exp\left({\sum }_n{\tilde{\lambda }}_n{r}_o^n\right)\mathrm{d}{r}_o}={\int }_{{\mathcal{R}}_o}\pi \left(r_o\right)r_o^n\mathrm{d}{r}_o={\tilde{\beta }}_n,$$

(53)

which states that the $n^{\mathrm{th}}$ moment of the noncentrality parameter $r_o$, as weighted by the Gibbs prior

$$\pi \left(r_o\right)=\frac{exp\left({\sum }_n{\tilde{\lambda }}_n{r}_o^n\right)}{{\int }_{{\mathcal{R}}_o}exp\left({\sum }_n{\tilde{\lambda }}_n{r}_o^n\right)\mathrm{d}{r}_o},$$

(54)

is, within the factor $c_n$, equal to the $n^{\mathrm{th}}$ coefficient of the Bayes factor expansion (51) on the orthogonal polynomial basis in random sample space. In practical terms, the determination of the Gibbs prior in terms of a small number of polynomial expansion coefficients ${\tilde{\lambda }}_n$ for the entropic dual convex results in a substantial reduction in the computing time needed to find the Gibbs potential minimum.

Once obtained, the Gibbs prior can be used to modelize the empirical density $\upsilon \left(r\right)$ according to Equation (48), in which the central distribution $\upsilon \left(r\right|0)$ stands for the weight function $w\left(r\right)$ for the ultraspherical distributions listed in (42). A generic density $\upsilon \left(r\right)$ can thus be expanded as

$$\begin{array}{ccc}\hfill \upsilon \left(r\right)& =& {\int }_{{\mathcal{R}}_o}\upsilon \left(r\right|r_o)\pi \left(r_o\right)\mathrm{d}{r}_o=w\left(r\right){\int }_{{\mathcal{R}}_o}T\left(r\right|r_o)\pi \left(r_o\right)\mathrm{d}{r}_o\hfill \\ & =& w\left(r\right){\int }_{{\mathcal{R}}_o}\left[\sum _{n=0}^{\infty }{c}_n{P}_n\left(r\right)r_o^n\right]\pi \left(r_o\right)\mathrm{d}{r}_o=w\left(r\right)\sum _{n=0}^{\infty }{c}_n\left[{\int }_{{\mathcal{R}}_o}{r}_o^n\pi \left(r_o\right)\mathrm{d}{r}_o\right]P_n\left(r\right)\hfill \\ & =& w\left(r\right)\sum _{n=0}^{\infty }{c}_n{\tilde{\beta }}_n{P}_n\left(r\right)=w\left(r\right)\sum _{n=0}^{\infty }\left[\frac{1}{\parallel P_n{\parallel }^2}{\int }_{\mathcal{R}}{P}_n\left(r^{\prime }\right)BF\left(r^{\prime }\right)w\left(r^{\prime }\right)\mathrm{d}{r}^{\prime }\right]P_n\left(r\right)\hfill \end{array}$$

(55)

where we have invoked the optimality condition (53) in going from the second to the third line. This sequence of equalities provides us with two alternative ways,

$$\begin{array}{}\mathrm{(56a)}& BF\left(r\right)=\frac{\upsilon \left(r\right)}{w\left(r\right)}& ={\int }_{{\mathcal{R}}_o}T\left(r\right|r_o)\pi \left(r_o\right)\mathrm{d}{r}_o\hfill \mathrm{(56b)}& & =\sum _{n=0}^{\infty }{c}_n\left[{\int }_{{\mathcal{R}}_o}{r}_o^n\pi \left(r_o\right)\mathrm{d}{r}_o\right]P_n\left(r\right)=\sum _{n=0}^{\infty }\left[\frac{1}{\parallel P_n{\parallel }^2}{\int }_{\mathcal{R}}{P}_n\left(r^{\prime }\right)BF\left(r^{\prime }\right)w\left(r^{\prime }\right)\mathrm{d}{r}^{\prime }\right]P_n\left(r\right),\end{array}$$

to model the Bayes factors and the empirical densities: either one determines the Gibbs prior $\pi \left(r_o\right)$ via the optimization problem (52) and uses it to weigh the analytical orthogonal polynomial-generating functions $T\left(r\right|r_o)$ listed in (44); or, one directly determines the prior moments in parametric space in terms of the Bayes factor $BF\left(r\right)$ expansion coefficients in random variable space, as per the last equality above. The latter strategy amounts to the construction of a probability distribution in terms of its moments [12]. In that respect, recall the Hausdorff moment problem, which stipulates that the collection of all moments of a probability distribution on a bounded interval uniquely determine the distribution [19], and the Hamburger moment problem, which considers the uniqueness of solutions for the same problem on the unbounded real line [20].

Numerical exploration indicates that the convergence of the polynomial expansion for the ultraspherical noncentral t and F distributions’ generating functions is affected by the Gibbs phenomenon at the extremities of their finite noncentrality parameter ranges [21]. For this reason, when considering the noncentral t and F distributions, one can either solve for the Gibbs prior to weigh the corresponding analytic generating functions (44) in constructing the models, or use the large b limit for the distributions which converge with the noncentral normal N and ${\chi }^2$ distributions, respectively, together with the associated Hermite and Laguerre orthogonal polynomial families, to directly compute the prior moments in parametric space in terms of the Bayes factor expansion coefficients in random variable space, as according to (56). Both strategies are found to be unaffected by the Gibbs effect, as exemplified in the next section.

5. Genomics Examples

We first revisit, in this section, the NCBI Gene Expression Omnibus head and neck squamous cell carcinoma microarray dataset produced by [22], pertaining to 22 paired samples ($N=44,$${\nu }_2=42,$$b=41/2$) of normal versus cancerous tissues, and interrogating 11,302 genes via 2-sample t-tests. In Le Blanc [10], a maximal entropy prior equal to the central ultraspherical t distribution (1) was postulated. In Le Blanc [8], the entropic dual convex $\lambda \left(x\right)$ was computed on the full random variable range $-1<x<1$. Here, the Bayes factor is provided with two models. The first model $\upsilon \left(x\right)=w_t\left(x\right)\int T\left(x\right|z)\pi \left(z\right)\mathrm{d}z$ is given in terms of the Gibbs prior $\pi \left(z\right)$ (54) and the analytical Gegenbauer polynomial-generating function $T\left(x\right|z)$ (11). We assessed the computational gains in our determination of the entropic convex dual $\lambda \left(x\right)$ via Equation (52): with the empirical distribution $\upsilon \left(x\right)$ binned in 200 bins, and using the $\mathrm{MATLAB}$^® optimization subroutine fminunc in a 64-bit Windows environment on a PC with an Intel(R) Core(TM) i9-9900K CPU @ 3.60 GHz processor and 64 GB of RAM, the determination times for the entropic convex dual $\lambda \left(x\right)$ in terms of a small number of Gegenbauer polynomial coefficients were curtailed by a factor ranging from 100 to 1000, compared to previous determination times of $\lambda \left(x\right)$ on the full random variable range $-1<x<1$, with elapsed times as short as a few tenths of a second. In the second model, according to (56) and with the Hermite polynomials standing for the Gegenbauer polynomials in the large b limit, the Bayes factor $BF\left(r\right)$ expansion in random variable space on the Hermite polynomials provides us with the Gibbs prior moments $\int \pi \left(\delta \right){\delta }^n\mathrm{d}\delta$ in parametric space. As can be seen in Figure 2, the convergence of the two models with the empirical distribution is rapidly achieved with ten coefficients or less.

Next, we study a genome-wide association study (GWAS) dataset. A GWAS is an observational study assessing a genome-wide set of genetic variants in different individuals, and seeking to identify statistically significant variant-trait associations. Such studies commonly focus on associations between single-nucleotide polymorphisms (SNPs) and traits. We retrieved the GenoMICC EUR vs. UK biobank controls dataset from the GenOMICC (Genetics Of Mortality In Critical Care) GWAS, comparing 2244 critically ill patients with COVID-19 from UK intensive care units with European ancestry-matched control individuals selected from the large population-based cohort of the UK Biobank [23]. A logistic regression model was used for each of the 4,380,209 SNPs individually tested for statistical significance. We computed the empirical density of all the statistical test-associated ${\chi }_{{\nu }_1=1}^2$ r-statistics, and used it to compute the first eight terms of the direct polynomial expansion (56) on Laguerre polynomials for the Bayes Factor $BF\left(p\right(r\left)\right)$. As illustrated in Figure 3, the accrual of the successive polynomial expansion terms allows for an incrementally better fit of the p-value empirical density, which strongly deviates from the NHST null hypothesis $U(0,1)$ in the low p-value range where associations are detected. In the Bayesian framework, NHST statistical significance is replaced by the strength of Bayesian evidence, as assessed by magnitude of $BF\left(p\right)$, which, in turn, allows for the computation of a local false discovery rate [11]

$$\mathrm{fdr}\left(p\right)=1/(1+BF(p\left)\right)$$

(57)

as is also illustrated in Figure 3.

6. Geophysics Examples

We begin this section by modeling, in Figure 4, Earth’s above-sea-level emerging land/ice latitudinal density [24]. Because an ideal Earth surface can be described by the sphere ${\mathcal{S}}^{{\nu }_2=2}$, we have that the parameter $b=({\nu }_2-1)/2=1/2$, with the weight factor $w_t^{(1/2)}\left(x\right)$ in (15) reducing to the constant $1/2.$ As a consequence, Equation (55) is equivalent to a standard expansion on Legendre polynomials—equivalently, Gegenbauer $C_n^{(1/2)}$ polynomials—except for the fact that our convention regarding the central distribution normalization $w_t^{(1/2)}\left(x\right)=1/2$ (rather than the usual weight $w\left(x\right)=1$) on the span $[-1,1]$ ensures that expansion (55) defines a normalized density.

We conclude with an example in which the translating factor (9) for the ultraspherical noncentral t distribution is kept as-is, instead of using its expansion (11) in terms of Gegenbauer polynomials, to argue that the noncentral spherical t distribution (8) on ${\mathcal{S}}^2$ can be used to describe the geometry of the gravity forces of a simple tide model. Consider Figure 5, which describes the gravitational pull of an ideal Moon of mass m on the thin water layer covering an ideal Earth of mass M and radius R at distance D from the Moon. The gravitational tidal force $\mathcal{T}$ per unit mass on point s on the Earth’s surface is given by

$$\mathcal{T}=\frac{Gm}{{(D-Rcos\theta )}^2+R^2{sin}^2\theta }=\frac{Gm}{D^2-2DRcos\theta +R^2}.$$

(58)

The horizontal and vertical components of that force are given by

$$\begin{array}{ccc}\hfill {\mathcal{T}}_h& =& +\frac{Gm}{D^2-2DRcos\theta +R^2}\times cos\phi =+Gm\frac{D-Rcos\theta }{{[D^2-2DRcos\theta +R^2]}^{3/2}},\hfill \end{array}$$

(59)

$$\begin{array}{ccc}\hfill {\mathcal{T}}_v& =& -\frac{Gm}{D^2-2DRcos\theta +R^2}\times sin\phi =-Gm\frac{Rsin\theta }{{[D^2-2DRcos\theta +R^2]}^{3/2}},\hfill \end{array}$$

(60)

respectively. Setting our physical units so that $Gm/D^2=1$, and defining $cos{\theta }_{\delta }=R/D$, we have that

$$\begin{array}{ccc}\hfill {\mathcal{T}}_h& =& \frac{1-cos\theta cos{\theta }_{\delta }}{{[1-2cos\theta cos{\theta }_{\delta }+{cos}^2{\theta }_{\delta }]}^{3/2}}=T_t^{({\nu }_2=2)}\left(\theta \right|\delta ),\hfill \end{array}$$

(61)

$$\begin{array}{ccc}\hfill {\mathcal{T}}_v& =& \frac{-sin\theta }{{[1-2cos\theta cos{\theta }_{\delta }+{cos}^2{\theta }_{\delta }]}^{3/2}}.\hfill \end{array}$$

(62)

The horizontal tidal force component ${\mathcal{T}}_h$ is thus simply given by the noncentral spherical t distribution translation factor $T_t^{({\nu }_2=2)}\left(\theta \right|\delta )$, as provided by Equation (4). One can integrate this horizontal tidal force over the entire spherical shell by weighing it with the central spherical distribution ${\upsilon }_t^{({\nu }_2=2)}\left(\theta \right|0)$ to account for the spherical geometry of the Earth. Since the integrand to this integral corresponds to the noncentral spherical t distribution ${\upsilon }_t^{({\nu }_2=2)}\left(\theta \right|\delta )=T_t^{({\nu }_2=2)}\left(\theta \right|\delta ){\upsilon }_t^{({\nu }_2=2)}\left(\theta \right|0)$ on ${\mathcal{S}}^2$, the tidal force integrates to one in our unit system. In order for the system to be stationary, an opposing force must be opposed to this integrated tidal force of one. The inertial centrifugal force originating from the joint rotation of the Moon and Earth around their center of mass plays such a role. It is evaluated to be given by $-Gm/D^2=-1$ in our unit system [25], which we add to the normalized horizontal tidal force ${\mathcal{T}}_h$ above. We have plotted the resulting force field in Figure 5, with the noncentrality parameter $cos{\theta }_{\delta }$ set to an unrealistic value of 0.2 to enhance the visualization of the geometrical distribution of the tidal forces. Unlike most illustrative tide diagrams with symmetrical bulges found in the literature, our tide diagram, together with its unsymmetrical equatorial water bulges, provides a more accurate visual depiction of the expected asymmetry of the tidal forces resulting from purely geometrical considerations.

7. Discussion

We have seen that the univariate noncentral distributions can be constructed in a modular fashion by multiplying their central distributions with specific translation factors. Using geometrical arguments, we have found that the translation factors for the ultraspherical noncentral t, normal, F, and ${\chi }^2$ distributions stand for generating functions for the Gegenbauer, Hermite, Jacobi, and Laguerre polynomial families, respectively, with their central distributions standing for the corresponding polynomial family-defining weights. These developments clearly link four of the most important classical continuous probability distributions with the powerful orthogonal polynomial formalism. To the best of our knowledge, the derivation of these translation factors and their identification as orthogonal polynomial family-generating functions has not been carried out before.

Jaynes’ maximal entropy prior is obtained through the unconstrained minimization of the Gibbs potential. In parametric Bayesian inference, the formal expression for the Gibbs potential comprises an integral of the product of the empirical distribution’s entropic convex dual $\lambda \left(r\right)$ times the parametric kernel, or, equivalently, the likelihood function $\rho \left(r\right|r_o).$ In the case of the ultraspherical noncentral distributions, the latter integral yields discretized expansion coefficients of the entropic convex duals on orthogonal polynomial bases. The determination of the entropic convex duals is thus reduced to the much simpler and computationally economical determination of a few low-order orthogonal polynomial coefficients. By invoking the moment problem and the duality principle, prior moments in parametric space are equated with Bayes factors expansion coefficients over orthogonal polynomial bases in random variable space. To the best of our knowledge, the expansion and discretization of convex duals over orthogonal polynomial bases has not been proposed before. In an approach which bears some similarities to our work, Alibrandi and Ricciardi [26] proposed a moment-based approach, the use of a discretized kernel set ${\left\{\rho \left(r\right|r_i)\right\}}_{i=1}^N$, and the use of Jaynes’ maximal entropy principle to determine the set’s weighting factors $\left\{{\pi }_i\right|{\sum }_{i=1}^N{\pi }_i=1\}$. While their kernel discretization procedure is an ad hoc procedure, our kernel discretization is principled since it is based on the powerful orthogonal polynomial formalism.

The machine learning community has begun to exploit the classical orthogonal polynomial formalism. A non-exhaustive review identifies the use of orthogonal polynomials in optical character recognition [27], in support vector machine kernel construction [28], and in polynomial-based iteration methods for symmetric linear systems [29,30]. In a similar vein, we hope to contribute the present formalism to both the statistical and the machine learning communities.