SOCR ≫

TCIU Website ≫

TCIU GitHub ≫

Spacekime Analytics (Time Complexity and Inferential Uncertainty)

Spacetime IID vs. Spacekime Sampling

1 IID Spacetime Sampling vs. Spacekime Sampling

For all natural spacetime processes, various population characteristics like the mean, variance, range, and quantiles can be estimated by collecting independent and identically distributed (IID) samples. These samples represent observed data that is traditionally used to obtain sample-driven estimates of the specific population characteristics via standard formulas like the sample arithmetic average, variance, range, quantiles, etc. The latter approximate their population counterparts and form the basis for classical parametric and non-parametric statistical inference.

Typically, reliable spacetime statistical inference is conditional on the distribution of the native process as well as a sample-size reflecting the characteristics of the phenomenon. We will demonstrate that spacekime analytics can be equally effective with measuring a single spacetime observation and having a reasonable estimates of the unobserved process kime-phases.

Without loss of generality, suppose we have a pair of cohorts \(A\) and \(B\) and we obtain a series of measurements \(\{X_{A,i}\}_{i=1}^{n_A}\) and \(\{X_{B,i}\}_{i=1}^{n_B}\), respectively. Obviously the relations between the cohorts could widely vary, from being samples of the same process, to being related or completely independent.

To allow us to examine the extreme cases of pairing (1) IID cohorts (\(A\) and \(B\)), and (2) independent but differently distributed cohorts (\(A\) and \(C\)). The latter case may be thought of as a split of the first cohort (\(A\)) into training (\(C\)) and testing (\(D\)) groups. This design allows us to compare the classical spacetime-derived population characteristics of cohort \(A\) to their spacekime-reconstructed counterparts obtained using a single random kime-magnitude observation from \(A\) and kime-phases estimates derived from cohorts \(B\), \(C\) or \(D\).

1.1 fMRI timeseries

The demonstration below is based on a functional magnetic resonance imaging (fMRI) data, which is a 4D hypervolume with intensities representing the blood oxygenation level dependence at a specific spacetime location \((x,y,z,t)\). For simplicity, we will only focus on two fixed spatial locations with varying intensity distributions.

library(EBImage)
require(brainR)
library(spatstat) 
library(ggplot2)
library(kSamples)
library(reshape2)
library(beanplot)
library(rstanarm)

fMRIURL <- "http://socr.umich.edu/HTML5/BrainViewer/data/fMRI_FilteredData_4D.nii.gz"
fMRIFile <- file.path(tempdir(), "fMRI_FilteredData_4D.nii.gz")
download.file(fMRIURL, dest=fMRIFile, quiet=TRUE)
fMRIVolume <- readNIfTI(fMRIFile, reorient=FALSE)
# dimensions: 64 x 64 x 21 x 180 ; 4mm x 4mm x 6mm x 3 sec 

fMRIVolDims <- dim(fMRIVolume); # fMRIVolDims
# time_dim <- fMRIVolDims[4]; time_dim ## 180

# 2. extract the time-corse of 1D mid-axial slice (3D) hypervolume
xA_fMRI_1D_x20_y20_z11 <- fMRIVolume[20, 20, 11, ]; # length(xA_fMRI_1D_x20_y20_z11)   #  180
# hist(xA_fMRI_1D_x20_y20_z11)
library(plotly)
plot_ly(x = ~xA_fMRI_1D_x20_y20_z11, type = "histogram") %>%
  layout(bargap=0.1)

xB_fMRI_1D_x30_y30_z13 <- fMRIVolume[30, 30, 13, ]; # length(xB_fMRI_1D_x30_y30_z13)   #  180
# hist(xB_fMRI_1D_x30_y30_z13)

# Now, combine your two 1D timeseries into one dataframe for joint hist plotting as densities.  
# First make a new column in each that will be 
# a variable to identify where they came from later.
xA_df <- as.data.frame(xA_fMRI_1D_x20_y20_z11)
colnames(xA_df) <- "value"; xA_df$cohort <- "xA"
xB_df <- as.data.frame(xB_fMRI_1D_x30_y30_z13)
colnames(xB_df) <- "value"; xB_df$cohort <- "xB"

# and combine into your new data frame vegLengths
xA_xB_df <- rbind(xA_df, xB_df)

# ggplot(xA_xB_df, aes(value, fill = cohort)) + 
#   geom_density(alpha = 0.5, size=1.2) +
#   theme(text = element_text(size=20)) +
#   xlim(c(10200, 12000))

density_xA <- density(xA_xB_df[ which(xA_xB_df$cohort=="xA"), ]$value)
density_xB <- density(xA_xB_df[ which(xA_xB_df$cohort=="xB"), ]$value)
df_xA <- as.data.frame(cbind(x=density_xA$x, y=density_xA$y))
df_xB <- as.data.frame(cbind(x=density_xB$x, y=density_xB$y))

plot_ly(df_xA, x = ~x, y = ~y, type = 'scatter', 
        mode = 'lines', name = "xA", fill = 'tozeroy') %>% 
  add_trace(x = ~df_xB$x, y = ~df_xB$y, type = 'scatter', 
        mode = 'lines', name = "xB", fill = 'tozeroy') %>% 
  layout(title="Cohort A and B Distributions",
         xaxis = list(title = 'value'), yaxis = list(title = 'Density'))

Clearly the intensities of cohorts \(A\) and \(B\) are independent and follow different distribution. We’ll split the first cohort (\(A\)) into training (\(C\)) and testing (\(D\)) subgroups. Then we will:

• transform all four cohorts into Fourier k-space,
• iteratively randomly sample single observations from cohort \(C\),
• reconstruct the data into spacetime using alternative kime-phase estimates derived from cohorts \(B\) and \(D\), and
• compute and compare the classical spacetime-derived population characteristics of cohort \(A\) to their counterparts obtained using a single \(C\) kime-radial measurements paired with \(B\) and \(D\) kime-phases.

# Generic function to Transform Data to/from k-space (Space/Fourier domain)
kSpaceTransform <- function(data, inverse = FALSE, reconPhases = NULL) {
  # ForwardFT (rawData, FALSE, NULL)
  # InverseFT(magnitudes, TRUE, reconPhasesToUse) or InverseFT(FT_data, TRUE, NULL)
  FT_data <- array(complex(), length(data))
  mag_FT_data <- array(complex(), length(data))
  phase_FT_data <- array(complex(), length(data))
  IFT_reconPhases_data <- array(complex(), length(data))

  if (inverse == FALSE | is.null(reconPhases)) {
      FT_data <- fft(data, inverse)
      X2 <- FT_data
      mag_FT_data <- sqrt(Re(X2)^2+Im(X2)^2) 
      phase_FT_data <- atan2(Im(X2), Re(X2)) 
  }
  else {  # for IFT synthesis using user-provided Phases, typically from kime-phase aggregators
      Real <- data * cos(reconPhases)  
      Imaginary <- data * sin(reconPhases) 
      IFT_reconPhases_data <- 
          Re(fft(Real+1i*Imaginary, inverse = TRUE)/length(data))
  }

    ######### Test the FT-IFT analysis-synthesis back-and-forth transform process 
    #         to confirm calculations
    # X2 <- FT_data[ , 1]; mag_FT_data[ , 1] <- sqrt(Re(X2)^2+Im(X2)^2); 
    # phase_FT_data[ , 1] <- atan2(Im(X2), Re(X2)); 
    # Real2 = mag_FT_data[ , 1] * cos(phase_FT_data[ , 1])
    # Imaginary2 = mag_FT_data[ , 1] * sin(phase_FT_data[ , 1])
    # man_hat_X2 = Re(fft(Real2 + 1i*Imaginary2, inverse = T)/length(X2))
    # ifelse(abs(man_hat_X2[5] - data[5, 1]) < 0.001, "Perfect Syntesis", "Problems!!!")
    #########
  
    if (inverse == FALSE | is.null(reconPhases)) {
      return(list("magnitudes"=mag_FT_data, "phases"=phase_FT_data))
      # Use kSpaceTransform$magnitudes & kSpaceTransform$phases to retrieve teh Mags and Phases
    }
    else {
      return(IFT_reconPhases_data)
      # Use Re(kSpaceTransform) to extract spacetime Real-valued reconstructed data
    }
}

# 1. Split the first cohort ($A$) into *training* ($C$) and *testing* ($D$) subgroups.
subset_int <- sample(length(xA_df$value),floor(length(xA_df$value)*0.8))  
# 80% training + 20% testing
xC_fMRI_train <- xA_df$value [subset_int]; # length(xC_fMRI_train) # 144
xD_test <- xA_df$value [-subset_int]; # length(xD_test)  # 36

# 2. Transform all four cohorts into Fourier k-space
# xA, xB, xC_fMRI_train; xD_test
xA <- xA_fMRI_1D_x20_y20_z11; # length(xA)   #  180
xB <- xB_fMRI_1D_x30_y30_z13; # length(xB)   #  180

ft_xA <- fft(xA); ft_xB <- fft(xB)
ft_xC_fMRI_train <- fft(xC_fMRI_train); ft_xD_test <- fft(xD_test); 

# Magnitudes and Phases: Phase  <- atan(Im(img_ff)/Re(img_ff))
mag_ft_xA <- sqrt(Re(ft_xA)^2+Im(ft_xA)^2)
mag_ft_xB <- sqrt(Re(ft_xB)^2+Im(ft_xB)^2)
mag_ft_xC_fMRI_train <- sqrt(Re(ft_xC_fMRI_train)^2+Im(ft_xC_fMRI_train)^2)
mag_ft_xD_test <- sqrt(Re(ft_xD_test)^2+Im(ft_xD_test)^2)

phase_ft_xA <- atan2(Im(ft_xA), Re(ft_xA))
phase_ft_xB <- atan2(Im(ft_xB), Re(ft_xB))
phase_ft_xC_fMRI_train <- atan2(Im(ft_xC_fMRI_train), Re(ft_xC_fMRI_train))
phase_ft_xD_test <- atan2(Im(ft_xD_test), Re(ft_xD_test))

# Double-Check FT-IFT==I   ImplicitlyInvert the FT (IFT)
fftinv <- function( x ) { fft( x, inverse=TRUE ) / length( x ) }
# head(Re(fftinv(ft_xA))); head(xA)

# 3. Iteratively randomly sample single observations from cohort $C$,
N <- 30 # to 30 simulations
# take a random sample of size N (without replacement) from $C$
N_sampleIndx <- sample(1:length(xC_fMRI_train), N, replace=FALSE)
xC_fMRI_sampleN <- xC_fMRI_train[N_sampleIndx]
ft_xC_fMRI_sampleN_mag <- mag_ft_xC_fMRI_train[N_sampleIndx]

# 4. reconstruct the $C$ data into spacetime using a single ft_xC_fMRI_sampleN_mag value and alternative kime-phase estimates derived from cohorts $B$ and $D$ 
# for each ft_xC_fMRI_sampleN_mag[i] value, use $B$ and $D$ phases to reconstruct ift_ft_xC_fMRI_sampleN_PhaseB ift_ft_xC_fMRI_sampleN_PhaseD
ift_ft_xC_fMRI_1sampleN_PhaseB <- 
      array(dim=c(length(xC_fMRI_train), length(xC_fMRI_sampleN)))
ift_ft_xC_fMRI_1sampleN_PhaseD <-
      array(dim=c(length(xC_fMRI_train),length(xC_fMRI_sampleN)))
ift_ft_xC_fMRI_1sampleN_PhaseC <-
      array(dim=c(length(xC_fMRI_train),length(xC_fMRI_sampleN)))
ift_ft_xC_fMRI <- array(dim=length(xC_fMRI_train))
# dim(ift_ft_xC_fMRI_1sampleN_PhaseB) # [1] Time=144 Samples_N=30

for (i in 1:N) {
  ift_ft_xC_fMRI_1sampleN_PhaseB[ , i] <- 
    Re(kSpaceTransform(rep(ft_xC_fMRI_sampleN_mag[i], length(xC_fMRI_train)),
                       TRUE, phase_ft_xB[1:length(xC_fMRI_train)]))
  ift_ft_xC_fMRI_1sampleN_PhaseD[ , i] <- 
    Re(kSpaceTransform(rep(ft_xC_fMRI_sampleN_mag[i],
                           length(phase_ft_xD_test)), TRUE,
                       phase_ft_xD_test[1:length(phase_ft_xD_test)]))
  ift_ft_xC_fMRI_1sampleN_PhaseC[ , i] <- 
    Re(kSpaceTransform(rep(ft_xC_fMRI_sampleN_mag[i], length(xC_fMRI_train)),
                       TRUE, phase_ft_xC_fMRI_train[1:length(xC_fMRI_train)]))
}
ift_ft_xC_fMRI <- Re(kSpaceTransform(mag_ft_xC_fMRI_train, TRUE,
                       phase_ft_xC_fMRI_train[1:length(xC_fMRI_train)]))
# head(xC_fMRI_train) == head(ift_ft_xC_fMRI)

# 5.  compute and compare the *classical spacetime-derived* population characteristics of cohort $A$ to their counterparts obtained using a single $C$ kime-radial measurements paired with $B$ and $D$ kime-phases.
# Data = xC_fMRI_train, ift_ft_xC_fMRI_1sampleN_PhaseB, ift_ft_xC_fMRI_1sampleN_PhaseD
# length(xC_fMRI_train) == length(ift_ft_xC_fMRI_1sampleN_PhaseB[ , 1])

summary(scale(xC_fMRI_train))

##        V1          
##  Min.   :-3.23209  
##  1st Qu.:-0.67117  
##  Median : 0.01454  
##  Mean   : 0.00000  
##  3rd Qu.: 0.63308  
##  Max.   : 3.48507

summary(scale(ift_ft_xC_fMRI_1sampleN_PhaseC[ , 15]))

##        V1          
##  Min.   :-2.85429  
##  1st Qu.:-0.70335  
##  Median :-0.02874  
##  Mean   : 0.00000  
##  3rd Qu.: 0.66837  
##  Max.   : 3.07460

summary(scale(ift_ft_xC_fMRI_1sampleN_PhaseB[ , 15]))

##        V1          
##  Min.   :-3.32577  
##  1st Qu.:-0.56299  
##  Median :-0.03926  
##  Mean   : 0.00000  
##  3rd Qu.: 0.60848  
##  Max.   : 2.56786

summary(scale(ift_ft_xC_fMRI_1sampleN_PhaseD[ , 15]))

##        V1          
##  Min.   :-2.09096  
##  1st Qu.:-0.62007  
##  Median :-0.04813  
##  Mean   : 0.00000  
##  3rd Qu.: 0.59282  
##  Max.   : 2.02428

1.1.1 Data Visualization

# Plot all histograms as densities
ift_ft_xC_fMRI_1sampleN_PhaseC_df <- 
  as.data.frame(scale(ift_ft_xC_fMRI_1sampleN_PhaseC))
# colnames(xA_df) <- "value"; xA_df$cohort <- "xA"
#  xA_scale_df <- as.data.frame(scale(xA_df$value))
#  colnames(xA_scale_df) <- "value"; xA_scale_df$cohort <- "xA"
xB_df <- as.data.frame(scale(xB_fMRI_1D_x30_y30_z13))
colnames(xB_df) <- "value"; xB_df$cohort <- "xB"

# 
# and combine into your new data frame Lengths
xA_xB_df <- rbind(xA_df, xB_df)

# ggplot(xA_xB_df, aes(value, fill = cohort)) + 
#   geom_density(alpha = 0.5, size=1.2) +
#   theme(text = element_text(size=20)) +
#   xlim(c(10200, 12000))

density_xA <- density(xA_df$value)
df_xA <- as.data.frame(cbind(x=density_xA$x, y=density_xA$y))

plot_ly(df_xA, x = ~x, y = ~y, type = 'scatter', 
        mode = 'lines', name = "xA", fill = 'tozeroy') %>% 
  layout(title="Cohort A Distribution",
         xaxis = list(title = 'value'), yaxis = list(title = 'Density'))

# length(xC_fMRI_train); dim(ift_ft_xC_fMRI_1sampleN_PhaseB)
# dim(ift_ft_xC_fMRI_1sampleN_PhaseC); dim(ift_ft_xC_fMRI_1sampleN_PhaseD)

# Compute the averages accross all N=30 experiments for the B, C & D reconstructions
ift_ft_xC_fMRI_1sampleN_PhaseB_avg <- apply(ift_ft_xC_fMRI_1sampleN_PhaseB, 1, mean)
ift_ft_xC_fMRI_1sampleN_PhaseC_avg <- apply(ift_ft_xC_fMRI_1sampleN_PhaseC, 1, mean)
ift_ft_xC_fMRI_1sampleN_PhaseD_avg <- apply(ift_ft_xC_fMRI_1sampleN_PhaseD, 1, mean)

# Plot 4 density curves (orig=xC_fMRI and 3 reconstructions from B, C and D)
xC_fMRI_train_scale_df <- as.data.frame(scale(xC_fMRI_train))
colnames(xC_fMRI_train_scale_df) <- "value"; xC_fMRI_train_scale_df$series <- "xC_fMRI_original"

ift_ft_xC_fMRI_1sampleN_PhaseB_avg_scale_df <- as.data.frame(scale(ift_ft_xC_fMRI_1sampleN_PhaseB_avg))
colnames(ift_ft_xC_fMRI_1sampleN_PhaseB_avg_scale_df) <- "value"
ift_ft_xC_fMRI_1sampleN_PhaseB_avg_scale_df$series <- "SK_PhaseB"

ift_ft_xC_fMRI_1sampleN_PhaseC_avg_scale_df <- as.data.frame(scale(ift_ft_xC_fMRI_1sampleN_PhaseC_avg))
colnames(ift_ft_xC_fMRI_1sampleN_PhaseC_avg_scale_df) <- "value"
ift_ft_xC_fMRI_1sampleN_PhaseC_avg_scale_df$series <- "SK_PhaseC"

ift_ft_xC_fMRI_1sampleN_PhaseD_avg_scale_df <- as.data.frame(scale(ift_ft_xC_fMRI_1sampleN_PhaseD_avg))
colnames(ift_ft_xC_fMRI_1sampleN_PhaseD_avg_scale_df) <- "value"
ift_ft_xC_fMRI_1sampleN_PhaseD_avg_scale_df$series <- "SK_PhaseD"

# and combine into your new data frame vegLengths
xC_fMRI_SK_Phases_B_C_D_df <- rbind(xC_fMRI_train_scale_df, ift_ft_xC_fMRI_1sampleN_PhaseB_avg_scale_df, 
             ift_ft_xC_fMRI_1sampleN_PhaseC_avg_scale_df, ift_ft_xC_fMRI_1sampleN_PhaseD_avg_scale_df)

# library(ggplot2)
# ggplot(xC_fMRI_SK_Phases_B_C_D_df, aes(value, fill = series)) + 
#   geom_density(aes(color=series, linetype = series), alpha=0.4, size=1.2) +  # position = "stack"
#   theme(text = element_text(size=20)) +
#   scale_fill_manual( values = c("yellow", "red", "blue", "green")) +
#   geom_line(data=xC_fMRI_train_scale_df, stat = "density", color="purple", lty=4, lwd=2) +
#   ## guides(color = guide_legend(order=1)) +
#   theme(axis.title.x=element_blank(),axis.text.x=element_blank(), axis.ticks.x=element_blank())
#   # theme(legend.position="bottom")

density_xC_fMRI_orig <- density(xC_fMRI_SK_Phases_B_C_D_df[ 
  which(xC_fMRI_SK_Phases_B_C_D_df$series=="xC_fMRI_original"), ]$value)
density_SK_PhaseB <- density(xC_fMRI_SK_Phases_B_C_D_df[ 
  which(xC_fMRI_SK_Phases_B_C_D_df$series=="SK_PhaseB"), ]$value)
density_SK_PhaseC <- density(xC_fMRI_SK_Phases_B_C_D_df[ 
  which(xC_fMRI_SK_Phases_B_C_D_df$series=="SK_PhaseC"), ]$value)
density_SK_PhaseD <- density(xC_fMRI_SK_Phases_B_C_D_df[ 
  which(xC_fMRI_SK_Phases_B_C_D_df$series=="SK_PhaseD"), ]$value)

df_xC_fMRI_orig <- as.data.frame(cbind(
  x=density_xC_fMRI_orig$x, y=density_xC_fMRI_orig$y))
df_SK_PhaseB <- as.data.frame(cbind(
  x=density_SK_PhaseB$x, y=density_SK_PhaseB$y))
df_SK_PhaseC <- as.data.frame(cbind(
  x=density_SK_PhaseC$x, y=density_SK_PhaseC$y))
df_SK_PhaseD <- as.data.frame(cbind(
  x=density_SK_PhaseD$x, y=density_SK_PhaseD$y))

plot_ly(df_xC_fMRI_orig, x = ~x, y = ~y, type = 'scatter', 
        mode = 'lines', name = "xC_fMRI_orig", fill = 'tozeroy') %>% 
  add_trace(x = ~df_SK_PhaseB$x, y = ~df_SK_PhaseB$y, type = 'scatter', 
        mode = 'lines', name = "SK_PhaseB", fill = 'tozeroy') %>% 
  add_trace(x = ~df_SK_PhaseC$x, y = ~df_SK_PhaseC$y, type = 'scatter', 
        mode = 'lines', name = "SK_PhaseC", fill = 'tozeroy') %>% 
  add_trace(x = ~df_SK_PhaseD$x, y = ~df_SK_PhaseD$y, type = 'scatter', 
        mode = 'lines', name = "SK_PhaseD", fill = 'tozeroy') %>% 
  layout(title="Cohort Distribuitions - xC_fMRI_orig, SK_PhaseB, SK_PhaseC and SK_PhaseD", 
         legend = list(orientation = 'h'),
         xaxis = list(title = 'value'), yaxis = list(title = 'Density'))

1.1.1.1 Violin Plot

# ggplot(xC_fMRI_SK_Phases_B_C_D_df,aes(x=series, y=value, fill=series)) + 
#   geom_violin(trim=FALSE) +
#   geom_boxplot(width=0.1) + 
#   theme_bw()

xC_fMRI_SK_Phases_B_C_D_df %>%
  plot_ly(x = ~series, y = ~value , split = ~series, type = 'violin',
    box = list(visible = T), meanline = list(visible = T)) %>%
  layout(xaxis = list(title = "series"), 
         yaxis = list(title = "density", zeroline = F))

1.1.1.2 Empirical CDF

# ggplot(xC_fMRI_SK_Phases_B_C_D_df,aes(x=value, color=series)) +
#   stat_ecdf(size = 0.5)

df <- dplyr::arrange(xC_fMRI_SK_Phases_B_C_D_df, value)

pl <- ggplot(df, aes(x=value, color=series)) +
  stat_ecdf(size = 0.5)

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

ggplotly(pl) 

1.1.1.3 Beanplot

boxplot(xC_fMRI_SK_Phases_B_C_D_df$value ~ xC_fMRI_SK_Phases_B_C_D_df$series)
beanplot(xC_fMRI_SK_Phases_B_C_D_df$value ~ xC_fMRI_SK_Phases_B_C_D_df$series,
         border = "#CAB2D6",
         col = c("#CAB2D6", "#33A02C", "#B2DF8A"),
         side="second",add = T)

1.2 IID Simulation

The following simulation example generates two mixture distribution random samples each of \(n=10,000\) observations, \(\{X_{A,i}\}_{i=1}^{n_A}\), where \(X_{A,i} = 0.3U_i + 0.7V_i\), \(U_i \sim N(0,1)\) and \(V_i \sim N(5,3)\), and \(\{X_{B,i}\}_{i=1}^{n_B}\), where \(X_{B,i} = 0.4P_i + 0.6Q_i\), \(P_i \sim N(20,20)\) and \(Q_i \sim N(100,30)\).

n=10000
mu1 <- 0; mu2 <- 5
sig1 <- 1; sig2 <- 3
weight <- 0.7   

mixedDistFunc <- function (n, weight, mu1, mu2, sig1, sig2) {
    set.seed(1234); U <- rnorm(n, mean=mu1, sd = sig1)
    set.seed(1234); V <- rnorm(n,mean=mu2, sd = sig2)
    # randomly choose U or V
    set.seed(1234); wght <- rbinom(n, size=1, prob=weight) 
    X <- U*(1 - wght) + V*wght 
}

xA <- mixedDistFunc(n=n, weight, mu1, mu2, sig1, sig2)
hist(xA, freq = F)

# length(xA)   #  n=10,000

xB <- mixedDistFunc(n=n, 0.6, 20, 100, 20, 30)
hist(xB, freq = F)

# length(xB) 

# Now, combine your two univariate sets into one dataframe for joint hist plotting as densities.  
# First make a new column in each that will be 
# a variable to identify where they came from later.
xA_df <- as.data.frame(xA); colnames(xA_df)<-"value"; xA_df$cohort<-"xA"
xB_df <- as.data.frame(xB); colnames(xB_df)<-"value"; xB_df$cohort<-"xB"

# and combine into your new data frame vegLengths
xA_xB_df <- rbind(xA_df, xB_df)

Figure 5.2

ggplot(xA_xB_df, aes(value, fill = cohort)) + 
  geom_density(alpha = 0.5, size=1.2) +
  theme(text = element_text(size=20))

Clearly the intensities of cohorts \(A\) and \(B\) are independent and follow different distributions. We’ll split the first cohort (\(A\)) into training (\(C\)) and testing (\(D\)) subgroups, and then:

• transform all four cohorts into Fourier k-space,
• iteratively randomly sample single observations from cohort \(C\),
• reconstruct the data into spacetime using a single kime-magnitude value and alternative kime-phase estimates derived from cohorts \(B\), \(C\), and \(D\), and
• compute and compare the classical spacetime-derived population characteristics of cohort \(A\) to their spacekime counterparts obtained using a single \(C\) kime-magnitude paired with \(B\), \(C\), or \(D\) kime-phases.

# Generic function to Transform 1D Data to/from k-space (Space/Fourier domain)
kSpaceTransform <- function(data, inverse = FALSE, reconPhases = NULL) {
  # ForwardFT (rawData, FALSE, NULL)
  # InverseFT(magnitudes, TRUE, reconPhasesToUse) or InverseFT(FT_data, TRUE, NULL)
  FT_data <- array(complex(), length(data))
  mag_FT_data <- array(complex(), length(data))
  phase_FT_data <- array(complex(), length(data))
  IFT_reconPhases_data <- array(complex(), length(data))

  if (inverse == FALSE | is.null(reconPhases)) {
      FT_data <- fft(data, inverse)
      X2 <- FT_data
      mag_FT_data <- sqrt(Re(X2)^2+Im(X2)^2) 
      phase_FT_data <- atan2(Im(X2), Re(X2)) 
  }
  else {  # for IFT synthesis using user-provided Phases, typically from kime-phase aggregators
      Real <- data * cos(reconPhases)  
      Imaginary <- data * sin(reconPhases) 
      IFT_reconPhases_data <- 
          Re(fft(Real+1i*Imaginary, inverse = TRUE)/length(data))
  }

    ######### Test the FT-IFT analysis-synthesis back-and-forth transform process 
    #         to confirm calculations
    # X2 <- FT_data[ , 1]; mag_FT_data[ , 1] <- sqrt(Re(X2)^2+Im(X2)^2); 
    # phase_FT_data[ , 1] <- atan2(Im(X2), Re(X2)); 
    # Real2 = mag_FT_data[ , 1] * cos(phase_FT_data[ , 1])
    # Imaginary2 = mag_FT_data[ , 1] * sin(phase_FT_data[ , 1])
    # man_hat_X2 = Re(fft(Real2 + 1i*Imaginary2, inverse = T)/length(X2))
    # ifelse(abs(man_hat_X2[5] - data[5, 1]) < 0.001, "Perfect Syntesis", "Problems!!!")
    #########
  
    if (inverse == FALSE | is.null(reconPhases)) {
      return(list("magnitudes"=mag_FT_data, "phases"=phase_FT_data))
      # Use kSpaceTransform$magnitudes & kSpaceTransform$phases to retrieve teh Mags and Phases
    }
    else {
      return(IFT_reconPhases_data)
      # Use Re(kSpaceTransform) to extract spacetime Real-valued reconstructed data
    }
}

# 1. Split the first cohort ($A$) into *training* ($C$) and *testing* ($D$) subgroups.
subset_int <- sample(length(xA_df$value),floor(length(xA_df$value)*0.8))  
# 80% training + 20% testing
xC <- xA_df$value [subset_int]; # length(xC) # 8000
xD <- xA_df$value [-subset_int]; # length(xD)  # 2000

# 2. Transform all four cohorts into Fourier k-space
# xA, xB, xC (train), xD (test)
ft_xA <- fft(xA); ft_xB <- fft(xB)
ft_xC <- fft(xC); ft_xD <- fft(xD)

# Magnitudes and Phases: Phase  <- atan(Im(img_ff)/Re(img_ff))
mag_ft_xA <- sqrt(Re(ft_xA)^2+Im(ft_xA)^2)
mag_ft_xB <- sqrt(Re(ft_xB)^2+Im(ft_xB)^2)
mag_ft_xC <- sqrt(Re(ft_xC)^2+Im(ft_xC)^2)
mag_ft_xD <- sqrt(Re(ft_xD)^2+Im(ft_xD)^2)

phase_ft_xA <- atan2(Im(ft_xA), Re(ft_xA))
phase_ft_xB <- atan2(Im(ft_xB), Re(ft_xB))
phase_ft_xC <- atan2(Im(ft_xC), Re(ft_xC))
phase_ft_xD <- atan2(Im(ft_xD), Re(ft_xD))

# Double-Check FT-IFT==I   ImplicitlyInvert the FT (IFT)
fftinv <- function( x ) { fft( x, inverse=TRUE ) / length( x ) }
# head(Re(fftinv(ft_xA))); head(xA)

# 3. Iteratively randomly sample single observations from cohort $C$,
N <- 30 #  30 simulations
# take a random sample of size N (without replacement) from $C$
set.seed(1234); N_sampleIndx <- sample(1:length(xC), N, replace=FALSE)
xC_sampleN <- xC[N_sampleIndx]
ft_xC_mag <- mag_ft_xC[N_sampleIndx]

# 4. reconstruct the $C$ data into spacetime using a single ft_xC_sampleN_mag value and alternative kime-phase estimates derived from cohorts $B$ and $D$ 
# for each ft_xC_sampleN_mag[i] value, use $B$ and $D$ phases to reconstruct ift_ft_xC_sampleN_PhaseB ift_ft_xC_sampleN_PhaseD
ift_ft_xC_1sampleN_PhaseB <- 
      array(dim=c(length(xC), length(xC_sampleN)))
ift_ft_xC_1sampleN_PhaseD <-
      array(dim=c(length(xC), length(xC_sampleN)))
ift_ft_xC_1sampleN_PhaseC <-
      array(dim=c(length(xC), length(xC_sampleN)))
ift_ft_xC <- array(dim=length(xC))
# dim(ift_ft_xC_1sampleN_PhaseB) # [1] Size=8000 Samples_N=30

for (i in 1:N) {
  ift_ft_xC_1sampleN_PhaseB[ , i] <- 
    Re(kSpaceTransform(rep(ft_xC_mag[i], length(xC)),
                       TRUE, phase_ft_xB[1:length(xC)]))
  ift_ft_xC_1sampleN_PhaseD[ , i] <- 
    Re(kSpaceTransform(rep(ft_xC_mag[i], length(phase_ft_xD)), 
                       TRUE, phase_ft_xD[1:length(phase_ft_xD)]))
  ift_ft_xC_1sampleN_PhaseC[ , i] <- 
    Re(kSpaceTransform(rep(ft_xC_mag[i], length(xC)),
                       TRUE, phase_ft_xC[1:length(xC)]))
}
ift_ft_xC <- Re(kSpaceTransform(mag_ft_xC, TRUE,phase_ft_xC[1:length(xC)]))
# head(xC) - head(ift_ft_xC) # roundoff should be ~ 0

# 5.  compute and compare the *classical spacetime-derived* population characteristics of cohort $A$ to their counterparts obtained using a single $C$ kime-radial measurements paired with $B$ and $D$ kime-phases.
# Data = xC_train, ift_ft_xC_1sampleN_PhaseB, ift_ft_xC_1sampleN_PhaseD
# length(xC) == length(ift_ft_xC_1sampleN_PhaseB[ , 1])

summary(scale(xC))

##        V1          
##  Min.   :-2.38784  
##  1st Qu.:-0.88609  
##  Median :-0.03893  
##  Mean   : 0.00000  
##  3rd Qu.: 0.75821  
##  Max.   : 3.59925

summary(scale(ift_ft_xC_1sampleN_PhaseC[ , 15]))

##        V1          
##  Min.   :-2.52901  
##  1st Qu.:-0.76221  
##  Median :-0.05584  
##  Mean   : 0.00000  
##  3rd Qu.: 0.72999  
##  Max.   : 3.73114

summary(scale(ift_ft_xC_1sampleN_PhaseB[ , 15]))

##        V1           
##  Min.   :-3.798440  
##  1st Qu.:-0.636799  
##  Median : 0.009279  
##  Mean   : 0.000000  
##  3rd Qu.: 0.645119  
##  Max.   : 3.986702

summary(scale(ift_ft_xC_1sampleN_PhaseD[ , 15]))

##        V1          
##  Min.   :-2.66007  
##  1st Qu.:-0.79651  
##  Median :-0.08165  
##  Mean   : 0.00000  
##  3rd Qu.: 0.73477  
##  Max.   : 3.39448

1.2.1 Data Visualization

# Plot all histograms as densities
ift_ft_xC_1sampleN_PhaseC_df <- 
  as.data.frame(scale(ift_ft_xC_1sampleN_PhaseC))
# colnames(xA_df) <- "value"; xA_df$cohort <- "xA"
#  xA_scale_df <- as.data.frame(scale(xA_df$value))
#  colnames(xA_scale_df) <- "value"; xA_scale_df$cohort <- "xA"
xB_df <- as.data.frame(scale(xB))
colnames(xB_df) <- "value"; xB_df$cohort <- "xB"

# 
# and combine into your new data frame Lengths
xA_xB_df <- rbind(xA_df, xB_df)


ggplot(xA_xB_df, aes(value, fill = cohort)) + 
  geom_density(alpha = 0.5, size=1.2) +
  theme(text = element_text(size=20))