Let
\(f_X = \frac{1}{\sigma \sqrt{2\pi}} \cdot \exp{\left[-\frac{x^2}{2\sigma^2}\right]}\), then the maximum is attained at \(f_X(0) = \frac{1}{\sigma \sqrt{2\pi}}\).
Using a 217x217x11 cube we have for the \(X\) and \(Y\) dimensions let \(a = 217 / 4\) and for \(Z\) let \(a = 11 / 4\) and solve for \(\sigma_{X,Y}\) and \(\sigma_Z\):
\[ \frac{1}{\sigma \sqrt{2\pi}} \exp{\left[-\frac{a^2}{2\sigma^2}\right]} = \frac{1}{2\sigma \sqrt{2\pi}} \\ \Rightarrow a = \sigma \sqrt{2\ln{2}}, \] then \(\sigma_{X,Y}^2 \approx 2123.32\) and \(\sigma_{Z}^2 \approx 5.455\).
And we have \(f_{XYZ} = \phi(\cdot; \vec{0}, \Sigma)\) with \(\Sigma = \rm{diag}\left(2123.32, 2123.32, 5.455\right)\)
So the \(217 \times 217 \times 11\) array mask, \(A\), used to weight each cube will be defined as follows:
\(A_{i+109,j+109,k+6} = \int_i^{i+1}\int_j^{j+1}\int_k^{k+1} \phi(\cdot; \mu, \Sigma)dxdydz\) for \(i,j,k \in \{[-108,108] \times [-108,108] \times [-5,5]\}\)
f <- "collman14v2_fullCubes_20171101T1630.csv.h5"
h5ls(f)
loc <- h5read(f, name = "Locations")
colnames(loc) <- c("x", "y", "z")
chan <- h5read(f, name = "Channels")
cubes <- h5read(f, name = "collman14v2_cubes")
G <- readRDS("maskGaussian_mu0_sigma2123_2123_5-455.rds")
th <- theme(axis.text = element_blank(), axis.ticks = element_blank(),
axis.title.y = element_blank(), axis.title.x = element_blank(),
legend.position="bottom", legend.key.size = unit(1,'lines'),
panel.spacing = unit(0, "lines"))
png("meda_plots_gaussian/GaussianMask.png", width = 720, height = 140)
pdat <- melt(G)
ggplot(pdat,
aes(Var1,Var2, group = Var3, fill = value)) +
geom_raster() +
scale_y_reverse() +
facet_grid(. ~ Var3) +
#scale_fill_gradient(low = "black", high = "green") + th
scale_fill_gradient(low = "black", high = "white") + th
dev.off()
F0 <- matrix(NA, nrow = 1036, ncol = 11)
for(j in 1:dim(cubes)[5]){
for(i in 1:dim(cubes)[4]){
F0[i, j] <- sum(G * cubes[,,,i,j])
}
}
colnames(F0) <- chan
write.csv(F0, file = "collman14v2_GaussianCubes_20171101T1630.csv",
row.names = FALSE)
ccol <- c('blue', 'blue', 'blue', 'red', 'red',
'red', 'red', 'darkgreen', 'darkgreen', 'darkgreen',
'darkgreen')
sF0 <- scale(F0, center = TRUE, scale = TRUE)
set.seed(317)
Lg <- runAll(sF0, ccol = ccol)
Lg[[1]] <- mlocation(F0, ccol = ccol)
w = 720
h = 720
png("meda_plots_gaussian/d1heat.png", width = w, height = h)
p2 <- plot(Lg[[2]])
show(p2)
dev.off()
png("meda_plots_gaussian/mlocation.png", width = w, height = 0.5*h)
p1 <- plot(Lg[[1]])
show(p1)
dev.off()
png("meda_plots_gaussian/cumulativeVariance.png", width = w, height = h)
p3 <- plot(Lg[[3]])
show(p3)
dev.off()
png("meda_plots_gaussian/outliers.png", width = w, height = h)
p4 <- plot(Lg[[4]])
show(p4)
dev.off()
png("meda_plots_gaussian/cor.png", width = w, height = h)
p5 <- plot(Lg[[6]])
show(p5)
dev.off()
png("meda_plots_gaussian/pairhexGaussian.png", width = 2*w, height = 2*h)
pairhex(sF0)
dev.off()
png("meda_plots_gaussian/hmcClassificationsGaussian.png", width = 2*w, height = 2*h)
cr <- viridis(max(Lg[[7]]$dat$labels$col))
pairs(Lg[[7]]$dat$data, pch = 19, cex = 0.5, col = cr[Lg[[7]]$dat$labels$col])
dev.off()
png("meda_plots_gaussian/dendrograms.png", width = w, height = h)
p8 <- plotDend(Lg[[7]])
show(p8)
dev.off()
png("meda_plots_gaussian/stackMeans.png", width = w, height = 2*h)
p9 <- stackM(Lg[[7]], ccol = ccol, depth = 3, centered = TRUE)
show(p9)
dev.off()
png("meda_plots_gaussian/clusterMeans.png", width = w, height = h)
p10 <- clusterMeans(Lg[[7]], ccol = ccol)
show(p10)
dev.off()