Jak znaleźć macierz kowariancji wielokąta?

9

Wyobraź sobie, że masz wielokąt zdefiniowany przez zestaw współrzędnych (x1,y1)...(xn,yn) a jego środek masy wynosi (0,0). Można traktować wielokąt jako rozkład równomierny z granicą wielokąta. wprowadź opis zdjęcia tutaj

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.

Ingolifs
źródło
Zakładając „znajdź”, zakładam po prostu pobieranie próbek z wielokąta, a następnie obliczanie kowariancji próbek, czy nie masz na myśli?
Stephan Kolassa
Czy możesz również edytować swój post, aby uwzględnić współrzędne wielokąta, aby ludzie mogli się nim bawić?
Stephan Kolassa
1
@StephanKolassa Mam na myśli traktowanie wielokąta jako jednolitej dwuwymiarowej gęstości prawdopodobieństwa z granicą wielokąta. Jasne, możesz próbkować punkty, a limit byłby taki sam, ale szukam metody a priori. Obraz jest tylko ilustracją z farby, której użyłem. Dane ze świata rzeczywistego, których zamierzam użyć, to zarysy stanów i regionów.
Ingolifs,
1
Masz rację, że zwykłym terminem „macierz kowariancji” jest moment bezwładności lub drugi moment. Główne osie są zorientowane w swoich ośmiu kierunkach. Uruchomienie PCA na współrzędnych jest nieprawidłowe: jest równoznaczne z założeniem, że cała masa znajduje się w wierzchołkach. Najbardziej bezpośrednie metody obliczania barycenter - pierwsza chwila - są omówione w moim poście na gis.stackexchange.com/a/22744/664 . Drugi moment oblicza się w ten sam sposób z niewielkimi modyfikacjami. Kula wymaga specjalnych rozważań.
whuber
2
Działa to w drugą stronę: oblicz tensor bezwładności i znajdź na nim jego główne osie. Technika w twoim przypadku obejmuje Twierdzenie Greena, które pokazuje, że wymagane całki
μk,l(P)=Pxkyldxdy
można obliczyć całek konturu wokół w jednej postaci gdzieTakie postacie są łatwe do znalezienia, ponieważ każda odpowiednia kombinacja liniowa i zadziała. Całka konturu jest sumą całek na krawędziach. Pωdω=xkyldxdy.xkyl+1dxxk+1yldy
whuber

Odpowiedzi:

10

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,Pp(x,y).p

μ0,0(P)=Pp(x,y)dxdy.

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

μ1,0(P)=1μ0,0(P)Pxp(x,y)dxdy.

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

μk,l(P)=1μ0,0(P)P(xμ1,0(P))k(yμ0,1(P))lp(x,y)dxdy

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

I(P)=(μ2,0(P)μ1,1(P)μ1,1(P)μ0,2(P)).

PCA od otrzymuje się główne osie z są to wektory jednostkowe skalowane przez ich wartości własnych.I(P)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ływanieP,

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

Pdω=Pω
ω=M(x,y)dx+N(x,y)dyP
dω=(xN(x,y)yM(x,y))dxdy.

Na przykład, z i stałą ( tj. Jednolitą) gęstością możemy (przez inspekcję) wybrać jedną z wielu rozwiązania, takie jakdω=xkyldxdyp,

ω(x,y)=1l+1xkyl+1dx.

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 postaciuvt

tu+tw

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 iwvuuv.t0|vu|.xytdxdydt.t,kl.


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.

Rysunek przedstawiający osie wielokątów i głównych

#
# 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")
Whuber
źródło
+1 Wow, to świetna odpowiedź!
ameba
7

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)PA

C=[CXXCXYCXYCYY]

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:CXX=E[X2]XCYY=E[Y2]YCXY=E[XY]XY1AP

(1)CXX=1APx2dVCYY=1APy2dVCXY=1APxydV

Triangulacja

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

P=T1Tn

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

wprowadź opis zdjęcia tutaj

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

(2)CXX=1Ai=1nTix2dVCYY=1Ai=1nTiy2dVCXY=1Ai=1nTixydV

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(x1,y1),(x2,y2),(x3,y3)

vx=[x1x2x3]vy=[y1y2y3]1=[111]L=[100110111]

Następnie:

(3)Tx2dV=A6Tr(vxvxTL)Ty2dV=A6Tr(vyvyTL)TxydV=A12(1TvxvyT1+vxTvy)

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:vxivyiTi(3)(2)

(4)CXX=16i=1nTr(vxi(vxi)TL)CYY=16i=1nTr(vyi(vyi)TL)CXY=112i=1n(1Tvxi(vyi)T1+(vxi)Tvyi)

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 .CCPC(4)

user20160
źródło
2
+1 Można to uprościć, umożliwiając ukierunkowane trójkąty, eliminując w ten sposób potrzebę właściwej triangulacji. Zamiast tego możesz po prostu ustalić dowolne centrum i zsumować (podpisane) wartości nad trójkątami tak często się to robi, ponieważ jest znacznie mniej wybredne. Łatwo zauważyć, że takie podsumowanie jest zasadniczo tym samym, co zastosowanie twierdzenia Greena, ponieważ każdy termin w podsumowaniu jest ostatecznie funkcją krawędziPodejście to zostało zilustrowane w sekcji „Obszar” na quantdec.com/SYSEN597/GTKAV/section2/chapter_11.htm . OOPiPi+1:PiPi+1.
whuber
@whuber Ciekawe, dziękuję za zwrócenie na to uwagi
20160
Obie te odpowiedzi są dobre, choć trochę powyżej mojego wykształcenia. Gdy będę pewien, że je w pełni rozumiem, postaram się ustalić, kto dostanie nagrodę.
Ingolifs,