R: Wie kann ich aus dem geometrischen Raum außerhalb einer 4D-Hyperellipse abtasten?

Question

Originaltitel (unpräzise): Wie finden Sie alle Vektoren, die eine 4-unbekannte nichtlineare Ungleichungen erfüllen?

---

Die Frage ist, einen Raum C zu finden, so dass jede x 4x1 Vektor in C erfüllt:

t(x) %*% t(Q) %*% Q %*% x > a, 

in denen Q eine 100 x 4 Matrix ich schon weiß, und a eine positive Konstante ist .

Ich habe versucht, die Lösung von Paketen wie ellipse, rootSolve und bvpSolve zu finden. Aber ich kann keine geeignete Lösung finden.

Jede Idee oder Lösung wird herzlich willkommen geheißen.

Answer 1

Anmerkung: Der unten stehende Algorithmus kann verwendet werden, um von der Oberfläche/dem Verteiler einer 4D-Hyperellipse oder von deren innerem Raum zu tasten.

Ich habe den Titel Ihrer Frage geändert. Es ist unmöglich, alle Lösungen aufzulisten, obwohl dieser Raum eine einfache mathematische Darstellung hat. Wir können unser Bestes aus diesem Raum probieren.

---

Transformation von Ellipse

hier Sphäre einige Mathematik basiert auf Cholesky-Faktorisierung. Alternativ betrachten wir symmetrische Eigenzerlegung, und ich habe eine Demonstration/Vergleich zwischen diesen beiden bei Obtain vertices of the ellipse on an ellipse covariance plot (created by car::ellipse), mit schönen Zahlen für die Geometrie.

Da Q bekannt ist, R ist bekannt. Der folgende Code wird R R:

R <- chol(crossprod(Q)) 

y aus der Hypersphäre mit dem Radius größer ist als sqrt(a). Wenn wir y aus einem solchen Raum probieren können, können wir es dann Karte zu x durch eine Dreieckssystem lösen:

x <- backsolve(R, y) 

---

Probenahme von y

Wir n-sphere coordinates verwenden können, so Raum parametrisieren . Für 4D Raum, haben wir:

die folgenden R-Funktion Proben n y-Vektoren, die von solchen Raum. Wegen der endlichen Repräsentation von Gleitkommazahlen können wir nicht den Unendlichkeitsradius aber .Machine$double.xmax bestenfalls haben. Wir verwenden aber auch ein optionales Argument rmax, wenn wir einen eingeschränkteren Radius wünschen.

ry <- function (n, rmin, rmax = NA) { 
    if (is.na(rmax)) rmax <- .Machine$double.xmax 
    if (rmin > rmax) stop("larger `rmax` expected!") 
    r <- runif(n, rmin, rmax) 
    phi1 <- runif(n, 0, pi) 
    phi2 <- runif(n, 0, 2 * pi) 
    phi3 <- runif(n, 0, 2 * pi) 
    matrix(c(r * cos(phi1), 
      r * sin(phi1) * cos(phi2), 
      r * sin(phi1) * sin(phi2) * cos(phi3), 
      r * sin(phi1) * sin(phi2) * sin(phi3)), 
     nrow = 4L, byrow = TRUE, dimnames = list(paste0("y", 1:4), NULL)) 
    } 

einige Beispiele Versuchen:

## radius between 4 and 10 
set.seed(0); ry(5, 4, 10) 
#   [,1]  [,2]  [,3]  [,4]  [,5] 
#y1 7.5594886 -5.31049687 -6.1388372 -3.5991830 -3.728597 
#y2 5.1402945 0.47936481 0.4799181 -2.5085948 -6.480402 
#y3 0.2614002 -1.68833263 -0.1950092 -5.9975328 -4.213166 
#y4 -2.0859078 0.02440839 -0.9452077 0.3052708 3.954674 

## radius between 4 and "inf" 
set.seed(0); ry(5, 4) 
#    [,1]   [,2]   [,3]   [,4]   [,5] 
#y1 1.299100e+308 -4.531902e+307 -6.588856e+307 -4.983772e+307 -6.442420e+307 
#y2 8.833607e+307 4.090831e+306 5.150993e+306 -3.473640e+307 -1.119710e+308 
#y3 4.492167e+306 -1.440799e+307 -2.093047e+306 -8.304756e+307 -7.279678e+307 
#y4 -3.584637e+307 2.082977e+305 -1.014498e+307 4.227070e+306 6.833046e+307 

I gewählt als jede Spalte jeder Zeile als eine Probe zu verwenden, die später Matrixberechnung zu erleichtern.

---

Transforming y zu x

Wir nehmen nun an haben

set.seed(0); Q <- matrix(runif(10 * 4), 10L, 4L) 

Wir bekommen R

R <- chol(crossprod(Q)) 
#   [,1]  [,2]  [,3]  [,4] 
#[1,] 2.176848 1.420882 1.2517326 1.4481875 
#[2,] 0.000000 1.077816 0.1045581 0.4646328 
#[3,] 0.000000 0.000000 1.2284251 0.3961126 
#[4,] 0.000000 0.000000 0.0000000 0.9019771 

SieAngenommen haben, dann bilden wir y-x:

a <- 4 
set.seed(0); y <- ry(5, sqrt(a), 10) ## we set an `rmax` here 
x <- backsolve(R, y) ## each column is a sample of `x` 
#   [,1]  [,2]  [,3]  [,4]  [,5] 
#[1,] 0.7403534 -1.49866534 -2.2359350 2.0269516 2.948561 
#[2,] 5.5481682 0.41827601 0.7024109 -1.7606720 -7.288339 
#[3,] 0.9373905 -1.01984708 0.1430660 -4.4180688 -4.749419 
#[4,] -2.2616584 0.01995357 -0.8367956 0.2995693 4.299268 

---

prüfen können wir prüfen, ob die oben x unsere Anforderung nicht erfüllt.

z <- Q %*% x 
ax <- colSums(x^2) ## value of `diag(x'Q'Qx)` 
#[1] 84.15453 17.00795 24.77044 43.33361 85.85250 

Sie sind alle größer als 4