Question

Я хочу построить изоповерхность заданного c% -контура в трехмерной оценке плотности ядра. Затем я хочу узнать, какие точки находятся внутри этой 3d-формы.

Я покажу, что я подхожу к 2-й ситуации, чтобы проиллюстрировать свою проблему (код, имитируемый с R - Как найти точки в указанных c Contour и Как построить контурную линию, показывающую, где находятся 95% значений, в R и в ggplot2 ).

library(MASS)
library(misc3d)
library(rgl)
library(sp)

# Create dataset
set.seed(42)
Sigma <- matrix(c(15, 8, 5, 8, 15, .2, 5, .2, 15), 3, 3)
mv <- data.frame(mvrnorm(400, c(100, 100, 100),Sigma))

### 2d ###
# Create kernel density
dens2d <- kde2d(mv[, 1], mv[, 2], n = 40)
# Find the contour level defined in prob
dx <- diff(dens2d$x[1:2])
dy <- diff(dens2d$y[1:2])
sd <- sort(dens2d$z)
c1 <- cumsum(sd) * dx * dy 
prob <- .5
levels <- sapply(prob, function(x) { 
  approx(c1, sd, xout = 1 - x)$y
})

# Find which values are inside the defined polygon
ls <- contourLines(dens2d, level = levels)
pinp <- point.in.polygon(mv[, 1], mv[, 2], ls[[1]]$x, ls[[1]]$y)

# Plot it
plot(mv[, 1], mv[, 2], pch = 21, bg = "gray")
contour(dens2d, levels = levels, labels = prob,
        add = T, col = "red")
points(mv[pinp == 1, 1], mv[pinp == 1, 2], pch = 21, bg = "orange")

Итак, 50% контур определяется с помощью прибл., Контур создается с использованием contourLines, а затем point.in.polygon находит точки внутри этого контура.

Я хочу сделать то же самое, но в 3d ситуация. Вот что мне удалось:

### 3d ###
# Create kernel density
dens3d <- kde3d(mv[,1], mv[,2], mv[,3], n = 40)

# Find the contour level defined in prob
dx <- diff(dens3d$x[1:2])
dy <- diff(dens3d$y[1:2])
dz <- diff(dens3d$z[1:2])
sd3d <- sort(dens3d$d)
c3d <- cumsum(sd3d) * dx * dy * dz
levels <- sapply(prob, function(x) {
  approx(c3d, sd3d, xout = 1 - x)$y
})

# Find which values are inside the defined polygon
# # No idea

# Plot it
points3d(mv[,1], mv[,2], mv[,3], size = 2)
box3d(col = "gray")
contour3d(dens3d$d, level = levels, x = dens3d$x, y = dens3d$y, z = dens3d$z, #exp(-12)
          alpha = .3, color = "red", color2 = "gray", add = TRUE)
title3d(xlab = "x", ylab = "y", zlab = "z")

Итак, я не ушел далеко.

Я понимаю, что способ определения уровня в 3-м случае неверен, и я предполагаю, что проблема лежит в пределах c3d <- cumsum(sd3d) * dx * dy * dz, но я, честно говоря, не знаю, как действовать.

И, как только 3D-контур будет правильно определен , Я был бы признателен за любые советы о том, как приблизиться к тому, какие точки находятся в пределах этого контура.

Большое спасибо!

Редактировать : Основываясь на предложении от пользователя 2554330, я буду отредактируйте мой вопрос, чтобы добавить тестовый код, сравнивающий его или ее предложение с тем, который я разместил здесь. (Я понимаю, что этой цели использования контура в качестве вывода для новых точек данных не было в первоначальном вопросе, и я извиняюсь за эту поправку.)

Кроме того, я немного поспешил в своем комментарии ниже. Насколько хорошо эти два подхода работают в 2D-случае, зависит от размера выборки. При выборке n = 48 или около того подход user2554330 охватывает около 69% населения (тогда как подход, который я опубликовал, захватывает около 79%), но при выборке n = 400 или около того подход user2554330 захватывает около 79% (против 83%). ).

# Load libraries
library(MASS)
library(misc3d)
library(rgl)
library(sp)
library(oce)
library(akima)

# Create dataset
set.seed(42)
tn <- 1000 # number in pop
Sigma <- matrix(c(15, 8, 5, 8, 15, .2, 5, .2, 15), 3, 3)
mv <- data.frame(mvrnorm(tn, c(100, 100, 100),Sigma)) # population

prob <- .8 # rather than .5
simn <- 100 # number of simulations
pinp <- rep(NA, simn)
cuts <- pinp
sn <- 48 # sample size, at n = 400 user2554330 performs better

### 2d scenario
for (isim in 1:simn) {

  # Sample
  smv <- mv[sample(1:tn, sn), ]

  # Create kernel density
  dens2d <- kde2d(smv[, 1], smv[, 2], n = 40,
                  lims = c(min(smv[, 1]) - abs(max(smv[, 1]) - min(smv[, 1])) / 2,
                           max(smv[, 1]) + abs(max(smv[, 1]) - min(smv[, 1])) / 2,
                           min(smv[, 2]) - abs(max(smv[, 2]) - min(smv[, 2])) / 2,
                           max(smv[, 2]) + abs(max(smv[, 2]) - min(smv[, 2])) / 2))



  # Approach based on https://stackoverflow.com/questions/30517160/r-how-to-find-points-within-specific-contour
  # Find the contour level defined in prob
  dx <- diff(dens2d$x[1:2])
  dy <- diff(dens2d$y[1:2])
  sd <- sort(dens2d$z)
  c1 <- cumsum(sd) * dx * dy 
  levels <- sapply(prob, function(x) { 
    approx(c1, sd, xout = 1 - x)$y
  })
  # Find which values are inside the defined polygon
  ls <- contourLines(dens2d, level = levels)
  # Note below that I check points from "population"
  pinp[isim] <- sum(point.in.polygon(mv[, 1], mv[, 2], ls[[1]]$x, ls[[1]]$y)) / tn



  # Approach based on user2554330
  # Find the estimated density at each observed point
  sdatadensity<- bilinear(dens2d$x, dens2d$y, dens2d$z, 
                    smv[,1], smv[,2])$z
  # Find the contours
  levels2 <- quantile(sdatadensity, probs = 1- prob, na.rm = TRUE)
  # Find within
  # Note below that I check points from "population"
  datadensity <- bilinear(dens2d$x, dens2d$y, dens2d$z, 
                    mv[,1], mv[,2])$z
  cuts[isim] <- sum(as.numeric(cut(datadensity, c(0, levels2, Inf))) == 2, na.rm = T) / tn

}

summary(pinp)
summary(cuts)

> summary(pinp)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0030  0.7800  0.8205  0.7950  0.8565  0.9140 
> summary(cuts)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.5350  0.6560  0.6940  0.6914  0.7365  0.8120

Я также пытался увидеть, как подход user2554330 работает в 3D-ситуации, с помощью кода ниже:

# 3d scenario
for (isim in 1:simn) {

  # Sample
  smv <- mv[sample(1:tn, sn), ]

  # Create kernel density
  dens3d <- kde3d(smv[,1], smv[,2], smv[,3], n = 40,
                  lims = c(min(smv[, 1]) - abs(max(smv[, 1]) - min(smv[, 1])) / 2,
                           max(smv[, 1]) + abs(max(smv[, 1]) - min(smv[, 1])) / 2,
                           min(smv[, 2]) - abs(max(smv[, 2]) - min(smv[, 2])) / 2,
                           max(smv[, 2]) + abs(max(smv[, 2]) - min(smv[, 2])) / 2,
                           min(smv[, 3]) - abs(max(smv[, 3]) - min(smv[, 3])) / 2,
                           max(smv[, 3]) + abs(max(smv[, 3]) - min(smv[, 3])) / 2))

  # Approach based on user2554330
  # Find the estimated density at each observed point
  sdatadensity <- approx3d(dens3d$x, dens3d$y, dens3d$z, dens3d$d, 
                          smv[,1], smv[,2], smv[,3])
  # Find the contours
  levels <- quantile(sdatadensity, probs = 1 - prob, na.rm = TRUE)
  # Find within
  # Note below that I check points from "population"
  datadensity <- approx3d(dens3d$x, dens3d$y, dens3d$z, dens3d$d, 
                          mv[,1], mv[,2], mv[,3])
  cuts[isim] <- sum(as.numeric(cut(datadensity, c(0, levels, Inf))) == 2, na.rm = T) / tn

}

summary(cuts)

> summary(cuts)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.1220  0.1935  0.2285  0.2304  0.2620  0.3410

Я бы предпочел определить контур так, чтобы указанная вероятность (близка к) вероятность получения будущих точек данных, взятых из той же совокупности, даже когда выборка n относительно мала (т. е. <50). </p>

user2554330 · Answer 1 · 01 февраля 2020

Вместо того, чтобы пытаться определить, какие точки находятся внутри контура, я бы попытался оценить плотность в каждой точке и раскрасить точки в соответствии с тем, как это значение сравнивается с уровнем контура. Это может привести к другому решению для нескольких точек вблизи границы, но должно быть довольно близко.

Чтобы выполнить эту оценку, вы можете использовать функцию oce::approx3d для оценки плотности.

Другая вещь, которую я бы сделал, - это выбрать контур на основе квантилей наблюдаемых плотностей, а не пытаться моделировать трехмерный интеграл от оцененной плотности.

Вот код для всего этого:

library(MASS)
library(misc3d)
library(rgl)
library(oce)
#> Loading required package: testthat
#> Loading required package: gsw

# Create dataset
set.seed(42)
Sigma <- matrix(c(15, 8, 5, 8, 15, .2, 5, .2, 15), 3, 3)
mv <- data.frame(mvrnorm(400, c(100, 100, 100),Sigma))

### 3d ###
# Create kernel density
dens3d <- kde3d(mv[,1], mv[,2], mv[,3], n = 40)

# Find the estimated density at each observed point
datadensity <- approx3d(dens3d$x, dens3d$y, dens3d$z, dens3d$d, 
                        mv[,1], mv[,2], mv[,3])

# Find the contours
prob <- .5
levels <- quantile(datadensity, probs = prob, na.rm = TRUE)

# Plot it
colours <- c("gray", "orange")
cuts <- cut(datadensity, c(0, levels, Inf))
for (i in seq_along(levels(cuts))) {
  gp <- as.numeric(cuts) == i
  spheres3d(mv[gp,1], mv[gp,2], mv[gp,3], col = colours[i], radius = 0.2)
}
box3d(col = "gray")
contour3d(dens3d$d, level = levels, x = dens3d$x, y = dens3d$y, z = dens3d$z, #exp(-12)
          alpha = .1, color = "red", color2 = "gray", add = TRUE)
title3d(xlab = "x", ylab = "y", zlab = "z")

А вот график, который был создан:

Создайте% -контур в 3D-плотности ядра и найдите точки в этом контуре

Пожалуйста, войдите или зарегистрируйтесь чтобы ответить на этот вопрос.

1 Ответ

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

Создайте% -контур в 3D-плотности ядра и найдите точки в этом контуре

Пожалуйста, войдите или зарегистрируйтесь чтобы ответить на этот вопрос.

1 Ответ

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

Похожие темы