¿Cómo encontrar la matriz de covarianza de un polígono?

Question

Imagina que tienes un polígono definido por un conjunto de coordenadas $(x_1,y_1)...(x_n,y_n)$ y su centro de masa está en $(0,0)$ . Puede tratar el polígono como un distribución uniforme con un límite poligonal.

Busco un método que encuentre el matriz de covarianza de un polígono .

Sospecho que la matriz de covarianza de un polígono está estrechamente relacionada con la _segundo momento del área_ pero no estoy seguro de que sean equivalentes. Las fórmulas que se encuentran en el artículo de la wikipedia que he enlazado parecen (es una suposición, no me queda especialmente claro en el artículo) referirse a la inercia rotacional alrededor de los ejes x, y y z en lugar de los ejes principales del polígono.

(Por cierto, si alguien puede indicarme cómo calcular los ejes principales de un polígono, también me sería útil)

Resulta tentador limitarse a realizar el PCA en las coordenadas pero al hacerlo se encuentra con el problema de que las coordenadas no están necesariamente repartidas uniformemente por el polígono, y por tanto no son representativas de la densidad del mismo. Un ejemplo extremo es el contorno de Dakota del Norte, cuyo polígono está definido por un gran número de puntos que siguen el río Rojo, además de sólo dos puntos más que definen el borde occidental del estado.

Answer 1

Hagamos primero un análisis.

Supongamos que dentro del polígono $\mathcal{P}$ su densidad de probabilidad es una función proporcional $p(x,y).$ Entonces la constante de proporcionalidad es la inversa de la integral de $p$ sobre el polígono,

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

El baricentro del polígono es el punto de coordenadas medias, calculado como sus primeros momentos. El primero es

$$\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.$$

El tensor de inercia puede representarse como la matriz simétrica de segundos momentos calculada después de trasladar el polígono para poner su baricentro en el origen: es decir, la matriz de segundos momentos centrales

$$\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$$

donde $(k,l)$ van desde $(2,0)$ a $(1,1)$ a $(0,2).$ El propio tensor también conocido como matriz de covarianza--es

$$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})}.$$

Un PCA de $I(\mathcal{P})$ produce el ejes principales de $\mathcal{P}:$ son los vectores propios unitarios escalados por sus valores propios.

---

A continuación, vamos a ver cómo hacer los cálculos. Como el polígono se presenta como una secuencia de vértices que describen su límite orientado $\partial\mathcal P,$ es natural invocar

Teorema de Green: $$\iint_{\mathcal{P}} \mathrm{d}\omega = \oint_{\partial\mathcal{P}}\omega$$ donde $\omega = M(x,y)\mathrm{d}x + N(x,y)\mathrm{d}y$ es una forma única definida en una vecindad de $\mathcal{P}$ y $$\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.$$

Por ejemplo, con $\mathrm{d}\omega = x^k y^l \mathrm{d}x\mathrm{d}y$ y constante ( es decir , uniforme) densidad $p,$ podemos (por inspección) seleccionar una de las muchas soluciones, como $$\omega(x,y) = \frac{-1}{l+1}x^k y^{l+1}\mathrm{d}x.$$

El sentido de esto es que la integral del contorno sigue los segmentos de línea determinados por la secuencia de vértices. Cualquier segmento de línea desde el vértice $\mathbf{u}$ al vértice $\mathbf{v}$ puede parametrizarse mediante una variable real $t$ en la forma

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

donde $\mathbf{w} \propto \mathbf{v}-\mathbf{u}$ es la dirección normal unitaria de $\mathbf{u}$ a $\mathbf{v}.$ Los valores de $t$ por lo tanto, van desde $0$ a $|\mathbf{v}-\mathbf{u}|.$ Con esta parametrización $x$ y $y$ son funciones lineales de $t$ y $\mathrm{d}x$ y $\mathrm{d}y$ son funciones lineales de $\mathrm{d}t.$ Así, el integrando de la integral de contorno sobre cada arista se convierte en un función polinómica de $t,$ que es fácilmente evaluable para pequeños $k$ y $l.$

---

Aplicación de este análisis es tan sencillo como codificar sus componentes. En el nivel más bajo necesitaremos una función que integre un polinomio de una forma sobre un segmento de línea. Las funciones de nivel superior las agregarán para calcular los momentos brutos y centrales para obtener el baricentro y el tensor de inercia, y finalmente podemos operar sobre ese tensor para encontrar los ejes principales (que son sus vectores propios escalados). La página web R El código siguiente realiza este trabajo. No tiene ninguna pretensión de eficacia: sólo pretende ilustrar la aplicación práctica del análisis anterior. Cada función es sencilla y las convenciones de nomenclatura son paralelas a las del análisis.

El código incluye un procedimiento para generar polígonos válidos cerrados, simplemente conectados y no autointerseccionados (deformando aleatoriamente puntos a lo largo de un círculo e incluyendo el vértice inicial como punto final para crear un bucle cerrado). A continuación hay unas cuantas sentencias para trazar el polígono, mostrar sus vértices, adosar el baricentro y trazar los ejes principales en rojo (el más grande) y azul (el más pequeño), creando un sistema de coordenadas orientado positivamente y centrado en el polígono.

#
# 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")

Answer 2

Edit: No me di cuenta de que whuber ya había contestado. Voy a dejar esto como un ejemplo de otro (tal vez menos elegante) enfoque del problema.

La matriz de covarianza

Dejemos que $(X,Y)$ sea un punto aleatorio de la distribución uniforme en un polígono $P$ con el área $A$ . La matriz de covarianza es:

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

donde $C_{XX} = E[X^2]$ es la varianza de $X$ , $C_{YY} = E[Y^2]$ es la varianza de $Y$ et $C_{XY} = E[XY]$ es la covarianza entre $X$ y $Y$ . Esto supone una media cero, ya que el centro de masa del polígono está situado en el origen. La distribución uniforme asigna una densidad de probabilidad constante $\frac{1}{A}$ a cada punto de $P$ Así que..:

$$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}$$

Triangulación

En lugar de intentar integrar directamente sobre una región complicada como $P$ podemos simplificar el problema dividiendo $P$ en $n$ subregiones triangulares:

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

En su ejemplo, una posible partición es la siguiente:

Hay varias formas de producir una triangulación (véase aquí ). Por ejemplo, se puede calcular el Triangulación de Delaunay de los vértices, entonces descarta las aristas que caen fuera $P$ (ya que puede ser no convexo como en el ejemplo).

Integrales sobre $P$ puede dividirse en sumas de integrales sobre los triángulos:

$$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}$$

Un triángulo tiene unos límites agradables y sencillos, por lo que estas integrales son más fáciles de evaluar.

Integrar sobre triángulos

Hay varias formas de integrar sobre triángulos. En este caso, he utilizado un truco que implica mapear un triángulo en el cuadrado de la unidad. La transformación a coordenadas baricéntricas podría ser una mejor opción.

Aquí están las soluciones para las integrales anteriores, para un triángulo arbitrario $T$ definido por los vértices $(x_1,y_1), (x_2,y_2), (x_3,y_3)$ . Déjalo:

$$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]$$

Entonces:

$$\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}$$

Poner todo en orden

Dejemos que $v_x^i$ y $v_y^i$ contienen las coordenadas x/y de los vértices de cada triángulo $T_i$ como en el caso anterior. Enchufe $(3)$ en $(2)$ para cada triángulo, observando que los términos de área se cancelan. Esto da la solución:

$$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}$$

Ejes principales

Los ejes principales vienen dados por los vectores propios de la matriz de covarianza $C$ como en el PCA. A diferencia del ACP, tenemos una expresión analítica para $C$ en lugar de tener que estimarlo a partir de puntos de datos muestreados. Nótese que los vértices en sí no son una muestra representativa de la distribución uniforme en $P$ por lo que no se puede tomar simplemente la matriz de covarianza muestral de los vértices. Pero, $C$ *es* una función relativamente simple de los vértices, como se ve en $(4)$ .