Como encontrar a matriz de covariância de um polígono?

9

Imagine que você tem um polígono definido por um conjunto de coordenadas e seu centro de massa é . Você pode tratar o polígono como uma distribuição uniforme com um limite poligonal. (x1,y1)...(xn,yn)(0,0)insira a descrição da imagem aqui

Estou atrás de um método que encontrará a matriz de covariância de um polígono .

Suspeito que a matriz de covariância de um polígono esteja intimamente relacionada ao segundo momento da área , mas, se são equivalentes, não tenho certeza. As fórmulas encontradas no artigo da wikipedia que eu vinculei parecem (um palpite aqui, não está especialmente claro para mim neste artigo) referir-se à inércia rotacional em torno dos eixos x, ye z em vez dos eixos principais do polígono.

(Aliás, se alguém puder me indicar como calcular os eixos principais de um polígono, isso também seria útil)

É tentador apenas executar o PCA nas coordenadas , mas isso ocorre com o problema de que as coordenadas não estão necessariamente distribuídas uniformemente ao redor do polígono e, portanto, não são representativas da densidade do polígono. Um exemplo extremo é o contorno de Dakota do Norte, cujo polígono é definido por um grande número de pontos após o rio Vermelho, além de apenas mais dois pontos que definem a borda ocidental do estado.

Ingolifs
fonte
Por "achado", presumo que simplesmente a amostragem do polígono e o cálculo da covariância das amostras não é o que você tem em mente?
Stephan Kolassa
Além disso, você pode editar sua postagem para incluir coordenadas para o polígono, para que as pessoas possam brincar com ela?
101319 Stephan Stephanassass
11
@StephanKolassa Quero dizer, tratar o polígono como uma densidade de probabilidade bivariada uniforme com limite poligonal. Claro, você pode amostrar pontos e o limite seria a mesma coisa, mas estou procurando um método a-priori. A imagem é apenas uma ilustração da tinta que eu usei. Os dados do mundo real que pretendo usar são os contornos de estados e regiões.
Ingolifs 04/07/19
11
Você está certo de que o termo usual para "matriz de covariância" é momento de inércia ou segundo momento. Os principais eixos estão orientados em suas auto-orientações. A execução do PCA nas coordenadas está incorreta: equivale a supor que toda a massa está localizada nos vértices. Os métodos mais diretos de computação do baricentro - o primeiro momento - são discutidos no meu post em gis.stackexchange.com/a/22744/664 . Os segundos momentos são calculados da mesma maneira com pequenas modificações. Considerações especiais são necessárias na esfera.
whuber
2
μk,l(P)=Pxkyldxdy
ω d ω = x k y l d x d y . x k y l + 1 d x x k + 1 y l d yPωdω=xkyldxdy.xkyl+1dxxk+1yldy

Respostas:

10

Vamos fazer algumas análises primeiro.

Suponha que dentro do polígono sua densidade de probabilidade seja a função proporcional Então a constante de proporcionalidade é o inverso da integral de sobre o polígono,Pp(x,y).p

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

O baricentro do polígono é o ponto das coordenadas médias, computadas como seus primeiros momentos. O primeiro é

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

O tensor inercial pode ser representado como a matriz simétrica de segundos momentos calculada após a tradução do polígono para colocar seu baricentro na origem: ou seja, a matriz dos segundos momentos centrais

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

onde varia de a a O próprio tensor - também conhecido como matriz de covariância - é(k,l)(2,0)(1,1)(0,2).

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

Um PCA de produz os eixos principais de esses são os autovetores unitários dimensionados por seus autovalores.I(P)P :P:


Em seguida, vamos descobrir como fazer os cálculos. Como o polígono é apresentado como uma sequência de vértices que descreve seu limite orientado é natural invocarP,

Teorema de Green: onde é uma forma única definida em uma vizinhança de e

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

Por exemplo, com densidade constante ( ou seja , uniforme) podemos (por inspeção) selecionar um dos muitos soluções, comodω=xkyldxdyp,

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

O ponto disso é que a integral do contorno segue os segmentos de linha determinados pela sequência de vértices. Qualquer segmento de linha do vértice ao vértice pode ser parametrizado por uma variável real no formatouvt

tu+tw

onde é a direção normal da unidade de paraOs valores de portanto, variam de a Sob este parametrização e são funções lineares de e e são funções lineares do Assim, o integrando da integral de contorno sobre cada aresta se torna uma função polinomial de que é facilmente avaliada para pequeno ewvuuv.t0|vu|.xytdxdydt.t,kl.


A implementação dessa análise é tão direta quanto a codificação de seus componentes. No nível mais baixo, precisaremos de uma função para integrar uma forma polinomial sobre um segmento de linha. Funções de nível superior as agregam para calcular os momentos brutos e centrais para obter o baricentro e o tensor inercial, e finalmente podemos operar nesse tensor para encontrar os eixos principais (que são seus autovetores em escala). O Rcódigo abaixo executa este trabalho. Não faz pretensões de eficiência: destina-se apenas a ilustrar a aplicação prática da análise anterior. Cada função é direta e as convenções de nomenclatura são paralelas às da análise.

Incluído no código está um procedimento para gerar polígonos válidos, fechados, simplesmente interceptados e sem auto-interseção (deformando pontos aleatoriamente ao longo de um círculo e incluindo o vértice inicial como seu ponto final para criar um loop fechado). A seguir, algumas instruções para plotar o polígono, exibir seus vértices, unir o baricentro e plotar os eixos principais em vermelho (maior) e azul (menor), criando um sistema de coordenadas orientadas positivamente para o polígono.

Figura mostrando o polígono e os eixos principais

#
# 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
fonte
+1 Uau, esta é uma ótima resposta!
ameba
7

Edit: Não percebeu que whuber já tinha respondido. Vou deixar isso como um exemplo de outra abordagem (talvez menos elegante) para o problema.

A matriz de covariância

Deixe que ser um ponto aleatório a partir da distribuição uniforme sobre um polígono com área . A matriz de covariância é:(X,Y)PA

C=[CXXCXYCXYCYY]

onde é a variação de , é a variação de e é a covariância entre e . Isso assume média zero, uma vez que o centro de massa do polígono está localizado na origem. A distribuição uniforme atribui densidade de probabilidade constante a cada ponto em , portanto:CXX=E[X2]XCYY=E[Y2]YCXY=E[XY]XY1AP

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

Triangulação

Em vez de tentar integrar diretamente sobre uma região complicada como , podemos simplificar o problema particionando em sub-regiões triangulares:PPn

P=T1Tn

No seu exemplo, um possível particionamento é semelhante a este:

insira a descrição da imagem aqui

Existem várias maneiras de produzir uma triangulação (veja aqui ). Por exemplo, você pode calcular a triangulação de Delaunay dos vértices e descartar as arestas que ficam fora de (já que pode ser não-convexo, como no exemplo).P

As integrais sobre podem então ser divididas em somas de integrais sobre os triângulos:P

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

Um triângulo possui limites simples e agradáveis, portanto é mais fácil avaliar essas integrais.

Integrando sobre triângulos

Existem várias maneiras de integrar sobre triângulos. Nesse caso, usei um truque que envolve o mapeamento de um triângulo para o quadrado da unidade. Transformar em coordenadores bariêntricos pode ser uma opção melhor.

Aqui estão soluções para as integrais acima, para um triângulo arbitrário definido pelos vértices . Deixei:T(x1,y1),(x2,y2),(x3,y3)

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

Então:

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

Juntando tudo

Deixe e contenham as coordenadas x / y dos vértices para cada triângulo , como acima. Conecte em para cada triângulo, observando que os termos da área são cancelados. Isso fornece a solução:vxivyiTi(3)(2)

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

Eixos principais

Os eixos principais são dados pelos vetores próprios da matriz de covariância , assim como no PCA. Ao contrário do PCA, temos uma expressão analítica para , em vez de precisar calculá-la a partir de pontos de dados amostrados. Observe que os próprios vértices não são uma amostra representativa da distribuição uniforme em ; portanto, não se pode simplesmente pegar a matriz de covariância da amostra dos vértices. Mas, * é * uma função relativamente simples dos vértices, como visto em .CCPC(4)

user20160
fonte
2
+1 Isso pode ser simplificado permitindo triângulos orientados , eliminando assim a necessidade de uma triangulação adequada. Em vez disso, você pode simplesmente estabelecer um centro arbitrário e somar os valores (assinados) sobre os triângulos é assim que geralmente é feito porque é muito menos exigente. É fácil ver que esse somatório é essencialmente o mesmo que aplicar o Teorema de Green, porque cada termo no somatório é, em última análise, uma função da bordaEssa abordagem é ilustrada na seção "Área" em quantdec.com/SYSEN597/GTKAV/section2/chapter_11.htm . O P i P i + 1 : P i P i + 1 .OOPiPi+1:PiPi+1.
whuber
@whuber Interessante, obrigado por apontar isso
user20160
Ambas as respostas são boas, embora um pouco acima do meu nível de escolaridade. Quando tiver certeza de que os entendo completamente, tentarei descobrir quem recebe a recompensa.
Ingolifs 11/07/19