## Examples of Probability Estimation on Manifold via Diffusion
("parameterized-manifold.R")
("diffusion.R")
("ridge.R")
("data.R")
(data.table)
(rgl)
## Circle with doubled density on the upper half -----------------------------------------------
## Data generation
N <- (train = 2000L, validate = 2000L)
DT <- data.table( = (((N)), times = N))
avgRelErr <- 0.03
(420L)
DT, s := (.N)]
DT, relErr := (.N) / $absZ * avgRelErr]
DT, x := (3 * * s) * (1 + relErr)]
DT, y := (3 * * s) * (1 + relErr)]
DTcircle <- data.table::data.table(s = (0, 300) / 300)
DTcircle, x := (2 * * s)]
DTcircle, y := (2 * * s)]
## Optimal Parzen window
X <- (DT[type == "train", .(x, y)])
Y <- (DT[type == "validate", .(x, y)])
DTparzen <- data.table( = avgRelErr * (0.2, 10, = 0.2))
DTparzen, anll := purrr::map_dbl(, , Y = Y, X = X)]
## DTparzen[anll < 0.3, plot(sigma, anll, type = 'b')]
## DTparzen[, abline(h = min(anll), col = "red")]
(, avgRelErr * (0.2, 10), Y = Y, X = X) # sigma = 0.025; ANLL = -0.003340218
## Compare with Silverman's Rule of Thumb for 1d Gaussian distributions.
h <- ((X), (X), (X, 2, ))
("Oversmoothing factors:", h / DTparzen[which.min(anll), ], "\n") # 8.3, 8.0
## Ridge finding
X <- (X)
outX <- (X)
Xs <- X, !outX]
Xs <- X, -(X)]
(llog <- (Xs, Xs, d = 1, h = 0.025, minCos = 0.05)) #(2k, 6s)
M <- l$M-(l$M), ]
## M <- M[, -outlier(M)]
(l1 <- (M, X, mult = 1)) #(4k, 0.7s)
<- 0.1
(l2 <- (M, l1$p, = , mult = 60, = 4)) #(4k, 31s)
DTmesh <- data.table::data.table(x = M[1, ], y = M[2, ], pEps = l1$p, fEps = l1$f,
= l2$, ft = l2$ft, pInf = l2$peron, fInf = l2$fPeron)
DTmesh, s := (y > 0, acos(x / (x^2 + y^2)), 2 * - acos(x / (x^2 + y^2))) / (2*)]
## Visualizations
(l$updatedPoints, names.arg = (l$updatedPoints),
xlab = "Iteration", ylab = "Updated pionts", main = "Ridge Estimation")
## Mesh plot
svglite::svglite("circle-scms.svg", = 10, = 10)
DT[type == "train", (x, y, = 1, = 19, = ("red", .2),
main = "Circle with doubled density on the upper half")]
DT[type == "validate", (x, y, = 19, = ("blue", .2))]
## DT[, rgl::plot3d(x, y, z, asp = FALSE)]
## DTwire[, rgl::lines3d(x, y, z, col = "red", lwd = 3)]
(M[1,], M[2,], = 19, = ("blue", .3))
(DT$x, DT$y, l$M[1,], l$M[2,], = (0, .5))
DTcircle, (x, y, = 'l', lwd = 2, = "white")]
(h = 0, lty = 2)
("bottomright", = ("data", "mesh", "manifold"),
= ("red", "blue", (.8)), = 19)
()
## Viz: CDF
svglite::svglite("circle-CDF.svg", = 10, = 10)
DTmesh[order(s)], (s, (pInf), = 's', lwd = 3, = "red",
= (0,1), = (0,1), xaxs = 'i', yaxs = 'i',
main = "CDF along arc length")]
(0, 0, 1/2, 2/3); (1, 1, 1/2, 2/3)
DTmesh[order(s)], (s[1], pInf[1], s[.N], 1)]
DTmesh[order(s)], (s, .I/.N, = 's', lwd = 1)]
DTmesh[order(s)], (s, (pEps), = 's', lwd = 1, = "magenta")]
DTmesh[order(s)], (s, (), = 's', lwd = 2, = "blue")]
("bottomright", = ("t = ", ("epsilon", , "Inf")),
= ("magenta", "blue", "red"), lty = 1, lwd = 2)
()
## Viz: PDF
svglite::svglite("circle-PDF.svg", = 10, = 10)
DTmesh, (s[order(s)], fEps[order(s)] / (fEps), = 'l', lwd = 1, = "magenta",
bty = 'n', main = "PDF along arc length (smoothed, normalized)")]
DTmesh, (s[order(s)], fInf[order(s)] / (fInf[order(s)]), lwd = 3, = "red")]
DTmesh, (s[order(s)], ft[order(s)] / (ft[order(s)]), lwd = 3, = "blue")]
(v = (0, 1/2, 1), h = (0, 1/2, 1), lty = 2)
("bottomright", = ("t = ", ("epsilon", , "Inf")),
= ("magenta", "blue", "red"), lty = 1, lwd = 2)
()
## Sphere with density sin(phi) -----------------------------------------------------------------
## Data generation
N <- 4000L
(420L)
## Sphere: sin(phi) on manifold
DT <- (N, (theta, phi) ( * phi))
DT, ("x", "y", "z") := (theta, phi)]
## Plots: parameter space vs. embedding space.
DT, (theta, phi, = 19, = (0, .3))]
(0,0,1,1, = (0, .2), border = )
DT, rgl::plot3d(x, y, z, aspect = ,
= (("red", "blue"))(100)[100 * phi + 1], alpha = .5)]
## Ambient noise
avgRelErr <- 0.03
<- avgRelErr / $absZ / (3 - 1)
(420L)
DT, ("x", "y", "z") := (.SD, (x) x + * (.N)), .SDcols = ("x", "y", "z")]
## Data plot
DT, rgl::plot3d(x, y, z, = (x^2+y^2+z^2 > 1, "red", "blue"), alpha = .5, aspect = )]
rgl::view3d(theta = 0, phi = 0, zoom = 0.6)
rgl::rgl.postscript("sphere.svg", "svg")
## 3D data plot
## htmlwidgets::saveWidget(rgl::rglwidget(width = 1000, height = 1000), "sphere-data.html")
## Computation: manifold estimation
X <- ((DT, .(x, y, z)]))
## (1) Data as initial mesh
## if (length(ol <- outlier(X, 6))) Xs <- X[, -ol]
(l <- (X, X, d = 2, epsilon = 0.2, maxIter = 30)) #(4k, 17s)
M <- l$M-(l$M), ]
(iqr <- ((((M[seq(3),]^2)), (1/4, 3/4))))
(mbias <- (((M[seq(3),]^2))) - 1)
## (IQR, median bias): (5.77e-3, 2.73%)
(l <- (X, X, d = 2, h = (4 * , 3), epsilon = 0.2, maxIter = 30)) #(4k, 20s)
## (mult, IQR, median bias): (5, 5.77e-3, 1.65%); (4, 6.76e-3, 1.0%); (3, 9.2e-3, 5.0e-3)
r <- ((M^2))
d <- (M)
(d, r, xlab = "minimum distance to set", ylab = "radius")
(h = (r) + (r) * (-3, 3, 4), v = (d) + (d) * (-3, 3))
## Remove outliers
M <- M, -(((r - (r)) / (r)) > 3)]
## (2) Mesh with uniform distribution on sphere
DTmesh <- (N, (theta, phi) 1)
DTmesh, ("x", "y", "z") := (theta, phi)]
M <- ((DTmesh, .(x, y, z)]))
## (3) Regular mesh on sphere
## TODO...
## Viz: 3D mesh plot
rgl::plot3d(M[1,], M[2,], M[3,], = "green", alpha = .3)
## htmlwidgets::saveWidget(rgl::rglwidget(width = 1000, height = 1000), "sphere-mesh.html")
## Viz: mesh radius
svglite::svglite("sphere-radius-ECDF.svg", = 8, = 8)
DT, .(r = (x^2+y^2+z^2))
][order(r)], (r, = (r < 1, "blue", "red"), lwd = 2, bty = 'n',
xlab = "radius", main = "ECDF of radius")]
(((M[seq(3),]^2)), add = , = "green", lwd = 2)
(v = 1, lty = 2)
("bottomright", = ("Data (outside sphere)", "Data (inside sphere)", "Mesh"),
= ("red", "blue", "green"), lty = 1, lwd = 2)
()
## Relative standard deviation
((((M[seq(3),]^2))) / DT, ((x^2+y^2+z^2))]) #0.19
## Computation: intial probability estimation
(l1 <- (M, X, mult = 1/2)) #(4k, 1s)
## Computation: probability diffusion
## (mult = 3, scale = 1); (mult = 2, scale = 10); smooth = 4
mult <- 2L
<- 10L
<- 3L
mp <- (M)
(l2 <- (M, l1$p, = , mult = mult, = , = )) #(4k, 30s)
DTmesh <- data.table::data.table(x = M[1, ], y = M[2, ], z = M[3, ], p0 = l1$p, f0 = l1$f,
= l2$, ft = l2$ft, pInf = l2$peron, fInf = l2$fPeron)
DTmesh, ("theta", "phi") := (x, y, z)[c("theta", "phi")]]
## Viz: PDF
rgl::par3d(userMatrix = ((1,4))[c(2,3,1,4),] + 1e-4)
DTmesh, rgl::plot3d(theta, phi, fInf, alpha = .3)]
rgl::open3d()
DTmesh, rgl::plot3d(theta, phi, ft, alpha = .3)]
## Viz: Comparison
svglite::svglite("sphere-PDF.svg", = 10, = 10)
DTmesh[theta > 4, (phi, ft, = 'n',
= (0, (ft)*1.02), yaxs = 'i', bty = 'n',
main = "PDF along zenith angle (smoothed)")]
DTmesh[theta > 4, (phi, f0 / (f0) * (ft), = 19, = (0, 0.2))]
DTmesh[theta > 4, (phi, ft, = 19, = ("blue", 0.2))]
DTmesh[theta > 4, (phi, fInf, = 19, = ("red", 0.2))]
(DTmesh[theta > 4 & phi > 0.4 * & phi < 0.6 * , (ft)] * (x),
0, , lwd = 2, = "blue", add = )
(h = DTmesh[theta > 4, (fInf)], lwd = 2, = "red")
("topright", = ("t = ", (0, ( * mult / mp, digits = 3), "Inf")),
= ("gray", "blue", "red"), = 19)
()
## Ghanem's two-circle data ------------------------------------------------------------
DT <- data.table::fread("two-circles.csv")
DT, (x, y, = 1, = 19, = (0, 0.3))]
## Density Ridge
X <- ((DT, .(x, y)]))
M0 <- X
outliers <- (X)
Xs <- ((outliers) == 0) X X, -outliers]
## Convergence criterion
(l <- (M0, Xs, d = 1)) #(231, 2.5s)
(l <- (M0, Xs, d = 1, epsilon = 0.055)) #(231, 1.8s)
(l <- (M0, Xs, d = 1, epsilon = 0.45)) #(231, 0.083s)
M <- l$M-(l$M),]
((M), = 19, = "blue")
## Density kernel bandwidth
## Reducing on Silverman's rule of thumb: h is too large; h * 0.1 is too small.
h <- ((X), (X), (X, 1, ))
(l <- (M0, Xs, d = 1, h = h * 0.2)) #(231, 2.5s)
M <- l$M-(l$M), ]
((M), = 19, = "blue")
((0, 2), (0, 2), circles = h, inches = , add = )
((0, 2), (0, 2), circles = h * 0.2, inches = , add = )
## Sampling
(42L)
mErr <- (((M0)), = 2) * h * 0.1
M1 <- M0 + mErr
(l <- (M1, Xs, d = 1, h = h * 0.2)) #(231, 2.5s)
M <- l$M-(l$M), ]
## Plot
svglite::svglite("twocircles-mssm.svg", = 24, = 12)
DT, (x, y, = 1, = 19, = (0, 0.3))]
((M1), = 19, = "blue")
((M), = 19, = "red")
(M1[1,], M1[2,], M[1,], M[2,])
()
## Single point perturbation
DTp <- DT[c(20, 115, 142, 230), .(x0 = x, y0 = y)]
(42L)
DTp <- DTp[rep((.N), each = 100)
], ("xErr", "yErr") := (2, (.N) * (h*0.166), simplify = )]]
DTp, x := x0 + xErr], y := y0 + yErr]]
M1 <- ((DTp, .(x, y)]))
(l <- (M1, Xs, d = 1, h = h * 0.4)) #(, s)
M <- l$M-(l$M), ]
svglite::svglite("twocircles-mssm-noise-contraction.svg", = 24, = 12)
(DTp, .(x, y)], DT, (x = M[1,], y = M[2,])), (x, y, = 1, = 'n')]
DT, (x, y, = 1, = 19, = (0, 0.3))]
(M1[1,], M1[2,], M[1,], M[2,], = "gray75")
DTp, (x, y, = 19, = ("red", alpha.f = 0.3))]
((M), = 19, = ("blue", alpha.f = 0.3))
DTp, (x0, y0, = 19, = "black")]
()
