Bagaimana menemukan matriks kovarians poligon?

Bayangkan Anda memiliki poligon yang didefinisikan oleh seperangkat koordinat $(x_1,y_1)...(x_n,y_n)$ dan pusat massanya adalah pada $(0,0)$ . Anda dapat memperlakukan poligon sebagai distribusi seragam dengan batas poligon.

Saya mencari metode yang akan menemukan matriks kovarian poligon .

Saya menduga bahwa matriks kovarians poligon terkait erat dengan momen kedua dari area , tetapi apakah mereka setara, saya tidak yakin. Rumus yang ditemukan dalam artikel wikipedia yang saya tautkan tampaknya (tebakan di sini, tidak terlalu jelas bagi saya dari artikel) untuk merujuk pada inersia rotasi di sekitar sumbu x, y dan z daripada sumbu utama poligon.

(Kebetulan, kalau ada yang bisa mengarahkan saya ke cara menghitung sumbu utama poligon, itu juga akan berguna bagi saya)

Ini tergoda untuk hanya melakukan PCA pada koordinat , tetapi hal itu berjalan ke masalah bahwa koordinat tidak selalu tersebar merata di sekitar poligon, dan karena itu tidak mewakili kepadatan poligon. Contoh ekstrem adalah garis besar North Dakota, yang poligonnya didefinisikan oleh sejumlah besar titik mengikuti sungai Merah, ditambah hanya dua titik lagi yang menentukan tepi barat negara bagian.

Mari kita lakukan analisis terlebih dahulu.

Misalkan dalam poligon kepadatan probabilitasnya adalah fungsi proporsional Maka konstanta proporsionalitas adalah kebalikan dari integral atas poligon,$\mathcal{P}$$p(x,y).$$p$

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

The barycenter dari poligon adalah titik koordinat rata-rata, dihitung sebagai momen pertama mereka. Yang pertama adalah

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

The tensor inersia dapat direpresentasikan sebagai array simetris saat kedua dihitung setelah menerjemahkan poligon untuk menempatkan barycenter pada titik asal: yaitu, matriks momen kedua pusat

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

di mana berkisar dari hingga hingga Tensor itu sendiri - alias matriks kovarians - adalah$(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 menghasilkan sumbu utama dari ini adalah vektor satuan eigen yang diukur oleh nilai eigennya.$I(\mathcal{P})$P :$\mathcal{P}:$

---

Selanjutnya, mari kita cari tahu cara melakukan perhitungan. Karena poligon disajikan sebagai urutan simpul yang menggambarkan batas berorientasi adalah wajar untuk memanggil$\partial\mathcal P,$

Teorema Green: mana adalah satu-bentuk yang didefinisikan di lingkungan dan

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

Misalnya, dengan dan kepadatan konstan ( mis. , Seragam) kita dapat (dengan inspeksi) memilih salah satu dari banyak solusi, seperti$\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.$$

Intinya adalah bahwa integral kontur mengikuti segmen garis yang ditentukan oleh urutan simpul. Segmen baris apa pun dari vertex ke vertex dapat diparameterisasi dengan variabel nyata dalam formulir$\mathbf{u}$$\mathbf{v}$$t$

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

di mana adalah arah normal satuan dari keNilai karena itu berkisar dari hingga Di bawah parameterisasi ini dan adalah fungsi linear dari dan dan adalah fungsi linear dari Dengan demikian integran dari integral kontur atas setiap tepi menjadi fungsi polinom dari yang mudah dievaluasi untuk kecil dan$\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.$

---

Menerapkan analisis ini sama mudahnya dengan mengkodekan komponen-komponennya. Pada level terendah kita akan memerlukan fungsi untuk mengintegrasikan satu-bentuk polinomial pada segmen garis. Fungsi tingkat yang lebih tinggi akan mengumpulkan ini untuk menghitung momen mentah dan pusat untuk mendapatkan barycenter dan tensor inersia, dan akhirnya kita dapat beroperasi pada tensor tersebut untuk menemukan sumbu utama (yang merupakan vektor eigen yang diskalakan). The Rkode di bawah melakukan pekerjaan ini. Itu tidak membuat pretensi efisiensi: ini dimaksudkan hanya untuk menggambarkan penerapan praktis dari analisis sebelumnya. Setiap fungsi mudah dan konvensi penamaan sejajar dengan analisis.

Termasuk dalam kode adalah prosedur untuk menghasilkan poligon tertutup, hanya terhubung, non-self-berpotongan yang valid (dengan mendeformasi secara acak titik-titik di sepanjang lingkaran dan termasuk titik awal sebagai titik terakhir untuk membuat loop tertutup). Berikut ini adalah beberapa pernyataan untuk memplot poligon, menampilkan simpulnya, berdampingan dengan barycenter, dan memplot sumbu utama dalam warna merah (terbesar) dan biru (terkecil), menciptakan sistem koordinat poligon-sentris berorientasi positif.

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

Sunting: Tidak menyadari bahwa whuber sudah menjawab. Saya akan meninggalkan ini sebagai contoh dari pendekatan lain (mungkin kurang elegan) untuk masalah ini.

Matriks kovarians

Mari menjadi titik acak dari distribusi seragam pada poligon dengan luas . Matriks kovarians adalah:$(X,Y)$$P$$A$

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

di mana adalah varian , adalah varian , dan adalah kovarians antara dan . Ini mengasumsikan nol rata-rata, karena pusat massa poligon terletak di titik asal. Distribusi seragam memberikan kepadatan probabilitas konstan ke setiap titik dalam , jadi:$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}$$

Triangulasi

Alih-alih mencoba mengintegrasikan langsung ke wilayah rumit seperti , kita dapat menyederhanakan masalah dengan mempartisi ke dalam subregional segitiga:$P$$P$$n$

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

Dalam contoh Anda, satu kemungkinan partisi tampak seperti ini:

Ada berbagai cara untuk menghasilkan triangulasi (lihat di sini ). Misalnya, Anda bisa menghitung triangulasi simpul Delaunay , kemudian membuang sisi-sisi yang berada di luar (karena ini mungkin bukan cembung seperti pada contoh).$P$

Integral over kemudian dapat dibagi menjadi jumlah integral atas segitiga:$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}$$

Segitiga memiliki batas yang bagus dan sederhana sehingga integral ini lebih mudah untuk dievaluasi.

Mengintegrasikan lebih dari segitiga

Ada berbagai cara untuk berintegrasi dengan segitiga. Dalam hal ini, saya menggunakan trik yang melibatkan pemetaan segitiga ke unit square. Mengubah ke koordinat barycentric mungkin menjadi pilihan yang lebih baik.

Berikut adalah solusi untuk integral di atas, untuk segitiga sembarang didefinisikan oleh simpul . Membiarkan:$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]$$

Kemudian:

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

Menyatukan semuanya

Biarkan dan berisi koordinat x / y dari simpul untuk setiap segitiga , seperti di atas. Masukkan ke dalam untuk setiap segitiga, mencatat bahwa ketentuan daerah dibatalkan. Ini memberikan solusinya:$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}$$

Kapak utama

Sumbu utama diberikan oleh vektor eigen dari matriks kovarian , seperti pada PCA. Tidak seperti PCA, kami memiliki ekspresi analitik untuk , daripada harus memperkirakannya dari titik data sampel. Perhatikan bahwa simpul itu sendiri bukan sampel yang representatif dari distribusi seragam pada , jadi orang tidak bisa begitu saja mengambil matriks kovarians sampel dari simpul tersebut. Tapi, * adalah * fungsi simpul yang relatif sederhana, seperti yang terlihat pada .$C$$C$$P$$C$$(4)$