Jak znaleźć macierz kowariancji wielokąta?

Wyobraź sobie, że masz wielokąt zdefiniowany przez zestaw współrzędnych $(x_1,y_1)...(x_n,y_n)$ a jego środek masy wynosi $(0,0)$. Można traktować wielokąt jako rozkład równomierny z granicą wielokąta.

Poszukuję metody, która znajdzie macierz kowariancji wielokąta .

Podejrzewam, że macierz kowariancji wielokąta jest ściśle związana z drugim momentem pola , ale nie jestem pewien, czy są one równoważne. Wzory znajdujące się w artykule w Wikipedii, który połączyłem, wydają się (przypuszczam, że nie jest to dla mnie szczególnie jasne z artykułu) odwoływać się do bezwładności obrotowej wokół osi x, y zamiast z głównych osi wielokąta.

(Nawiasem mówiąc, jeśli ktoś może wskazać mi, jak obliczyć główne osie wielokąta, byłoby to również przydatne)

Kuszące jest po prostu wykonanie PCA na współrzędnych , ale powoduje to, że współrzędne niekoniecznie są równomiernie rozmieszczone wokół wielokąta, a zatem nie są reprezentatywne dla gęstości wielokąta. Skrajnym przykładem jest zarys Północnej Dakoty, której wielokąt jest zdefiniowany przez dużą liczbę punktów wzdłuż rzeki Czerwonej, a także tylko dwa kolejne punkty określające zachodnią krawędź stanu.

Najpierw zróbmy analizę.

Załóżmy, że w obrębie wielokąta jego gęstość prawdopodobieństwa jest funkcją proporcjonalną Zatem stała proporcjonalności jest odwrotnością całki nad wielokątem,$\mathcal{P}$$p(x,y).$$p$

$$\mu_{0,0}(\mathcal{P})=\iint_{\mathcal P} p(x,y) \mathrm{d}x\,\mathrm{d}y.$$

Środka ciężkości wielokąta jest punktem średnich współrzędnych, obliczony jako pierwszych chwil. Pierwszy to

$$\mu_{1,0}(\mathcal{P})=\frac{1}{\mu_{0,0}(\mathcal{P})} \iint_{\mathcal P} x\,p(x,y)\mathrm{d}x\,\mathrm{d}y.$$

Bezwładnościowy napinacz może być przedstawiony jako symetryczny tablicy drugich momentów obliczane po przeliczeniu wielokąta umieścić jej środka ciężkości w punkcie początkowym, to znaczy w matrycy centralnych drugich momentów

$$\mu^\prime_{k,l}(\mathcal{P}) = \frac{1}{\mu_{0,0}(\mathcal{P})} \iint_{\mathcal P} \left(x - \mu_{1,0}(\mathcal{P})\right)^k\,\left(y - \mu_{0,1}(\mathcal{P})\right)^l\,p(x,y)\mathrm{d}x\,\mathrm{d}y$$

gdzie wynoszą od do do Sam tensor - inaczej macierz kowariancji - jest$(k,l)$$(2,0)$$(1,1)$$(0,2).$

$$I(\mathcal{P}) = \pmatrix{\mu^\prime_{2,0}(\mathcal{P}) & \mu^\prime_{1,1}(\mathcal{P}) \\ \mu^\prime_{1,1}(\mathcal{P}) & \mu^\prime_{0,2}(\mathcal{P})}.$$

PCA od otrzymuje się główne osie z są to wektory jednostkowe skalowane przez ich wartości własnych.$I(\mathcal{P})$$\mathcal{P}:$

---

Następnie sprawdźmy, jak wykonać obliczenia. Ponieważ wielokąt jest przedstawiany jako sekwencja wierzchołków opisujących jego zorientowaną granicę naturalne jest wywoływanie$\partial\mathcal P,$

Twierdzenie Greena: gdzie to jedna forma zdefiniowana w sąsiedztwie i

$$\iint_{\mathcal{P}} \mathrm{d}\omega = \oint_{\partial\mathcal{P}}\omega$$ $\omega = M(x,y)\mathrm{d}x + N(x,y)\mathrm{d}y$$\mathcal{P}$ $$\mathrm{d}\omega = \left(\frac{\partial}{\partial x}N(x,y) - \frac{\partial}{\partial y}M(x,y)\right)\mathrm{d}x\,\mathrm{d}y.$$

Na przykład, z i stałą ( tj. Jednolitą) gęstością możemy (przez inspekcję) wybrać jedną z wielu rozwiązania, takie jak$\mathrm{d}\omega = x^k y^l \mathrm{d}x\mathrm{d}y$$p,$

$$\omega(x,y) = \frac{-1}{l+1}x^k y^{l+1}\mathrm{d}x.$$

Chodzi o to, że całka konturu podąża za segmentami linii wyznaczonymi przez sekwencję wierzchołków. Każdy segment linii od wierzchołka do wierzchołka można sparametryzować za pomocą zmiennej rzeczywistej w postaci$\mathbf{u}$$\mathbf{v}$$t$

$$t \to \mathbf{u} + t\mathbf{w}$$

gdzie to normalny kierunek jednostki od doWartości wahają się zatem od do Pod tą parametryzacją i są liniowymi funkcjami i a są liniowymi funkcjami Tak więc podcałkową całki konturu na każdej krawędzi zostaje funkcja wielomianowa od , która jest łatwo ocenione dla małych i$\mathbf{w} \propto \mathbf{v}-\mathbf{u}$$\mathbf{u}$$\mathbf{v}.$$t$$0$$|\mathbf{v}-\mathbf{u}|.$$x$$y$$t$$\mathrm{d}x$$\mathrm{d}y$$\mathrm{d}t.$$t,$$k$$l.$

---

Wdrożenie tej analizy jest tak proste, jak kodowanie jej komponentów. Na najniższym poziomie potrzebujemy funkcji do zintegrowania wielomianowej formy jednoczęściowej na segmencie linii. Funkcje wyższego poziomu agregują je, aby obliczyć momenty surowe i centralne w celu uzyskania barycentrum i tensora bezwładnościowego, a na koniec możemy działać na tym tensorze, aby znaleźć główne osie (które są jego skalowanymi wektorami własnymi). Poniższy Rkod wykonuje tę pracę. Nie ma żadnych pretensji do wydajności: ma on jedynie zilustrować praktyczne zastosowanie powyższej analizy. Każda funkcja jest prosta, a konwencje nazewnictwa są zbieżne z konwencjami analizy.

Kod zawiera procedurę generowania prawidłowych zamkniętych, po prostu połączonych, nie przecinających się wielokątów (przez losowe deformowanie punktów wzdłuż koła i dołączenie początkowego wierzchołka jako jego ostatniego punktu w celu utworzenia zamkniętej pętli). Poniżej znajduje się kilka instrukcji do wykreślenia wielokąta, wyświetlenia jego wierzchołków, przylegania do centrum środka i wykreślenia głównych osi w kolorze czerwonym (największym) i niebieskim (najmniejszym), tworząc układ współrzędnych zorientowany dodatnio na wielokąt.

#
# Integrate a monomial one-form x^k*y^l*dx along the line segment given as an 
# origin, unit direction vector, and distance.
#
lintegrate <- function(k, l, origin, normal, distance) {
  # Binomial theorem expansion of (u + tw)^k
  expand <- function(k, u, w) {
    i <- seq_len(k+1)-1
    u^i * w^rev(i) * choose(k,i)
  }
  # Construction of the product of two polynomials times a constant.
  omega <- normal[1] * convolve(rev(expand(k, origin[1], normal[1])), 
                                expand(l, origin[2], normal[2]),
                                type="open")
  # Integrate the resulting polynomial from 0 to `distance`.
  sum(omega * distance^seq_along(omega) / seq_along(omega))
}
#
# Integrate monomials along a piecewise linear path given as a sequence of
# (x,y) vertices.
#
cintegrate <- function(xy, k, l) {
  n <- dim(xy)[1]-1 # Number of edges
  sum(sapply(1:n, function(i) {
    dv <- xy[i+1,] - xy[i,]               # The direction vector
    lambda <- sum(dv * dv)
    if (isTRUE(all.equal(lambda, 0.0))) {
      0.0
    } else {
      lambda <- sqrt(lambda)              # Length of the direction vector
      -lintegrate(k, l+1, xy[i,], dv/lambda, lambda) / (l+1)
    }
  }))
}
#
# Compute moments of inertia.
#
inertia <- function(xy) {
  mass <- cintegrate(xy, 0, 0)
  barycenter = c(cintegrate(xy, 1, 0), cintegrate(xy, 0, 1)) / mass
  uv <- t(t(xy) - barycenter)   # Recenter the polygon to obtain central moments
  i <- matrix(0.0, 2, 2)
  i[1,1] <- cintegrate(uv, 2, 0)
  i[1,2] <- i[2,1] <- cintegrate(uv, 1, 1)
  i[2,2] <- cintegrate(uv, 0, 2)
  list(Mass=mass,
       Barycenter=barycenter,
       Inertia=i / mass)
}
#
# Find principal axes of an inertial tensor.
#
principal.axes <- function(i.xy) {
  obj <- eigen(i.xy)
  t(t(obj$vectors) * obj$values)
}
#
# Construct a polygon.
#
circle <- t(sapply(seq(0, 2*pi, length.out=11), function(a) c(cos(a), sin(a))))
set.seed(17)
radii <- (1 + rgamma(dim(circle)[1]-1, 3, 3))
radii <- c(radii, radii[1])  # Closes the loop
xy <- circle * radii
#
# Compute principal axes.
#
i.xy <- inertia(xy)
axes <- principal.axes(i.xy$Inertia)
sign <- sign(det(axes))
#
# Plot barycenter and principal axes.
#
plot(xy, bty="n", xaxt="n", yaxt="n", asp=1, xlab="x", ylab="y",
     main="A random polygon\nand its principal axes", cex.main=0.75)
polygon(xy, col="#e0e0e080")
arrows(rep(i.xy$Barycenter[1], 2), 
       rep(i.xy$Barycenter[2], 2),
       -axes[1,] + i.xy$Barycenter[1],     # The -signs make the first axis .. 
       -axes[2,]*sign + i.xy$Barycenter[2],# .. point to the right or down.
       length=0.1, angle=15, col=c("#e02020", "#4040c0"), lwd=2)
points(matrix(i.xy$Barycenter, 1, 2), pch=21, bg="#404040")

Edycja: Nie zauważyłem, że whuber już odpowiedział. Zostawię to jako przykład innego (być może mniej eleganckiego) podejścia do problemu.

Macierz kowariancji

Niech za losową z rozkładu równomiernego na wieloboku z obszaru . Macierz kowariancji to:$(X,Y)$$P$$A$

$$C = \begin{bmatrix} C_{XX} & C_{XY} \\ C_{XY} & C_{YY} \end{bmatrix}$$

gdzie to wariancja , to wariancja , a to kowariancja między i . Zakłada to średnią zerową, ponieważ środek masy wielokąta znajduje się na początku. Rozkład równomierny przypisuje stałą gęstość prawdopodobieństwa do każdego punktu w , więc:$C_{XX} = E[X^2]$$X$$C_{YY} = E[Y^2]$$Y$$C_{XY} = E[XY]$$X$$Y$$\frac{1}{A}$$P$

$$C_{XX} = \frac{1}{A} \underset{P}{\iint} x^2 dV \quad C_{YY} = \frac{1}{A} \underset{P}{\iint} y^2 dV \quad C_{XY} = \frac{1}{A} \underset{P}{\iint} x y dV \tag{1}$$

Triangulacja

Zamiast próbować bezpośrednio zintegrować skomplikowany region, taki jak , możemy uprościć problem, dzieląc na trójkątnych podregionów:$P$$P$$n$

$$P = T_1 \cup \cdots \cup T_n$$

W twoim przykładzie jeden z możliwych partycjonowania wygląda następująco:

Istnieją różne sposoby uzyskania triangulacji (patrz tutaj ). Na przykład, możesz obliczyć triangulację wierzchołków Delaunaya , a następnie odrzucić krawędzie, które wypadają poza (ponieważ może nie być wypukłe jak w przykładzie).$P$

Całki powyżej można następnie podzielić na sumy całek w trójkątach:$P$

$$C_{XX} = \frac{1}{A} \sum_{i=1}^n \underset{T_i}{\iint} x^2 dV \quad C_{YY} = \frac{1}{A} \sum_{i=1}^n \underset{T_i}{\iint} y^2 dV \quad C_{XY} = \frac{1}{A} \sum_{i=1}^n \underset{T_i}{\iint} x y dV \tag{2}$$

Trójkąt ma ładne, proste granice, więc te całki są łatwiejsze do oceny.

Całkowanie przez trójkąty

Istnieją różne sposoby integracji nad trójkątami. W tym przypadku zastosowałem lewę, która polega na odwzorowaniu trójkąta na kwadrat jednostki. Przekształcenie w barycentryczne współrzędne może być lepszym rozwiązaniem.

Oto rozwiązania dla całek powyżej, dla dowolnego trójkąta zdefiniowanego przez wierzchołki . Pozwolić:$T$$(x_1,y_1), (x_2,y_2), (x_3,y_3)$

$$v_x = \left[ \begin{smallmatrix} x_1 \\ x_2 \\ x_3 \end{smallmatrix} \right] \quad v_y = \left[ \begin{smallmatrix} y_1 \\ y_2 \\ y_3 \end{smallmatrix} \right] \quad \vec{1} = \left[ \begin{smallmatrix} 1 \\ 1 \\ 1 \end{smallmatrix} \right] \quad L = \left[ \begin{smallmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 1 \end{smallmatrix} \right]$$

Następnie:

$$\underset{T}{\iint} x^2 dV = \frac{A}{6} \text{Tr}(v_x v_x^T L) \quad \underset{T}{\iint} y^2 dV = \frac{A}{6} \text{Tr}(v_y v_y^T L) \quad \underset{T}{\iint} x y dV = \frac{A}{12} (\vec{1}^T v_x v_y^T \vec{1} + v_x^T v_y) \tag{3}$$

Składając wszystko w całość

Niech i zawierają X / Y współrzędne wierzchołków każdego trójkąta , jak opisano powyżej. Podłącz do dla każdego trójkąta, zwracając uwagę, że warunki obszaru anulują się. To daje rozwiązanie:$v_x^i$$v_y^i$$T_i$$(3)$$(2)$

$$C_{XX} = \frac{1}{6} \sum_{i=1}^n \text{Tr} \big( v_x^i (v_x^i)^T L \big) \quad C_{YY} = \frac{1}{6} \sum_{i=1}^n \text{Tr} \big( v_y^i (v_y^i)^T L \big) \quad C_{XY} = \frac{1}{12} \sum_{i=1}^n \big( \vec{1}^T v_x^i (v_y^i)^T \vec{1} + (v_x^i)^T v_y^i \big) \tag{4}$$

Główne osie

Główne osie podane są przez wektory własne macierzy kowariancji , podobnie jak w PCA. W przeciwieństwie do PCA, mamy analityczne wyrażenie dla , zamiast konieczności szacowania go z próbkowanych punktów danych. Zauważ, że same wierzchołki nie są reprezentatywną próbką z jednorodnego rozkładu na , więc nie można po prostu pobrać przykładowej macierzy kowariancji wierzchołków. Ale * jest * stosunkowo prostą funkcją wierzchołków, jak widać w .$C$$C$$P$$C$$(4)$