SOCR ≫ | TCIU Website ≫ | TCIU GitHub ≫ |
The domain of the kime variable (\(k\)) is the complex plane parameterized by pairs of Descartes Cartesian coordinates, conjugate-pairs coordinates, or polar coordinates:
\[\mathbb{C \cong}\mathbb{R}^{2} = \left\{ z = (\ x\ ,y\ )\mid\ x,\ y\mathbb{\in R\ } \right\} \cong\] \[\left\{ \left( \ z\ ,\overline{z}\ \right)\mid\ z\mathbb{\in C,\ }z = x + iy,\ \ \overline{z}\ = x - iy,\ \ x,y\mathbb{\in R\ } \right\} \cong\]
\[\{\ k = r\ e^{i\varphi} = r\left( \cos\varphi + i\sin\varphi \right)\ \mid \ r \geq 0,\ - \pi \leq \varphi < \pi\ \}.\]
The Wirtinger derivative of a continuously differentiable function (\(f\)) of kime (\(k\)), \(f(k)\), and its conjugate are defined as first-order linear partial differential operators:
In Cartesian coordinates:
\(f'(z) = \frac{\partial f(z)}{\partial z} = \frac{1}{2}\left( \frac{\partial f}{\partial x} - i\frac{\partial f}{\partial y} \right)\) and \(f'\left( \overline{z} \right) = \frac{\partial f\left( \overline{z} \right)}{\partial\overline{z}} = \frac{1}{2}\left( \frac{\partial f}{\partial x} + i\frac{\partial f}{\partial y} \right).\)
In Conjugate-pair basis: \(df = \partial f + \overline{\partial}f = \frac{\partial f}{\partial z}dz + \frac{\partial f}{\partial\overline{z}}d\overline{z}.\)
In Polar kime coordinates:
\[f'(k) = \frac{\partial f(k)}{\partial k} = \frac{1}{2}\left( \cos\varphi\frac{\partial f}{\partial r} - \frac{1}{r}\sin\varphi\frac{\partial f}{\partial\varphi} - i\left( \sin\varphi\frac{\partial f}{\partial r} + \frac{1}{r}\cos\varphi\frac{\partial f}{\partial\varphi} \right) \right) = \frac{e^{- i\varphi}}{2}\left( \frac{\partial f}{\partial r} - \frac{i}{r}\ \frac{\partial f}{\partial\varphi} \right)\]
and
\[f'\left( \overline{k} \right) = \frac{\partial f\left( \overline{k} \right)}{\partial\overline{k}} = \frac{1}{2}\left( \cos\varphi\frac{\partial f}{\partial r} - \frac{1}{r}\sin\varphi\frac{\partial f}{\partial\varphi} + i\left( \sin\varphi\frac{\partial f}{\partial r} + \frac{1}{r}\cos\varphi\frac{\partial f}{\partial\varphi} \right) \right) = \frac{e^{i\varphi}}{2}\left( \frac{\partial f}{\partial r} + \frac{i}{r}\ \frac{\partial f}{\partial\varphi} \right).\]
Notes:
$|
\[\begin{matrix} x = r\ \cos\varphi \\ y = r\ \sin\varphi \\ \end{matrix}\].$ , $|
\[\begin{matrix} r = \sqrt{x^{2} + y^{2}} \\ \varphi = \arctan\left( \frac{y}{x} \right) = \arctan{(y,x)} \\ \end{matrix}\].$, $|
\[\begin{matrix} \frac{\partial}{\partial x} = \cos\varphi\frac{\partial}{\partial r} - \frac{1}{r}\ \sin\varphi\frac{\partial}{\partial\varphi} \\ \frac{\partial}{\partial y} = \sin\varphi\frac{\partial}{\partial r} + \frac{1}{r}\ \cos\varphi\frac{\partial}{\partial\varphi} \\ \end{matrix}\].$ , see Korn GA, Korn TM. Mathematical handbook for scientists and engineers: definitions, theorems, and formulas for reference and review, Courier Corporation; 2000.
\[(x,y) \rightarrow \left( \frac{1}{2}\left( z + \overline{z} \right),\ \ \ \frac{1}{2i}\left( z - \overline{z} \right) \right),\]
\[\frac{\partial}{\partial z} = \frac{\partial}{\partial x}\frac{\partial x}{\partial z} + \frac{\partial}{\partial y}\frac{\partial y}{\partial z}.\]
Therefore, \(\frac{\partial x}{\partial z} = \frac{1}{2}\frac{\partial\left( z + \overline{z} \right)}{\partial z} = \frac{1}{2}\) and \(\frac{\partial y}{\partial z} = \frac{1}{2i}\frac{\partial\left( z - \overline{z} \right)}{\partial z} = \frac{1}{2i} = - \frac{i}{2}.\)
Similarly,
\(\frac{\partial x}{\partial\overline{z}} = \frac{1}{2}\frac{\partial\left( z + \overline{z} \right)}{\partial\overline{z}} = \frac{1}{2}\) and \(\frac{\partial y}{\partial\overline{z}} = \frac{1}{2i}\frac{\partial\left( z - \overline{z} \right)}{\partial\overline{z}} = - \frac{1}{2i} = \frac{i}{2}\).
This explains the Cartesian coordinate derivatives:
\(f'(z) = \frac{\partial f(z)}{\partial z} = \frac{1}{2}\left( \frac{\partial f}{\partial x} - i\frac{\partial f}{\partial y} \right)\) and \(f'\left( \overline{z} \right) = \frac{\partial f\left( \overline{z} \right)}{\partial\overline{z}} = \frac{1}{2}\left( \frac{\partial f}{\partial x} + i\frac{\partial f}{\partial y} \right)\).
Below, we present the core principles of Wirtinger differentiation} and integration:
\[\left| \begin{matrix} x = \frac{1}{2}\left( z + \overline{z} \right) \\ y = \frac{1}{2i}\left( z - \overline{z} \right) \\ \end{matrix}\ . \right.\ \]
Wirtinger differentiation: The Wirtinger derivative of \(f\), \(df_{z}\) is an \(\mathbb{R}\)-linear operator on the tangent space \(T_{z}\mathbb{C \cong C}\), i.e., \(df_{z}\) is a differential 1-form on \(\mathbb{C}\). However, any such \(\mathbb{R}\)-linear operator (\(A\)) on \(C\) can be uniquely decomposed} as \(A = B + C\), where \(B\) is its complex-linear part* (i.e., \(B(iz) = iBz\)), and \(C\) is its complex-antilinear part (i.e., \(C(iz) = - iCz\)). The reverse (composition) mapping is \(Bz = \frac{1}{2}\left( Az - iA(iz) \right)\) and \(Cz = \frac{1}{2}\left( Az + iA(iz) \right)\).
For the Wirtinger derivative, this duality of the decomposition of \(\mathbb{R}\)-linear operators characterizes the conjugate partial differential operators \(\partial\) and \(\overline{\partial}\). That is, for all differentiable complex functions \(f\mathbb{:C \rightarrow C}\), the derivative can be uniquely decomposed as \(\mathbf{d}\mathbf{f}_{\mathbf{z}}\mathbf{= \partial f +}\overline{\mathbf{\partial}}\mathbf{f}\), where \(\partial\) is its complex-linear part} (\(\partial iz = i\partial z\)), and \(\overline{\partial}\) is its complex-antilinear part* (\(\overline{\partial}iz = - i\overline{\partial}z\)).
Applying the operators \(\frac{\mathbf{\partial}}{\mathbf{\partial z}}\) and \(\frac{\mathbf{\partial}}{\mathbf{\partial}\overline{\mathbf{z}}}\) to the identify function (\(z \rightarrow z = x + iy\)) and its complex-conjugate (\(z \rightarrow \overline{z} = x - iy\)) yields the natural derivatives: \(dz = dx + idy\) and \(d\overline{z} = dx - idy\). For each point in \(\mathbb{C}\), \(\{ dz,d\overline{z}\}\) represents a conjugate-pair basis for the \(\mathbb{C}\) cotangent space, with a dual basis of the partial differential operators:
\[\left\{ \frac{\partial}{\partial z}\ ,\ \frac{\partial}{\partial\overline{z}} \right\}\ .\]
Thus, for any smooth complex functions \(f\mathbb{:C \rightarrow C}\),
\[df = \partial f + \overline{\partial}f = \frac{\partial f}{\partial z}dz + \frac{\partial f}{\partial\overline{z}}d\overline{z}\ .\]
\[\lim_{\left| z_{i + 1} - z_{i} \right| \rightarrow 0}{\sum_{i = 1}^{n - 1}\left( f(z_{i})(z_{i + 1} - z_{i}) \right)} \cong \oint_{z_{a}}^{z_{b}}{f\left( z_{i} \right)dz}\ .\]
This assumes the path is a polygonal arc joining \(z_{a}\) to \(z_{b}\), via \(z_{1} = z_{a},\ z_{2},\ z_{3},\ ,\ z_{n} = z_{b}\), and we integrate the piecewise constant function \(f(z_{i})\) on the arc joining \(z_{i} \rightarrow z_{i + 1}\). Clearly the path \(z_{a} \rightarrow z_{b}\) needs to be defined and the limit of the generalized Riemann sums, as \(n \rightarrow \infty\), will yield a complex number representing the Wirtinger integral of the function over the path. Similarly, we can extend the classical area integrals, indefinite integral, and Laplacian:
Definite area integral: for \(\Omega \subseteq \mathbb{C}\), \(\int_{\Omega}^{}{f(z)dzd\overline{z}}\).
Indefinite integral: \(\int_{}^{}{f(z)dzd\overline{z}}\), \(df = \frac{\partial f}{\partial z}dz + \frac{\partial f}{\partial\overline{z}}d\overline{z}\),
The Laplacian in terms of conjugate pair coordinates is \(\nabla^{2}f \equiv \mathrm{\Delta}f = 4\frac{\partial}{dz}\frac{\partial f}{d\overline{z}} = 4\frac{\partial}{d\overline{z}}\frac{\partial f}{dz}\ .\)
More details about Wirtinger calculus of differentiation and integration are provided later.
In general, the classical concepts of derivative and rate of change are well defined for a kime-process \(f(t,\varphi)\) with respect to the (positive) real variable time (\(t \equiv r = |\kappa| = |r\ e^{i\varphi}|\)). Depending on the context, we interchangeably use \(\varphi, \phi, \psi, \theta\) and other lowercase Greek letters to represent the kime-phases. The smooth time rate of change of the process \(f\) is explicated as the classical partial derivative, \(\frac{\partial f}{\partial r}\).
Above we showed that the partial derivatives of a kime-process in Cartesian kime coordinates are also well defined, using Wirtinger differentiation. However, in the polar representation, this strategy may not apply for the second variable (kime-phase), \(\varphi\), which could be interpreted as a distribution. In essence, \(\varphi\) represents an unobservable random sampling variable from a symmetric distribution compactly-supported on \(\lbrack - \pi,\pi)\), or a periodic distribution.
Above we showed that in polar kime coordinates, analytic functions can be differentiated as follows
\[f'(k) = \frac{\partial f(k)}{\partial k} = \frac{1}{2}\left( \cos\varphi\frac{\partial f}{\partial r} - \frac{1}{r}\sin\varphi\frac{\partial f}{\partial\varphi} - i\left( \sin\varphi\frac{\partial f}{\partial r} + \frac{1}{r}\cos\varphi\frac{\partial f}{\partial\varphi} \right) \right) = \frac{e^{- i\varphi}}{2}\left( \frac{\partial f}{\partial r} - \frac{i}{r}\ \frac{\partial f}{\partial\varphi} \right)\]
and
\[f'\left( \overline{k} \right) = \frac{\partial f\left( \overline{k} \right)}{\partial\overline{k}} = \frac{1}{2}\left( \cos\varphi\frac{\partial f}{\partial r} - \frac{1}{r}\sin\varphi\frac{\partial f}{\partial\varphi} + i\left( \sin\varphi\frac{\partial f}{\partial r} + \frac{1}{r}\cos\varphi\frac{\partial f}{\partial\varphi} \right) \right) = \frac{e^{i\varphi}}{2}\left( \frac{\partial f}{\partial r} + \frac{i}{r}\ \frac{\partial f}{\partial\varphi} \right).\]
The problem with computing and interpreting \(f'(k)\), and in particular \(\frac{\partial f}{\partial\varphi}\), is that \(f(r,\varphi)\) may not be analytic. This is especially interesting since one plausible spacekime interpretation of the enigmatic kime-phase is that \(\varphi\) indexes the intrinsic stochastic sampling from some symmetric distribution, which corresponds to repeated sampling, or multiple IID random observations. This suggests that \(\varphi\) does not vary smoothly, but is rather a chaotic quantity (random variable) following some symmetric distribution \(\Phi_{}\), supported over \(\lbrack - \pi,\pi)\). In other words, the *kime-phase change} quantity \(\partial\varphi\) may be stochastic! Hence, we need to investigate approaches to:
Note that the differentiability of \(f(k)\) is not in question for the Cartesian representation of the kime domain, as Wirtinger derivatives explicate the multivalued differentiation of a multi-variable function, e.g., \(f'(k),\ \forall\kappa \equiv (x,y\mathbb{) \in C}\). Wirtinger differentiation provides a robust formulation even in the polar coordinate representation of the complex domain.
The problem arises in the spacekime-interpretation where we try to tie the analytical math formalism with the computational algorithms and the quantum-physics/statistical/experimental interpretation of the kime-phase as an intrinsically stochastic repeated sampling variability.
The following well-known lemma suggests that we can quantify the *distribution of the kime-phase changes} in the case when we merge the math formalism, computational modeling, and the statistical/experimental interpretation of the stochastic kime-phases.
We’ll use the standard convention where capital letters (\(\Gamma,\Phi,\Theta\)) denote random variables} and their lower-case counterparts (\(\gamma,\phi,\theta\)) denote the corresponding values} of these quantities. Let’s denote by \(P \equiv \Pr\) the probability functions with proper subindices indicating the specific corresponding random process.
Lemma: Suppose \(\phi,\ \theta \sim F_{\mathcal{D}}\) are IIDs associated with some distribution \(\mathcal{D}\), with CDF \(F_{\mathcal{D}}\) and PDF \(f_{\mathcal{D}}\). Then, the distribution of their difference \(\gamma = \phi - \theta \sim F_{\Gamma}\) is effectively the autocorrelation function
\[f_{\Gamma}(\gamma) = \ \ \int_{- \infty}^{\infty}{f_{\mathcal{D}}(\phi)f_{\mathcal{D}}(\phi - \gamma)d\phi}\ .\]
Proof: This lemma quantifies the phase-change distribution, not the actual phase-distribution. We’ll start with the most general case of two independent random variables, \(\Phi\bot\Theta\), following potentially different distributions, \(\Phi \sim F_{\Phi}\), \(\Theta \sim F_{\Theta}\).
Let’s compute the cumulative distribution of the difference random variable \(\Gamma = \Phi - \Theta\). In practice, \(\Gamma\) corresponds to the instantaneous change in the kime-phase between two random repeated experiments. Note that due to the exchangeability property associated with IID sampling {See SOCR on Overleaf}, the distribution of \(\Gamma = \Phi - \Theta\) is invariant with respond to the sequential order at which we record the pair of random observations. More specifically, if the random sample of kime-phases is \(\{\phi_{1},\ \phi_{2},\ \phi_{3},\ \cdots,\ \phi_{n},\ \cdots\}\), then for any index permutation \(\pi( \cdot )\), the distributions of
\[\Gamma_{1} = \Phi_{1} - \Phi_{2},\ \ \Gamma_{1}' = \Phi_{\pi(1)} - \Phi_{\pi(2)}\]
are the same, i.e., \(\Gamma_{1},\ \Gamma_{1}' \sim F_{\Gamma}\). Hence, the phase-change distribution, \(F_{\Gamma}\), is symmetric.
Let’s compute the CDF of \(\Gamma\), the difference between any two kime-phases,
\[F_{\Gamma}(\gamma) \equiv P_{\Gamma}(\Gamma \leq \gamma) = P_{\Gamma}(\Phi - \Theta \leq \gamma) = P_{\Phi}(\Phi \leq \gamma + \Theta)\ .\]
Jointly, \((\phi,\ \theta) \sim f_{(\Phi,\Theta)}\) and
\[\underbrace{\ \ F_{\Gamma}(\gamma)\ }_{difference\ CDF} = P_{\Phi}(\Phi \leq \gamma + \Theta) = \int_{- \infty}^{\infty}{\int_{- \infty}^{\gamma + \theta} {\underbrace{f_{\Phi,\Theta}(\phi,\ \theta)}_{joint\ PDF}\ d\phi d\theta}}\ \ \Longrightarrow\]
\[\underbrace{f_{\Gamma}(\gamma)}_{difference\ PDF} = \frac{d}{d\gamma}F_{\Gamma}(\gamma) = \frac{d}{d\gamma}\int_{- \infty}^{\infty}{\int_{- \infty}^{\gamma + \theta} {\underbrace{f_{\Phi,\Theta}(\phi,\ \theta)}_{joint\ PDF}\ d\phi d\theta}}\ .\]
For a fixed \(\theta_{o}\), we’ll make a variable substitution \(\phi = \nu + \theta_{o}\), so that \(d\phi = d\nu\). Hence,
\[f_{\Gamma}(\gamma) = \int_{- \infty}^{\infty}{\left( \frac{d}{d\gamma}\int_{- \infty}^{\gamma}{f_{\Phi,\Theta}(\nu + \theta,\ \theta)\ d\nu} \right)d\theta}\ .\]
Integrals are functions of their (lower and upper) limits and for any function \(g(\phi,\theta_{o})\),
\[\frac{d}{d\gamma}\int_{l(\gamma)}^{h(\gamma)}{g\left( \phi,\ \theta_{o} \right)\ d\phi} = g\left( h(\gamma),\ \theta_{o} \right)h'(\gamma) - g\left( l(\gamma),\ \theta_{o} \right)l'(\gamma).\]
In particular, \(\frac{d}{d\gamma}\int_{- \infty}^{\gamma}{g\left( \phi,\ \theta_{o} \right)\ d\phi} = g\left( \gamma,\ \theta_{o} \right)\).
Therefore, for any \(\theta\),
\[f_{\Gamma}(\gamma) = \int_{- \infty}^{\infty}{\underbrace{\left( \frac{d}{d\gamma}\int_{- \infty}^{\gamma}{f_{\Phi,\Theta}(\nu + \theta,\ \theta)\ d\nu} \right)}_{\theta\ is\ fixed}\ d\theta} = \int_{- \infty}^{\infty}{f_{\Phi,\Theta}(\gamma + \theta,\ \theta)d\theta}\ .\]
This derivation of the density of the kime-phase difference is rather general, \(f_{\Gamma}(\gamma) = f_{\Phi,\Theta}(\gamma + \theta,\ \theta)d\theta\). It works for any bivariate process, whether or not \(\phi,\theta\) are independent, correlated, or follow the same distribution.
In the special case when \(\phi,\ \theta \sim F_{\mathcal{D}}\) are IIDs, we can apply another change of variables transformation, \(\gamma + \theta = \phi\), to get \(d\theta = d\phi\) and
\[f_{\Gamma}(\gamma) = \int_{- \infty}^{\infty}{f_{\Phi,\Theta}(\phi,\phi - \gamma)d\phi}\ \ \underbrace{\ \ = \ \ }_{\phi\bot\theta}\ \ \int_{- \infty}^{\infty}{f_{\Phi}(\phi)\ f_{\Theta}(\phi - \gamma)\ d\phi} \underbrace{\ \ = \ \ }_{IIDs}\ \ \int_{- \infty}^{\infty}{f_{\mathcal{D}}(\phi)\ f_{\mathcal{D}}(\phi - \gamma)\ d\phi}\ .\]
In the most general case, we would like to compute the measure (i.e., size) of any Borel subset of the support \(\lbrack - \pi,\pi)\) of the distribution \(\mathcal{D}\). This will help with answering questions (2) and (3) from the above list. Of course, for any finite or countable Borel set, \(B = \left\{ \phi_{i} \right\}_{1}^{\alpha}\mathcal{\sim D}\), it’s measure will be trivial, \(\mu(B) = 0\). However, many regular (e.g., contiguous) subsets \(B' \subseteq \lbrack - \pi,\pi)\) will be non-trivial, i.e., \(0 < f_{\mathcal{D}}\left( B' \right) \leq 1\).
Recall that the entire observable universe is finite, since the age of the universe is \(\sim 13.8\) billion years. Due to the continuously faster universal expansion, from any point in spacetime the radius and the full diameter of the entire visible universe are 46 billion and 93 billion light-years, respectively, see article 1 and article 2. Hence, the *amount of observable data is always finite}. In statistical data-science terms, the number of observations and sample-sizes are always finite. However, theoretically, we can model them as infinitely increasing, as in the first two laws of probability theory, the central limit theorem (CLT) and the law of large numbers (LLN).
Now, let’s assume that we have a sequence of randomly sampled kime-phases \(\left\{ \phi_{i} \right\}_{1}^{\alpha}\mathcal{\sim D}\), where \(\alpha < \infty\). In other words, assume we have observed a finite number of repeated measurements corresponding to multiple observations (corresponding to different kime-phase directions) acquired/recorded at the same time} under identical experimental conditions. The phase-dispersion* is the variability observed in the recorded measurements, e.g., \(\sigma_{(\phi - \theta)}\), is directly related to the distribution of the kime-phase differences, \(F_{\Gamma} = F_{(\Phi - \Theta)}\).
Note that changes in spacetime (jointly spacetime or separately space and time alone), always permit analytic calculations, since the presence of significant intrinsic spatiotemporal correlations induce smooth process dynamics over 4D Minkowski spacetime. However, this analyticity may be broken in the 5D spacekime, especially in the kime-phase dimension, since kime-phase indexing of observations in inherently stochastic, rather than smooth or analytic.
In our special bivariate case above, we can assume \(\phi \equiv \phi_{i},\ \ \theta \equiv \phi_{j} \sim F_{\mathcal{D}}\) are a pair of IIDs associated with an a priori kime-phase distribution \(\mathcal{D}\) and corresponding to identical experimental conditions. For simplicity, we are only considering a pair of fixed kime-phase indices \(1 \leq i < j\). Since both kime-phase distributions coincide, \(\Phi \equiv \Theta \equiv \mathcal{D}\), the above Lemma shows that the PDF \(f_{\Gamma}(\gamma)\) of their difference (\(\gamma = \phi - \ \theta \equiv \phi_{i} - \phi_{j}\)) is the autocorrelation function
\[f_{\Gamma}(\gamma) = \int_{- \infty}^{\infty}{f_{\Phi}(\phi)f_{\Theta}(\phi - \gamma)d\phi} = \int_{- \infty}^{\infty}{f_{\mathcal{D}}(\phi)f_{\mathcal{D}} \left( \underbrace{\phi - \gamma}_{\theta} \right)d\phi} = \mathbb{E}(\phi,\ \theta) = \underbrace{R_{phase}(\phi_{i},\phi_{j})}_{autocorrelation}\ .\]
What if the kime-phase distribution is *infinitely supported}? How can we interpret the phase distribution \(\Phi\) over the finite bounded interval \(- \pi \leq \phi < \pi\)?
Does it make sense to interpret the phase distribution in terms of spherical coordinates as helical, instead of circular distribution?
Can we account for both compactly supported and infinitely supported (univariate) distributions by applying renormalization adjusting for self-interaction feedback, or a regularization. Clearly there are many different (injective or bijective) transformations that map the domain of the reals to \(\lbrack - \pi,\pi)\). For instance,
\[T_{o}:( - \pi,\pi)\mathbb{\overbrace{\longrightarrow}^{bijective} R},\ \ T_{o}(x) = \frac{x}{\sqrt{\pi^{2} - x^{2}}},\]
\[T_{1}:\left\lbrack - \frac{\pi}{\sqrt{2}},\frac{\pi}{\sqrt{2}} \right)\mathbb{\rightarrow R},\ \ T_{1}(x) = {cotan}\left( \arccos\frac{x}{\sqrt{\pi^{2} - x^{2}}} \right), \] see Wolfram Alpha, cotan(arccos(x/(Sqrt(Pi^{}2-x^{}2)))),
\[T_{2}\mathbb{:R \rightarrow \lbrack -}\pi,\pi),\ \ T_{2}(x) = 2\arctan(x),\] one-to-one injective, see Wolfram Alpha, 2*arctan(x),
\[T_{3}\mathbb{:R \rightarrow}\lbrack - \pi,\pi),\ T_{3}(x) = \ \ \frac{2\pi}{1 + e^{x}} - \pi,\] one-to-one injective, see Wolfram Alpha, (2*Pi)/(1+exp(x))-Pi.
Does it make sense to denote the density function of the kime-phase differences \(\gamma = \phi - \theta\), by \(f_{\Gamma}(\gamma) = f_{\Gamma}(\phi_{i} - \phi_{j})\) or by just by \(f_{\Gamma}(\gamma) = f_{\Gamma}(i - j)\)? In other words, can we drop the explicit phase reference and denote the density of the kime-phase differences simply by \(f_{\Gamma}(\gamma)\mathbb{= E}\left( \phi(i),\phi(j) \right)\mathbb{\equiv E}(i,j) = R_{phase}(i,j)\)?
At a given kime, \(k = te^{i\varphi}\), we only have an analytic representation of the distribution of the kime-phase difference, not an explicit analytic function estimate of the actual kime-phase change or its derivative (as the kime-phase is random). In this situation, how can we define the partial derivative (or a measure-theoretic/probabilistic rate of change) of the original process of interest, \(f'(k)\)? Specifically, even though we know how to compute \(\frac{\partial f}{\partial r}\), we also need to estimate \(\frac{\partial f}{\partial\varphi}\) but only having the probability density function of the kime-phase differences} \(\Delta\varphi = \varphi_{i} - \varphi_{j}\), \(1 \leq i < j\). In other words, we can work with the distribution* of \(\Delta\varphi\), but can’t estimate a specific value for this rate of change, \(\Delta\varphi\).
Above discussed the notion of analyticity, which reflects functions \(f:\mathbb{C}\to \mathbb{C}\) whose Taylor series expansion about \(z_o\) converges to the functional value \(f(z)\) in some neighborhood \(\forall z\in \mathcal{N}_{z_o}\subseteq \mathbb{C}\).
The (physical) measurability property is associated with precise ability to track the value of some physical observable of a quantum system. This typically involves some kind of apparatus for measuring specific properties so that during the process of acquiring this information, the system survives the measurement after the measurement process is complete and an immediate follow up measurement of the same quantity would yield the same outcome, or at least a highly stable result. In reality, few experiments yield such measurements where the system being experimentally studied is physically unchanged, not modified, or preserved by the measurement process itself. Consider an extreme example of a Geiger counter measurement counting the number of positrons, or a more simple photo-detector counting the number of photons in a single mode cavity field. The tracked particles (photons or positrons) counted (tracked) by the device (photo-detector or Geiger counter) are actually absorbed to create a pulse, current, or an electrical impulse. While we count (track) the number of particles in the field, we transform the signal (particle dynamics) into some form of electrical data. After the measurement is completed, the state of the measurable cavity field is left in a base vacuum state \(|0\rangle\), which has no memory of the state of the field before the experiment took place and we recorded the digital measurement. However, we may be able redesign the experiment by using the observed outcome measurements (reflecting the number of particles detected in the cavity during a fixed time period) and every time we then measure the number of particles in the cavity, we always maintain the same level of particles in the cavity, \(n\). With such feedback loop, the experimental procedure can maintain the cavity field in a constant state of \(n\) particles, i.e., we can assign the state \(|n\rangle\) to the cavity field.
The Von Neumann Measurement Postulate reflects the state of the system immediately after a measurement is made.
The Von Neumann measurement postulate may be expressed using projection operators \(\hat{P}_n = |q_n\rangle\langle q_n|\), where the state of the system after the measurement yields the result \(q_n\) is
\[\frac{\hat{P}_n|\psi\rangle}{\sqrt{\langle \psi |\hat{P}_n|\psi\rangle}}= \frac{\hat{P}_n|\psi\rangle}{\sqrt{|\langle {q}_n|\psi\rangle |^2}}\ .\]
Note that the denominator normalizes the state after the measurement to unity.
This postulate suggests that after performing a measurement and observing a particular outcome, if we immediately repeat the experiment and obtain another measurement we can expect to observe the same result, indicating the system is in the associated eigenstate. The fact that the outcome scalar value is stable upon repeated measurement indicating that the system is consistent and reliably yields the same instantaneous value after a measurement.
In statistical terms, the Von Neumann measurement postulate parallels the notion of repeated measurements corresponding to multiple (large) samples from the kime-phase distribution \(\{\phi_i\}_i \sim \Phi[-\pi,\pi)\). All measurable properties represent pooled sample statistics, such as the sample arithmetic average, which has an expected value equal to the phase distribution mean,
\[\mathbb{E}\left (\frac{1}{n}\sum_{i=1}^n {\phi_i}\right )=\mu_{\Phi}\ .\]
In TCIU Chapter 6 (Kime-Phases Circular distribution), the agreement between observed analyticity and Von Neumann measurability of real observables and the stochastic quantum mechanical nature of the enigmatic kime-phases.
This connection between quantum mechanical and relativistic properties is rooted in the fact that physical observable processes rely on phase-aggregators to pool quantum fluctuations. This process effectively smooths and denoises the actual quanta characteristics and yields highly polished classical physical measurements. This duality between quanta (single experimental data point) and relativistic observables (sample/distribution driven evidence) is related to the effort to define a kime operator whose eigenvalues are observable scalar measures of kime-phase aggregators (e.g., sample phase mean) and their corresponding eigenvectors are the state-vectors (e.g., system experimental states).
Let’s expand on the notions of kime-phase observability, operator eigenspectra, and kime operator. We’ll attempt to draw direct parallels between the quantum physics (statistical mechanics) and data science (statistical inference). Starting with a physical question (e.g., computing the momentum of a particle) or a statistical inference problem (e.g., estimate a population mean), we can run as experiment to collect data, which can be used to calculate the expectation value \(X\), a physical observable (energy, momentum, spin, position, etc.) or another statistic (dispersion, range, IQR, 5th percentile, CDF, etc.)
\[ \mathbb E(X)= \sum_i x_i\,p(x_i)\ ,\] where \(\{x_i\}_i\) are the possible outcome values and \(\{p(x_i)\}_i\) are their corresponding probabilities of those outcomes. Indeed, if the outcomes are continuous, the expectation is naturally formulated in terms of an integral.
We can associate a vector \(|x_i\rangle\) to each of these outcome states \(x_i\) and since we assume IID sampling (random observations), we can consider orthonormal states, \(\langle x_i | x_j\rangle =\delta_{ij}\). As quantum mechanics and matrix algebra (linear models) are linear, a superposition state of outcomes (solutions) \(\{|x_i\rangle \}_{i}\) represents another a valid outcome solution \(|\phi \rangle =\sum_i {\alpha_i |x_i\rangle }\). By the Born rule postulate, the interpretation of the superposition state is extends the probability of observing a base state \(|x_i\rangle\), \(p(x_i)=|\alpha_i|^2\), where composite state \(|\phi \rangle\) is normalized, \[\sum_i {\alpha_i^2}=1\ .\]
By orthonormality and using Dirac notation, \[\alpha_i = \langle x_i|\psi \rangle,\ \alpha_i^* = \langle \psi|x_i \rangle,\ \forall i\ .\] By linearity, the expectation value of \(X\) is
\[\begin{align} \mathbb E(X) &= \sum_i {|\alpha_i |^2 x_i} = \sum_i {\alpha_i^* \alpha_i x_i} = \sum_i {\langle \psi|x_i\rangle \langle x_i|\psi\rangle x_i } \\ &=\langle \psi|\left( \underbrace{\sum_i |x_i\rangle\langle x_i|x_i}_{Operator,\ \hat{X}}\right)|\psi\rangle \equiv \langle \psi| \hat{X} |\psi\rangle \end{align} \ .\]
My the spectral theorem, any self-adjoint (Hermitian) matrix \(\hat{X}\) can be written as \[\hat{X} = \sum_i {\overbrace{\underbrace{|\lambda_i\rangle}_{eigenvectors} \langle \lambda_i|}^{Projection\ along\ |\lambda_i\rangle} \underbrace{\lambda_i}_{eigenvalues}}\ .\] The eigenvectors \(\{|\lambda_i\rangle\}_i\) form an orthonormal basis suggesting that only linear combinations of these eigenvectors of \(\hat{X}\) give the observable outcomes of experimental measurements. For instance, all vectors \(\gamma\) orthogonal to \(|\lambda_i\rangle\) will be annihilated. When a state \(|\gamma\rangle \perp |\lambda_i\rangle,\ \forall i\), the state can’t be written as a linear combination of the eigenvectors and it does not contribute to the expectation value.
The expectation value of an operator \(\hat {\mathcal{K}}\) corresponds to the average outcome of many repeated measurements of the observable associated with that operator,
\[\langle\hat {\mathcal{K}}\rangle_{\psi} \equiv \langle\hat {\mathcal{K}}\rangle =\langle\psi| \hat {\mathcal{K}}|\psi\rangle = \overbrace{\sum_{\alpha}{ \underbrace{\nu_{\alpha}\ }_{obs.value}\ \underbrace{\ |\langle \psi |\phi_{\alpha}\rangle|^2}_{(transition)probability}}} ^{weighted-average-of-outcomes }\ ,\] where \(\{|\phi_{\alpha}\rangle\}_{\alpha}\) is a complete set of eigenvectors for the observable operator \(\hat {\mathcal{K}}\), i.e., \(\hat {\mathcal{K}}|\phi_{\alpha}\rangle = \nu_{\alpha}| \phi_{\alpha}\rangle\).
In other words, if an observable quantity of a system is measured
many times, the average or expectation value
of these
measurements is given by the corresponding operator \(\hat {\mathcal{K}}\) acting on the state of
the system, \(|\psi\rangle\). The
operator expectation is important because it ties the physical
interpretation (the observable) with the mathematical representation of
the quantum state (the operator). Each observable is associated with a
specific self-adjoint (Hermitian) operator, which ensures that the
observable scalar values (the eigenvalues of the operator) are
real.
The expectation value of the kime-phase operator \(\hat {\mathcal{K}}_{\varphi}\) reflects the average of the kime-phase distribution, \(\varphi\sim\Phi[-\pi,\pi)\). For most kime-phase distributions, \(\langle\hat {\mathcal{K}}\rangle_{\psi}=\mu_{\Phi}=0\), reflecting the symmetry f the phase distribution. The expectation value of the operator \(\hat {\mathcal{K}}_{\varphi}\) in a quantum state described by a normalized wave function \(\psi(x)\) is given by the integral:
\[\langle\hat {\mathcal{K}}\rangle_{\psi}=\int_{-\pi}^{\pi}{\psi^*(x)\ \hat {\mathcal{K}}\ \psi(x)\ dx}\ .\]
The wave function \(\psi(x)\) encodes all the possible kime-phase states the system could be in and their corresponding probabilities. Each value \(x\) in the domain of the wave function is a possible kime-phase observation drawn from the phase distribution \(\Phi\). The complex-valued wavefunction portrays the probability amplitudes correlating with these potential measurements. Bounded operators acting on these phase state vectors help draw random kime-phases or measure the observable phases.
Caution: In the simulation below, to emphasize the differences between different observable processes, we are artificially offsetting (shifting) the random samples, and their corresponding mean trajectories (observable patterns). These offsets are both in the vertical (space) dimension as well as in the radial (phase) space. These offsets are are completely artificial and just intended to enhance the interpretation of the kime-phase sampling process by avoiding significant overlaps of points and observable kime-dynamic trends.
library(animation)
library(circular)
library(plotly)
epsilon <- 0.1
sampleSize <- 1000 # total number of phases to sample for 3 different processes (x, y, z)
sizePerTime <- 100 # number of phases to use for each fixed time (must divide sampleSize)
circleUniformPhi <- seq(from=-pi, to=pi, length.out=sizePerTime)
oopt = ani.options(interval = 0.2)
set.seed(1234)
# sample the the kime-phases for all 3 different processes and the r time points
x <- rvonmises(n=sampleSize, mu=circular(pi/5), kappa=3)
y <- rvonmises(n=sampleSize, mu=circular(-pi/3), kappa=5)
z <- rvonmises(n=sampleSize, mu=circular(0), kappa=10)
r <- seq(from=1, to=sampleSize/sizePerTime, length.out=10)
# Define a function that renormalizes the kime-phase to [-pi, pi)
pheRenormalize <- function (x) {
out <- ifelse(as.numeric(x) <= pi, as.numeric(x)+pi, as.numeric(x)-pi)
return (out)
}
# transform Von Mises samples from [0, 2*pi) to [-pi, pi)
x <- pheRenormalize(x)
y <- pheRenormalize(y)
z <- pheRenormalize(z)
# vectorize the samples
vectorX = as.vector(x)
vectorY = as.vector(y)
vectorZ = as.vector(z)
# Starting phases, set the first phase index=1
plotX = c(vectorX[1])
plotY = c(vectorY[1])
plotZ = c(vectorZ[1])
# # iterate over all time points (outer loop) and all phases (inner loop)
# for (t in 1:length(r)) { # loop over time
# for (i in 2:100) { # loop over kime-phases
# plotX = c(plotX,vectorX[i])
# plotY = c(plotY,vectorY[i])
# plotZ = c(plotZ,vectorZ[i])
#
# # Try to "stack=T the points ....
# #r1 = sqrt((resx$x)^2+(resx$y)^2)/2;
# #resx$x = r1*cos(resx$data)
# #resx$x = r1*cos(resx$data)
#
# tempX = circular(plotX[[1]])
# tempY = circular(plotY[[1]])
# tempZ = circular(plotZ[[1]])
# resx <- density(tempX, bw=25, xaxt='n', yaxt='n')
# resy <- density(tempY, bw=25, xaxt='n', yaxt='n')
# resz <- density(tempZ, bw=25, xaxt='n', yaxt='n')
# res <- plot(resx, points.plot=TRUE, xlim=c(-1.1,1.5), ylim=c(-1.5, 1.5),
# main="Trivariate random sampling of\n kime-magnitudes (times) and kime-directions (phases)",
# offset=0.9, shrink=1.0, ticks=T, lwd=3, axes=F, stack=TRUE, bins=150)
# lines(resy, points.plot=TRUE, col=2, points.col=2, plot.info=res, offset=1.0, shrink=1.45, lwd=3, stack=T)
# lines(resz, points.plot=TRUE, col=3, points.col=3, plot.info=res, offset=1.1, shrink=1.2, lwd=3, stack=T)
# segments(0, 0, r[i]*cos(vectorX[i]), r[i]*sin(vectorX[i]), lwd=2, col= 'black')
# segments(0, 0, r[i]*cos(vectorY[i]), r[i]*sin(vectorY[i]), lwd=2, col= 'red')
# segments(0, 0, r[i]*cos(vectorZ[i]), r[i]*sin(vectorZ[i]), lwd=2, col= 'green')
# ani.pause()
# }
# }
####################################################
# pl_list <- list()
pl_scene <- plot_ly(type='scatter3d', mode="markers")
plotX <- list()
plotY <- list()
plotZ <- list()
plotX_df <- list() # need separate dataframes to store all time foliations
plotY_df <- list()
plotZ_df <- list()
for (t in 1:length(r)) { # loop over time
# loop over kime-phases
plotX[[t]] <- as.numeric(x[c(( (t-1)*length(r) + 1):((t-1)*length(r) + sizePerTime))])
plotY[[t]] <- as.numeric(y[c(( (t-1)*length(r) + 1):((t-1)*length(r) + sizePerTime))])
plotZ[[t]] <- as.numeric(z[c(( (t-1)*length(r) + 1):((t-1)*length(r) + sizePerTime))])
# Try to "stack=T the points ....
#r1 = sqrt((resx$x)^2+(resx$y)^2)/2;
#resx$x = r1*cos(resx$data)
#resx$x = r1*cos(resx$data)
tempX = circular(unlist(plotX[[t]]))
tempY = circular(unlist(plotY[[t]]))
tempZ = circular(unlist(plotZ[[t]]))
resx <- density(tempX, bw=25, xaxt='n', yaxt='n')
resy <- density(tempY, bw=25, xaxt='n', yaxt='n')
resz <- density(tempZ, bw=25, xaxt='n', yaxt='n')
# res <- plot(resx, points.plot=TRUE, xlim=c(-1.1,1.5), ylim=c(-1.5, 1.5),
# main="Trivariate random sampling of\n kime-magnitudes (times) and kime-directions (phases)",
# offset=0.9, shrink=1.0, ticks=T, lwd=3, axes=F, stack=TRUE, bins=150)
# pl_list[[t]]
unifPhi_df <- as.data.frame(cbind(t=t, circleUniformPhi=circleUniformPhi))
plotX_df[[t]] <- as.data.frame(cbind(t=t, plotX=unlist(plotX[[t]])))
plotY_df[[t]] <- as.data.frame(cbind(t=t, plotY=unlist(plotY[[t]])))
plotZ_df[[t]] <- as.data.frame(cbind(t=t, plotZ=unlist(plotZ[[t]])))
pl_scene <- pl_scene %>% add_trace(data=unifPhi_df, showlegend=FALSE,
x = ~((t-epsilon)*cos(circleUniformPhi)),
y = ~((t-epsilon)*sin(circleUniformPhi)), z=0,
name=paste0("Time=",t), line=list(color='gray'),
mode = 'lines', opacity=0.3) %>%
add_markers(data=plotX_df[[t]], x=~(t*cos(plotX)), y=~(t*sin(plotX)), z=0,
type='scatter3d', name=paste0("X: t=",t),
marker=list(color='green'), showlegend=FALSE,
mode = 'markers', opacity=0.3) %>%
add_markers(data=plotY_df[[t]], x=~((t+epsilon)*cos(plotY)),
y=~((t+epsilon)*sin(plotY)), z=0-epsilon, showlegend=FALSE,
type='scatter3d', name=paste0("Y: t=",t),
marker=list(color='blue'),
mode = 'markers', opacity=0.3) %>%
add_markers(data=plotZ_df[[t]], x=~((t+2*epsilon)*cos(plotZ)),
y=~((t+2*epsilon)*sin(plotZ)), z=0+epsilon, showlegend=FALSE,
type='scatter3d', name=paste0("Z: t=",t),
marker=list(color='red'),
mode = 'markers', opacity=0.3)
}
means_df <- as.data.frame(cbind(t = c(1:length(r)),
plotX_means=unlist(lapply(plotX, mean)),
plotY_means=unlist(lapply(plotY, mean)),
plotZ_means=unlist(lapply(plotZ, mean))))
pl_scene <- pl_scene %>%
# add averaged (donoised) phase trajectories
add_trace(data=means_df, x=~(t*cos(plotX_means)),
y=~(t*sin(plotX_means)), z=0,
type='scatter3d', showlegend=FALSE, mode='lines', name="Expected Obs X",
line=list(color='green', width=15), opacity=0.8) %>%
add_trace(data=means_df, x=~(t*cos(plotY_means)),
y=~(t*sin(plotY_means)), z=0-epsilon,
type='scatter3d', showlegend=FALSE, mode='lines', name="Expected Obs X",
line=list(color='blue', width=15), opacity=0.8) %>%
add_trace(data=means_df, x=~(t*cos(plotZ_means)),
y=~(t*sin(plotZ_means)), z=0+epsilon,
type='scatter3d', showlegend=FALSE, mode='lines', name="Expected Obs X",
line=list(color='red', width=15), opacity=0.8) %>%
add_trace(x=0, y=0, z=c(-2,2), name="Space", showlegend=FALSE,
line=list(color='gray', width=15), opacity=0.8) %>%
layout(title="Pseudo Spacekime (1D Space, 2D Kime) Kime-Phase Sampling and Foliation",
scene = list(xaxis=list(title="Kappa1"), yaxis=list(title="Kappa2"),
zaxis=list(title="Space"))) %>% hide_colorbar()
pl_scene
# pl_list %>%
# subplot(nrows = length(pl_list)) %>%
# layout(title="Integral Approximations by Riemann Sums", legend = list(orientation = 'h'))
In this simulation, the \(3\) spatial dimensions \((x,y,z)\in\mathbb{R}^3\) are compressed into \(1D\) along the vertical axis \(z\in\mathbb{R}^1\), the radial displacement represents the time dynamics, \(t\in \mathbb{R}^+\), and the angular scatters of three different processes, representing repeated measurements from \(3\) different circular distributions at a fixed spatiotemporal location, are shown in different colors. Mean-dynamics across time of the three different time-varying distributions are shown as smooth curves color-coded to reflect the color of their corresponding process distributions, samples, and sampling-distributions.
Data from repeated measurements are acquired in time-dynamic experiments tracking univariate or multivariate processes repeatedly over time. Hence, repeated measurement data are tracked over multiple points in time. The simple situation when a single response variable is measured over time (univariate case) is easier to describe the challenges, but the mode general multivariate response is more realistic in many data science applications.
In general, the responses over time may be heavily temporally correlated, yet still corrupted by noise modeled by the kime-phase distribution. Note that individual observations collected at points in time close together are likely to be more similar to one another compared to observations collected at distant times. However, even at the same time point, the phase distribution controls the level of expected dispersion/variability between multiple measurements obtained under identical experimental conditions (at the fixed point in time). Classical multivariate analysis treats observations collected at different time points as (potentially correlated) different variables.
Consider a completely randomized block design experiment to determine the effects of 5 surgical treatments on cardiovascular health outcomes where a group of \(35\) patients grouped into each of the \(5\) treatment groups of sizes \(\{7, 6, 9, 5, 8\}\). the health outcomes of every patient \(\{p_j\}_{j=1}^{35}\) are recorded at \(10\) consecutive time points, \(\{t_k\}_{k=0}^9\) following one of the available \(5\) clinical treatment regimens, \(\{R_i\}_{i=1}^5\). Note that the treatment group sizes are \(\{|R_1|=7,|R_2|=6,|R_3|=9,|R_4|=5,|R_5|=8 \}\).
There are many alternative strategies to analyze data from such repeated sample experiment. analysis of variance (ANOVA) and multivatiate analysis of variance (MANOVA) represent a pair of common parametric techniques.
Let’s make the following notations:
ANOVA provides a parametric statistical test (using the Fisher’s F distribution) to study interactions between treatment regimens and time, using the main effects of treatments and time. MANOVA facilitates testing for significant interaction effects between treatment regimens and time, and for main effects of treatment regimen;
The ANOVA (linear) model \[Y_{ijk}=\mu+\alpha_i+\beta_{j(i)}+\tau_k+(\alpha \tau)_{ik} + \epsilon_{ijk},\] assumes that the health outcome data \(Y_{ijk}\) for the \(i^{th}\) treatment regimen \(R_i\) for the \(j^{th}\) patient \(p_j\) at \(k^{th}\) time point \(t_k\) is a mixture of some (yet unknown) overall mean response \(\mu\), plus some the treatment effect \(\alpha_i\), the effect of the patient within that treatment regimen \(\beta_{j(i)}\), the effect of time \(\tau_k\), the effect of the interaction between treatment and time \((\alpha \tau)_{ik}\), and a residual error \(\epsilon_{ijk}\sim N(0,\sigma_{\epsilon}^2)\).
The explicit ANOVA assumptions include:
This ANOVA mixed-effects model includes a random effect of the patient and fixed effects for treatment regimen and time. The Fisher’s F-test is used to carry ANOVA using the following variance decomposition summary.
Source of Variation | Sum of Squares (SS) | Degrees of Freedom (df) | Mean Squares (MS) | F-Ratio |
---|---|---|---|---|
Between samples (Treatment) | \(SS_{treat}\) | a - 1 | \(\frac{SS_{treat}}{a-1}\) | \(\frac{SS_{treat}}{MS_{error(a)}}\) |
Within samples (Time) | \(SS_{time}\) | t - 1 | \(\frac{SS_{time}}{t-1}\) | \(\frac{SS_{time}}{MS_{error(b)}}\) |
Interaction (Treat*Time) | \(SS_{time}\) | (a-1)(t - 1) | \(\frac{SS_{treat\times times(b)}}{(a-1)(t-1)}\) | \(\frac{MS_{treat\times times(b)}}{MS_{error(b)}}\) |
Error(a) | \(SS_{error(a)}\) | N-a | \(\frac{SS_{error(a)}}{N-a}\) | |
Error(b) | \(SS_{error(b)}\) | (N-b)(t-1) | \(\frac{SS_{error(b)}}{(N-a)(t-1)}\) | |
Total | \(SS_{total}\) | Nt - 1 |
where,
Sum of Squares \(SS\) formulas are given below:
\[\begin{array}{lll}SS_{total}& =& \sum_{i=1}^{a}\sum_{j=1}^{n_i}\sum_{k=1}^{t}Y^2_{ijk}-Nt\bar{y}^2_{...}\\SS_{treat} &= &t\sum_{i=1}^{a}n_i\bar{y}^2_{i..} - Nt\bar{y}^2_{...}\\SS_{error(a)}& =& t\sum_{i=1}^{a}\sum_{j=1}^{n_i}\bar{y}^2_{ij.} - t\sum_{i=1}^{a}n_i\bar{y}^2_{i..}\\SS_{time}& =& N\sum_{k=1}^{t}\bar{y}^2_{..k}-Nt\bar{y}^2_{...}\\SS_{\text{treat x time}} &=& \sum_{i=1}^{a}\sum_{k=1}^{t}n_i\bar{y}^2_{i.k} - Nt\bar{y}^2_{...}-SS_{treat} -SS_{time}\end{array}\ .\]
The statistical inference based on the ANOVA (variance decomposition) model includes:
\[F = \overbrace{\frac{MS_{\text{treat x time}}}{MS_{error(b)}}}^{obs.\ test\ statistic} > \overbrace{F_{\underbrace{(a-1)(t-1)}_{df_1}, \underbrace{(N-a)(t-1)}_{df_2}, \underbrace{\alpha}_{signif.}}}^{F-distr.\ critical\ value}\]
\[F = \dfrac{MS_{treat}}{MS_{error(a)}} > F_{a-1, N-a, a}\ .\]
\[F = \frac{MS_{time}}{MS_{error(b)}} > F_{t-1, (N-a)(t-1), \alpha}\ .\]
\[\mathbf{Y}_{ij} = \left(\begin{array}{c}Y_{ij1}\\ Y_{ij2} \\ \vdots\\ Y_{ijt}\end{array}\right)\ .\]
Effectively, we treat the data at different time-points as different variables or features. One-way MANOVA assumptions include:
MANOVA inference includes:
\[\mathbf{Z}_{ij} = \left(\begin{array}{c}Z_{ij1}\\ Z_{ij2} \\ \vdots \\ Z_{ij, t-1}\end{array}\right) = \left(\begin{array}{c}Y_{ij2}-Y_{ij1}\\ Y_{ij3}-Y_{ij2} \\ \vdots \\Y_{ijt}-Y_{ij,t-1}\end{array}\right)\ .\]
The (random) vector \(\mathbf{Z}_{ij}\) is a function of the random data and we can compute the expected population mean for treatment \(i\), \(E(\mathbf{Z}_{ij}) = \boldsymbol{\mu}_{Z_i}\). The MANOVA test on these \(\mathbf{Z}_{ij}\)’s is based the null hypothesis \(H_o\colon \boldsymbol{\mu}_{Z_1} = \boldsymbol{\mu}_{Z_2} = \dots = \boldsymbol{\mu}_{Z_a}\).
The key factor exploited by linear models such as ANOVA/MANOVA is the powerful analytical decomposition of the observed process variation, which quantifies the type and amount of complementary (orthogonal) dispersion observed in data samples. Recall that \(N=\sum_i^a{n_i}\) and
\[\bar{y}_{...} = \frac{1}{Nt} \sum_{i=1}^{a}\sum_{j=1}^{n_i}\sum_{k=1}^{t} {Y_{ijk}}, \]
\[SS_{total} = \frac{1}{n-1} \sum_{i=1}^{a}\sum_{j=1}^{n_i}\sum_{k=1}^{t} {(Y^2_{ijk}- \bar{y})^2} = \sum_{i=1}^{a}\sum_{j=1}^{n_i}\sum_{k=1}^{t}Y^2_{ijk}-Nt\bar{y}^2_{...}\ .\]
Then we can expand the total variance as a sum of complementary factors
\[\begin{array}{lll} \underbrace{SS_{total}}_{\text{Overall Obs. Variance}} & =& \sum_{i=1}^{a}\sum_{j=1}^{n_i}\sum_{k=1}^{t}Y^2_{ijk}-Nt\bar{y}^2_{...}\\ \underbrace{SS_{treat}}_{\text{Obs. Var. due to Treatment}} &= & t\sum_{i=1}^{a}n_i\bar{y}^2_{i..} - Nt\bar{y}^2_{...} \\ \underbrace{SS_{time}}_{\text{Obs. Var. due to Time}} & =& N\sum_{k=1}^{t}\bar{y}^2_{..k}-Nt\bar{y}^2_{...} \\ \underbrace{SS_{treat \times time}}_{\text{Obs. Var. due to Treat*Time Interaction}} &=& \sum_{i=1}^{a}\sum_{k=1}^{t}n_i\bar{y}^2_{i.k} - Nt\bar{y}^2_{...}-SS_{treat} -SS_{time} \\ \underbrace{SS_{error(a)}}_{\text{Residual Var. (unaccounted by model)}} & =& t\sum_{i=1}^{a}\sum_{j=1}^{n_i}\bar{y}^2_{ij.} - t\sum_{i=1}^{a}n_i\bar{y}^2_{i..}\end{array}\ .\]
\[SS_{total} = \sum_{i=1}^{a}\sum_{j=1}^{n_i}\sum_{k=1}^{t} Y^2_{ijk}-Nt\bar{y}^2_{...} = \\ \underbrace{t\sum_{i=1}^{a}n_i\bar{y}^2_{i..} - Nt\bar{y}^2_{...}}_{SS_{treat}} + \underbrace{N\sum_{k=1}^{t}\bar{y}^2_{..k}-Nt\bar{y}^2_{...}}_{SS_{time}} +\\ \underbrace{\sum_{i=1}^{a}\sum_{k=1}^{t}n_i\bar{y}^2_{i.k} - Nt\bar{y}^2_{...}-SS_{treat} -SS_{time}}_{SS_{{treat \times time}}} +\\ \underbrace{t\sum_{i=1}^{a}\sum_{j=1}^{n_i}\bar{y}^2_{ij.} - t\sum_{i=1}^{a}n_i\bar{y}^2_{i..}}_{SS_{error(a)}}\ .\]
To draw the connection between observed variability decomposition and complex-time we will consider a study of single-subject task-based fMRI activation across multiple sessions. In this study, the researchers proposed a within-subject variance decomposition of a task-based fMRI using a large number of repeated measures, \(500\) trials of three different subjects undergoing \(100\) functional scans in \(9-10\) different sessions. In this case, the observed within-subject variance was segregated into \(4\) primary components - variance across-sessions, variance across-runs within a session, variance across-blocks within a run, and residual measurement (model-unaccounted error).
Within-subject variance decomposition:
That study reported that across \(16\) cortical networks, \(\sigma^2_{block}\) contributions to within-subject variance were larger in high-order cognitive networks compared to somatosensory brain networks. In addition, \(\sigma^2_{block} \gg \sigma^2_{session}\) in higher-order cognitive networks associated with emotion and interoception roles.
This spatiotemporal distribution factorization of the total observed within-subject variance illustrates the importance of identifying dominant variability components, which subsequently may guide prospective fMRI study-designs and data acquisition protocols.
All functional MRI runs had the same organization of blocks as shown in the diagram below.
# install.packages("DiagrammeR")
# See DiagrammeR Instructions:
# https://epirhandbook.com/en/diagrams-and-charts.html
library(DiagrammeR)
graphStr <- "# Mind commented instructions
digraph fMRI_Study_Design { # 'digraph' means 'directional graph', then the graph name
################# graph statement #################
graph [layout = dot, rankdir = LR, # TB = layout top-to-bottom
fontsize = 14]
################# nodes (circles) #################
node [shape = circle, fixedsize = true, width = 1.8]
Start [label = 'fMRI Run\nStart']
Begin [label = 'Begining\nRest period\n(30 seconds)']
Task [label = 'Task block\n(20 seconds)']
Rest [label = 'Rest block\n(20 seconds)', fontcolor = darkgreen]
End [label = 'Ending\nRest block\n(20 seconds)', fontcolor=darkgreen]
Total [label = 'Run Total\n(340 seconds)', fontcolor = darkgreen]
################ edges #######
Start -> Begin [label='initialize', fontcolor=red, color=red]
Begin -> Task [label='task', fontcolor=red, color=red]
Task -> Rest [label='Loop', fontcolor=red, color=red]
Rest -> End [label='wrap up', fontcolor=red, color=red]
End -> Total [label='Complete', fontcolor=red, color=red]
################# grouped edge #################
{Rest} -> Task [label = 'Repeat\nTask-Rest Block\n5 times', fontcolor=darkgreen,
color = darkgreen, style = dashed]
}
"
grViz(graphStr)
Suppose the complex-valued fMRI signal is \(Y_{session=i,run=j,block=k}\) and we fit a linear model decomposing the effect estimates
\[\underbrace{\hat{\beta}_{i(j(k))}}_{effect\\ estimate}=\underbrace{\alpha}_{overall\\ mean} + \underbrace{\theta_k}_{session\\ effect} + \underbrace{\xi_{j(k)}}_{run\\effect} + \underbrace{\eta_{i(j(k))}}_{block\\ effect} + \underbrace{\varepsilon_{i(j(k))}}_{measurement\\ error } \ .\]
This partitioning of the effects is in terms of the indices \(block=i\), \(run=j\), and \(day=k\), where parentheses \(()\) indicate the nesting structure between consecutive levels. We assume Gaussian distributions for all random effects and the residual error,
\[\theta_k\sim N(0,\sigma^2_{session}),\ \ \xi_{j(k)}\sim N(0,\sigma^2_{run}),\ \ \eta_{i(j(k))}\sim N(0,\sigma^2_{block}),\ \ \varepsilon_{i(j(k))}\sim N(0,\sigma^2_{i})\ .\]
Then, the variance decomposition of the effect estimate model is
\[Var(\hat{\beta}_{i(j(k))}) = \sigma^2_{session} + \sigma^2_{run} + \sigma^2_{block} + \sigma^2_{error}\ .\]
Question: Where do the random effects (and their estimates) come from?
Answer: As shown in this article, we start with a general linear model (GLM) \(Y = X\beta + \varepsilon\) representing the observed fMRI response (intensity) \(y\) at each voxel as a linear combination of some explanatory variables, columns in the design matrix \(X\). Each column in \(X\) corresponds to one effect that can be experimentally manipulated or that may confound the observed outcome \(y\).
The GLM expresses the response variable \(y\) as a mixture of explanatory variables and some residual error term \(\varepsilon\sim N(0,\Sigma)\). fMRI pre-processing protocols may filter the data with a convolution or residual forming matrix \(S\), leading to a generalized linear model that includes intrinsic serial correlations and applied extrinsic filtering. Different choices of \(S\) correspond to different estimation schemes
\[\underbrace{Y}_{observed\\ fMRI} = \overbrace{X}^{design\\ matrix} \underbrace{\beta}_{effects\\ vector} + \overbrace{\varepsilon}^{residual\\ error}\ .\]
\[S=\begin{cases} \Sigma^{-\frac{1}{2}}, & whitening \\ 1, & none \\ 1-X_oX_o^+, & adjustment \end{cases}\ .\]
\[SY = SX\beta + S\varepsilon, \ \ V=S\Sigma S',\ \ \underbrace{\hat{\beta}=(SX)^+ SY}_{\text {random effect estimates}},\ \ \frac{c'\hat{\beta}}{\sqrt{\hat{V}(c'\hat{\beta})}}\sim t_{df},\]
\[\hat{V}(c'\hat{\beta})=\hat{\sigma}^2 c'(SX)^+ V((SX)^+)',\\ R=I-SX(SX)^+, \ \ \hat{\sigma}^2=\frac{Y'RY}{trace(RV)}, df=\frac{trace(RV)^2}{trace(RVRV)}\ .\]
The effect parameter estimates are obtained using least squares based on the matrix-pseudo-inverse of the filtered design matrix, denoted by \(^+\). The most general effect of interest is specified by a vector of contrast weights \(c\) that give a weighted sum or compound of parameter estimates, referred to as a contrast. However, often we use unitary contrasts corresponding to a single variable in the design matrix \(X\), e.g., \(c=(1,0,\cdots,0)\) corresponds to the effect of the first covariate (first column in \(X\)). Once the effects are estimated (we have computed the LS estimates of \(\hat{\beta}_c\), we can compute the \(T_c\sim t_{df}\) statistic to assess the statistical significance of the effect.
In essence, fMRI analyses typically involve two phases:
\[({\text{First-level analysis}})\ \underbrace{Y_m=X_m \beta_m + \varepsilon_m}_{Individual\ Subject\ Model}, \ \underbrace{Y_m=\{Y_{m,t=0},Y_{m,t=1},\cdots, Y_{m,t=T}\}}_{fMRI\ series\ signal},\ \forall \ (subject)\ m,\\ ({\text{Second-level analysis}})\ \ \underbrace{\beta_m=X_{gm} \beta_g +\varepsilon_{gm}}_{Population\ Model}, \ \ \underbrace{\beta_g=\{\beta_{mo},\beta_{m1},\cdots, \beta_{mG}\}}_{subject-level-1\ regression\ coefficients\\ (effects)\ \forall EV\ column\ in\ X},\ \forall g\in G\ (group)\ .\]
This two-level fMRI analysis can also be represented into a General Linear Mixed Effects Model:
\[Y_m=\underbrace{\beta_o+X_{m1} \beta_1+\cdots+X_{mp} \beta_p}_{fixed-effects}+ \underbrace{\beta_{mo}+Z_{m1} \beta_{m1}+\cdots + Z_{mq} \beta_{mq}}_{random-effects} +\varepsilon_m\ , \forall \ (subject)\ m .\]
Again, \(Y_m=\{Y_{m,t=0},Y_{m,t=1},\cdots, Y_{m,t=T}\}\) is a column vector of observed fMRI time-series for individual \(m\) at a fixed 3D voxel location, \(\{X_{mj}\}\) is the vector of the \(p\) regressors included in the linear (LME) model, \(\{\beta_o,\cdots, \beta_p\}\) are the fixed-effect regression coefficients, which are identical across all subjects, the column vectors \(\{Z_{mq}\}_q\) are the random effect regressors with corresponding random effects coefficients \(\{\beta_{mq}\}_q\sim N(0,\Psi)\), and \(q\) is the total number of random effects included in the model. The random effects capture the variability across subjects for each of the regressors \(\{Z_{mq}\}_q\), \(\Psi\) is a \((q+1)\times (q+1)\) variance-covariance matrix, and the \(n_m\times 1\) vector of within-group errors \(\varepsilon_m=\{\varepsilon_{m1},\cdots,\varepsilon_{mn_m}\}\sim N(0,\sigma^2 V)\), where \(\sigma^2\) and \(V\) are the variance and the correlation matrix.
The two-tier fMRI analysis can be formulated into a generalized linear mixed-effects model as follows. For simplicity of notation, assume that the first level design matrix for subject \(m\), \(X_m\), has only four columns representing an intercept and \(3\) covariates, \(EV_1, EV_2, EV_3\), and there are \(T\) temporal volumes in the 4D fMRI spatio-temporal data
\[X_m=\begin{pmatrix} 1 & EV_{11} & EV_{21} & EV_{31} \\ 1 & EV_{12} & EV_{22} & EV_{32} \\ \vdots & \vdots & \vdots & \vdots \\ 1 & EV_{1T} & EV_{2T} & EV_{3T} \end{pmatrix}\equiv [X_{mo}\ X_{m1}\ X_{m2}\ X_{m3}\ ]\ ,\]
where \(\forall m,\ X_{mj},\ 0\leq j \leq 3\) is a column vector of ones (for the intercept, \(j=0\)), or \(EV_j, \ 1\leq j \leq 3\). The vector of regression coefficients corresponding to each explanatory variable is
\[\beta_m=\begin{pmatrix} \beta_{mo} \\ \beta_{m1} \\ \beta_{m2} \\ \beta_{m3} \end{pmatrix} ,\]
Suppose that in the second-level analysis we want to estimate the average population-wide effect of each EV for the whole group. This corresponds to using hte following design matrix
\[X_g=\begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}\equiv I_{4\times 4}\ .\]
Next we plug this the second-level population model \(\beta_m=X_{gm} \beta_g +\varepsilon_{gm}\) into the first-level model
\[\overbrace{Y_m=X_m \beta_m +\varepsilon_{m}}^{First \ Level\ (individual)}\\ =X_m \underbrace{\left ( X_{gm} \beta_g +\varepsilon_{gm}\right )}_{ Second\ Level\ (population)} +\varepsilon_{m}=\\ [X_{mo}\ X_{m1}\ X_{m2}\ X_{m3}\ ] \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \beta_{go} \\ \beta_{g1} \\ \beta_{g2} \\ \beta_{g3} \end{pmatrix} + [X_{mo}\ X_{m1}\ X_{m2}\ X_{m3}\ ] \begin{pmatrix} \varepsilon_{gmo} \\ \varepsilon_{gm1} \\ \varepsilon_{gm2} \\ \varepsilon_{gm3} \end{pmatrix} +\varepsilon_{m}=\\ \underbrace{\underbrace{X_{mo}\beta_{go}+X_{m1}\beta_{g1}+X_{m2}\beta_{g2}+X_{m3}\beta_{g3}}_{fixed\ effects}+ \underbrace{X_{mo}\varepsilon_{gmo}+X_{m1}\varepsilon_{gm1}+X_{m2}\varepsilon_{gm2}+ X_{m3}\varepsilon_{gm3}}_{random\ effects}+ \underbrace{\varepsilon_m}_{residual\ errors}}_{General\ LME} \ ,\] where \(\varepsilon_m\sim N(0,\sigma^2V)\), \(\varepsilon_g\sim N(0,R_g)\), and \(R_g\) is the variance-covariance matrix of the group-level errors.
The intrinsic process error \(\varepsilon_m\sim N(0,\sigma^2V)\) reflects classical variability in the underlying process distribution. The more subtle error-tensor \(\varepsilon_g\sim N(0,R_g)\), where \(\varepsilon_m\sim N(0,\sigma^2V)\), depends on the variance-covariance matrix of the group-level errors \(R_g\) and can be interpreted as repeated-measures variability due to the kime-phase distribution. Hence, the kime-phases can be used to model or track repeated-measure variations. Recall this kime-phase simulation we discussed earlier demonstrating alternative kime-phase distributions coupled with corresponding mean-phase-aggregators (solid radial lines) representing real observable (deterministic, not stochastic) measurable quantities. The latter avoid the intrinsic kime-phase distribution noise due to random sampling through the phase-space distribution \(\Phi_{[-\pi,\pi))}\) when taking real measurements depending de facto on some phase-aggregating kernels that smooth the intractable (unobservable) and intrinsically-random kime-phases.
In 1915, Albert Einstein proposed the general theory of relativity (GTR). The following year, in 1916, he proposed the following three tests to validate or falsify GTR. Some of these GTR hypotheses were confirmed using data acquired by 1919.
Similarly, to test, affirm, or reject the complex-time model, i.e., the existence of kime-phases and the interpretation of the spacekime representation of repeated longitudinal measurements, we need to propose testable hypotheses that can be used to determine the validity of the spacekime representation. Ideally, we need to express such hypotheses as if-then statements indicating that under certain conditions or actions (investigator controlled, independent variables, IVs), specific outcomes or results (dependent variables, DVs) are to be expected. We may need to propose experiments that can be carried to prove or disprove the existence of kime-phases, with substantial reproducibility.
In essence, there may be viable physics experiments or statistical-analytical experiments that can be formulated to test the spacekime representation.
For instance, potential physics experiments to consider include:
On the other hand, potential statistical and analytical experiments to test the validity of the kime-phase as a stochastic representation of repeated sampling may rely on tests confirming or disproving scientific inference improvements based on the spacekime-model compared to classical approaches. For example:
In quantum superposition particles are also distributions always maintaining simultaneously all possible basis states (realities) until the moment they’re sampled, drawn, measured, or detected, i.e., observed as scalars, vectors, or tensors. Similarly in statistics, data science, and AI, we often collect, sample, simulate, or model data as proxy of a physical phenomenon we are studying. The final data-driven scientific inference about the phenomenon is expected to be reproducible, consistent, unbiased, accurate, and reliable, reflecting the process distribution properties despite the fact that all data-driven inference is sample evidence-based, whereas the underlying process distributions are generalized functions. This statistical analytic-duality resembles the particle-wave duality in physics.
In applied sciences, sample statistics, such as the sample mean, variance, percentiles, range, inter-quartile range, etc., are always well-defined. Their values are not known at any given time until the (finite) sample \(\{x_i\}_{i=1}^n\) is collected (observations are recorded). For instance the sample mean statistic is the arithmetic average of all observations, \(\bar{x}_n=\frac{1}{n}\sum_{i-1}^n {x_i}\), yet, not knowing the value of the sample-mean prior to collecting the data is simply due to lack of evidence, since we have a closed form expression of how to obtain the sample average statistic estimating the population (physical system’s) expected mean response, \(\mu_X\ \underbrace{\leftarrow}_{{n\to\infty}}\ \bar{x}_n\). Quantum systems interact in ways that can be explained with superpositions of different discrete base-states, and quantum system measurements yield statistical results corresponding to any one of the possible states appearing at random.
Similarly, analytical systems interact in ways that can be explained with superpositions of different, discrete, and finite samples, and analytical system measurements yield statistics corresponding to any one of the specific sample-statistic outcomes, which vary between samples and experiments. However, the first two fundamental laws of probability theory governing all statistical inference (CLT and LLN) yield that this between sample variability of sample-statistics is expected to decrease rapidly with the increases of the sample size. For instance, IID samples \(\{x_i\}_{i=1}^n\sim \mathcal{D}(\mu_X,\sigma_X^2)\), drawn from most nice distributions \(\mathcal{D}\) with mean \(\mu_X\) and variance \(\sigma_X^2\), yield sample-means \(\bar{x}_n\) whose sampling distribution converges (in distribution) to a Normal distribution, an asymptotic limiting distribution with rapidly decaying variance as the sample size increases
\[\underbrace{\bar{x}_n }_{Sampling\\ distirbution}\ \ \ \overbrace{\underbrace{\longrightarrow}_{n\to\infty}}^{Convergence\\ in\ distribution} \underbrace{\mathcal{N}\left(\mu_X,\frac{\sigma_X^2}{n}\right)}_{Asymptotic\\ distribution}\ .\]
Let’s try to explicate the duality between sampling the kime-phase distribution and the complementary sampling of the process state space (spatiotemporal sampling).
Suppose the process \(X\sim F_X\), where \(F_X\) is an arbitrary distribution and let \(\sigma^2=Var(X)\). If \(\{X_1,\cdots ,X_n\}\sim F_X\) is an IID sample the CLT suggests that under mild assumptions \(\bar{X}\longrightarrow N(\mu,\sigma^2/n)\) as \(n\to\infty\). In repeated spatiotemporal sampling, \(\forall\ n\), we observe \(K\) (longitudinal) sample means \(\{\bar{X}_{n1},\cdots, \bar{X}_{nK}\}\), where \(\forall 1\leq k\leq K\) each IID sample \(\{X_{1k},\cdots ,X_{nk}\}\sim F_X\) (across the kime-phase space) allows us to compute the \(k^{th}\) sample mean \(\bar{X}_{nk}\).
Observe that computing \(\bar{X}_{nk}\) from the \(k^{th}\) sample \(\{X_{1k},\cdots ,X_{nk}\}\) is identical to directly sampling \(\bar{X}_{nk}\) from the distribution \(F_{\bar{X}_n}\). Let’s reflect on this set up which clearly involves \(2\) independent sampling strategies.
Note that the distribution \(F_{\bar{X}_n}\) need not be Normal, it may be, but in general it won’t be Normal. For instance suppose our spatiotemporal sampling scheme involves exponential distribution, i.e., \(\{X_{1k_o},\cdots , X_{nk_o}\}\sim f_X\equiv Exp(\lambda)\equiv\Gamma(1,\lambda)\) and \(\sum_{i=1}^n{X_{ik_o}}\sim \Gamma\left(n,\lambda\right )\),
\[f_{Exp{(\lambda)}}(x)=\lambda e^{-\lambda x}\ \ ; \ \mathbb{E}(X_{Exp{(\lambda)}})=\frac{1}{\lambda} \ \ ; \ Var(X_{Exp{(\lambda)}})= \frac{1}{\lambda^2}\ ; \\ f_{\Gamma(shape=\alpha,scale=\lambda)}(x)=\frac{x^{\alpha-1}\lambda^{\alpha}}{\Gamma(\alpha)}e^{-\lambda x}\ \ ; \ \mathbb{E}(X_{\Gamma(\alpha,\lambda)})=\alpha\lambda \ \ ; \ Var(X_{\Gamma(\alpha,\lambda)})= \alpha\lambda^2\ \ .\]
Based on this sample, the sampling distribution of \(\bar{X}_{nk_o}\equiv \frac{1}{n}\sum_{i=1}^n{X_{ik_o}}\sim Exp\left (\lambda\right )\equiv\Gamma\left(n,n\lambda\right )\), since a sum of exponential IID \(Exp(\lambda)\) variables is gamma distributed, \(\sum_{i=1}^n{X_{ik_o}}\sim \Gamma\left(n,\lambda\right )\) and scaling a random variable \(Y=cX\) yields \(F_Y(x)=\frac{1}{c}f_X(\frac{x}{c})\), in our case the constant multiplier \(c=\frac{1}{n}\), and \(Y=\frac{1}{n}\sum_{i=1}^n{X_i}\).
Of course, as \(n\to\infty\), \(\Gamma\left(n,\frac{\gamma}{n}\right )\to N(0,1)\), yet for any fixed \(n\), the distribution is similar, but not identical to normal.
Prior to data acquisition, \(\bar{X}_{nk_o}\) is a random variable, once the observed data values are plugged in, its a constant. Hence, the sample mean random variable \(\bar{X}_{nk_o}\) based on \(\{X_{1k_o},\cdots , X_{nk_o}\}\sim F_X\), and the random variable \(Y\sim F_{\bar{X}_{nk_o}}\) represent exactly the same random variable. In other words, drawing \(K\) samples of IID observations \(\{X_{1k_o},\cdots , X_{nk_o}\}\sim F_X\) and computing \(\bar{X}_{nk_o}=\frac{1}{n}\sum_{i=1}^n{X_{ik_o}}\) is equivalent to drawing \(K\) samples directly from \(F_{\bar{X}_{nk_o}}\).
Below is an example of a 3D Brownian motion/Wiener process as an example of a stochastic random walk process. In this example, the Wiener process is intentionally disturbed by random Poisson noise, which leads to occasional rapid and stochastic disruptions.
library(plotly)
# 3D Wiener Process
N=500 # Number of random walk steps
# Define the X, Y, and Z, coordinate displacements independently
# xdis = rnorm(N, 0 , 1)
# ydis = rnorm(N, 0 , 2)
# zdis = rnorm(N, 0 , 3)
# xdis = cumsum(xdis)
# ydis = cumsum(ydis)
# zdis = cumsum(zdis)
# To use simulated 3D MVN Distribution
dis <- mixtools::rmvnorm(N, mu=c(0,0,0),
sigma=matrix(c(1,0,0, 0,2,0, 0,0,3), ncol=3))
# aggregate the displacements to get the actual 3D Cartesian Coordinates
xdis = cumsum(dis[,1])
ydis = cumsum(dis[,2])
zdis = cumsum(dis[,3])
# add Poisson noise
at = rpois(N, 0.1)
for(i in c(1:N)) {
if(at[i] != 0) {
xdis[i] = xdis[i]*at[i]
ydis[i] = ydis[i]*at[i]
zdis[i] = ydis[i]*at[i]
}
}
# plot(xdis, ydis, type="l",
# main ="Brownian Motion in Two Dimension with Poisson Arrival Process",
# xlab="x displacement", ylab = "y displacement")
plot_ly(x=xdis, y=ydis, z=zdis, type="scatter3d", mode="markers+lines",
text=~c(1:N), hoverinfo='text',
marker=list(color='gray'), showlegend=F) %>%
# Extentuate the starting and ending points
add_markers(x=xdis[1], y=ydis[1], z=zdis[1], marker=list(size=20,color="green"),
text=paste0("Starting Node 1")) %>%
add_markers(x=xdis[N], y=ydis[N], z=zdis[N],
marker=list(size=20,color="red"), text=paste0("Ending Node ", N)) %>%
layout(title="3D Brownian Motion",
scene=list(xaxis=list(title="X displacement"),
yaxis=list(title="Y displacement"),
zaxis=list(title="Z displacement")))
# 1D Wiener Process
# dis = rnorm(N, 0, 1);
# at = rpois(N,1)
# for(i in 1:N) {
# if(at[i] != 0){
# dis[i]= dis[i]*at[i]
# }
# }
# dis = cumsum(dis)
# plot(dis, type= "l",
# main= "Brownian Motion in One Dimension with Poisson Arrival Process",
# xlab="time", ylab="displacement")
# ub = 20; lb = -20
# xdis = rnorm(N, 0 ,1)
# xdis1 = rep(1,N)
# xdis1[1] = xdis[1]
# for(i in c(1:(N-1))){
# if(xdis1[i] + xdis[i+1] > ub) { xdis1[i+1] <- ub }
# else if(xdis1[i] + xdis[i+1] < lb) { xdis[i+1] = lb }
# else { xdis1[i+1] = xdis1[i] + xdis[i+1] }
# }
#
# plot(xdis1, type="l",main="Brownian Motion with bound in 1-dim", xlab="displacement",ylab="time")
# Compute the row Euclidean differences
df <- data.frame(cbind(xdis, ydis, zdis))
rowEuclidDistance <- dist(df)
plot_ly(z=as.matrix(rowEuclidDistance), type="heatmap") %>%
layout(title="Heatmap of Euclidean Distances between Consecutive Steps of Wiener Process")