Как получить область эллипса из двумерных нормальных распределенных данных?

12

У меня есть данные, которые выглядят так:

фигура

Я попытался применить нормальное распределение (оценка плотности ядра работает лучше, но мне не нужна такая большая точность), и это работает довольно хорошо. Плотность графика составляет эллипс.

Мне нужно получить эту функцию эллипса, чтобы решить, находится ли точка в области эллипса или нет. Как это сделать?

R или Mathematica код приветствуются.

matejuh
источник

Ответы:

18

Corsario предлагает хорошее решение в комментарии: используйте функцию плотности ядра для проверки на включение в набор уровней.

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

library(mvtnorm) # References rmvnorm()
set.seed(17)
p <- rmvnorm(1000, c(250000, 20000), matrix(c(100000^2, 22000^2, 22000^2, 6000^2),2,2))

Эллипсы определяются первым и вторым моментами данных:

center <- apply(p, 2, mean)
sigma <- cov(p)

Формула требует инверсии дисперсионно-ковариационной матрицы:

sigma.inv = solve(sigma, matrix(c(1,0,0,1),2,2))

Функция «высота» эллипса является отрицательной величиной логарифма двумерной нормальной плотности :

ellipse <- function(s,t) {u<-c(s,t)-center; u %*% sigma.inv %*% u / 2}

журнал(2πйе(Σ))

Чтобы проверить это , давайте нарисуем некоторые его контуры. Это требует создания сетки точек в направлениях x и y:

n <- 50
x <- (0:(n-1)) * (500000/(n-1))
y <- (0:(n-1)) * (50000/(n-1))

Вычислите функцию высоты в этой сетке и постройте ее:

z <- mapply(ellipse, as.vector(rep(x,n)), as.vector(outer(rep(0,n), y, `+`)))
plot(p, pch=20, xlim=c(0,500000), ylim=c(0,50000), xlab="Packets", ylab="Flows")
contour(x,y,matrix(z,n,n), levels=(0:10), col = terrain.colors(11), add=TRUE)

Контурный сюжет

Очевидно, это работает. Следовательно, тест для определения того, находится ли точка внутри эллиптического контура на уровне является(s,T)с

ellipse(s,t) <= c

Mathematica выполняет свою работу аналогичным образом: вычисляет матрицу дисперсии-ковариации данных, инвертирует ее, создает ellipseфункцию, и все готово.

Whuber
источник
Спасибо всем, особенно @whuber. Это именно то, что мне нужно.
matejuh
Btw. Есть ли простое решение для контуров оценки плотности ядра? Потому что, если я хочу быть более строгим, мои данные выглядят так: github.com/matejuh/doschecker_wiki_images/raw/master/… соответственно. github.com/matejuh/doschecker_wiki_images/raw/master/...
matejuh
Я не могу найти простое решение в R. Рассмотрим использование функции SmoothKernelDistribution в Mathematica 8.
whuber
2
Соответствует ли уровни уровню доверия? Я так не думаю. Как я могу это сделать, пожалуйста?
matejuh
Это требует нового вопроса, потому что вам нужно указать, на что вы надеетесь, и, судя по вашим графикам, есть сомнения относительно того, являются ли такие эллипсы адекватными описаниями данных.
whuber
9

Сюжет прост с ellipse()функцией mixtoolsпакета для R:

library(mixtools)
library(mvtnorm) 
set.seed(17)
p <- rmvnorm(1000, c(250000, 20000), matrix(c(100000^2, 22000^2, 22000^2, 6000^2),2,2))
plot(p, pch=20, xlim=c(0,500000), ylim=c(0,50000), xlab="Packets", ylab="Flows")
ellipse(mu=colMeans(p), sigma=cov(p), alpha = .05, npoints = 250, col="red") 

введите описание изображения здесь

Стефан Лоран
источник
5

Первый подход

Вы можете попробовать этот подход в Mathematica.

Давайте сгенерируем некоторые двумерные данные:

data = Table[RandomVariate[BinormalDistribution[{50, 50}, {5, 10}, .8]], {1000}];

Затем нам нужно загрузить этот пакет:

Needs["MultivariateStatistics`"]

И сейчас:

ellPar=EllipsoidQuantile[data, {0.9}]

дает вывод, который определяет 90% -ый доверительный эллипс. Значения, полученные из этого вывода, имеют следующий формат:

{Ellipsoid[{x1, x2}, {r1, r2}, {{d1, d2}, {d3, d4}}]}

x1 и x2 указывают точку, в которой эллипс в центре, r1 и r2 задают радиусы полуоси, а d1, d2, d3 и d4 задают направление выравнивания.

Вы также можете построить это:

Show[{ListPlot[data, PlotRange -> {{0, 100}, {0, 100}}, AspectRatio -> 1],  Graphics[EllipsoidQuantile[data, 0.9]]}]

Общая параметрическая форма эллипса:

ell[t_, xc_, yc_, a_, b_, angle_] := {xc + a Cos[t] Cos[angle] - b Sin[t] Sin[angle],
    yc + a Cos[t] Sin[angle] + b Sin[t] Cos[angle]}

И вы можете построить это так:

ParametricPlot[
    ell[t, ellPar[[1, 1, 1]], ellPar[[1, 1, 2]], ellPar[[1, 2, 1]], ellPar[[1, 2, 2]],
    ArcTan[ellPar[[1, 3, 1, 2]]/ellPar[[1, 3, 1, 1]]]], {t, 0, 2 \[Pi]},
    PlotRange -> {{0, 100}, {0, 100}}]

Вы можете выполнить проверку на основе чисто геометрической информации: если евклидово расстояние между центром эллипса (ellPar [[1,1]]) и вашей точкой данных больше, чем расстояние между центром эллипса и границей эллипс (очевидно, в том же направлении, в котором находится ваша точка), то эта точка данных находится вне эллипса.

Второй подход

Этот подход основан на плавном распределении ядра.

Вот некоторые данные, распространяемые аналогично вашим данным:

data1 = RandomVariate[BinormalDistribution[{.3, .7}, {.2, .3}, .8], 500];
data2 = RandomVariate[BinormalDistribution[{.6, .3}, {.4, .15}, .8], 500];
data = Partition[Flatten[Join[{data1, data2}]], 2];

Мы получаем плавное распределение ядра по этим значениям данных:

skd = SmoothKernelDistribution[data];

Мы получаем числовой результат для каждой точки данных:

eval = Table[{data[[i]], PDF[skd, data[[i]]]}, {i, Length[data]}];

Мы фиксируем порог и выбираем все данные, которые выше этого порога:

threshold = 1.2;
dataIn = Select[eval, #1[[2]] > threshold &][[All, 1]];

Здесь мы получаем данные, которые выходят за пределы региона:

dataOut = Complement[data, dataIn];

И теперь мы можем построить все данные:

Show[ContourPlot[Evaluate@PDF[skd, {x, y}], {x, 0, 1}, {y, 0, 1}, PlotRange -> {{0, 1}, {0, 1}}, PlotPoints -> 50],
ListPlot[dataIn, PlotStyle -> Darker[Green]],
ListPlot[dataOut, PlotStyle -> Red]]

Точки зеленого цвета - это точки выше порога, а точки красного цвета - точки ниже порога.

введите описание изображения здесь

VLC
источник
Спасибо, ваш второй подход очень помогает мне с распределением ядра. Я программист, а не статистик, и я новичок в Mathmatica и R, поэтому я очень ценю вашу помощь. В вашем втором подходе для меня ясно, как проверить одну точку, где она лежит. Но как это сделать при первом подходе? Я полагаю, что мне нужно сравнить мою точку зрения с определением эллипсоида. Можете ли вы предоставить, пожалуйста, как? Теперь я должен надеяться, что в R есть такие же определения, потому что мне нужно использовать его в RinRuby ...
matejuh
@matejuh Я только добавил еще несколько строк о первом подходе, который может направить вас к решению.
VLC
2

ellipseФункция в ellipseпакете для R будет генерировать эти эллипсы ( на самом деле многоугольник приближая эллипс). Вы могли бы использовать этот эллипс.

ellipseχ2

Грег Сноу
источник
1

Я нашел ответ по адресу: /programming/2397097/how-can-a-data-ellipse-be-superimpposed-on-a-ggplot2-scatterplot

#bootstrap
set.seed(101)
n <- 1000
x <- rnorm(n, mean=2)
y <- 1.5 + 0.4*x + rnorm(n)
df <- data.frame(x=x, y=y, group="A")
x <- rnorm(n, mean=2)
y <- 1.5*x + 0.4 + rnorm(n)
df <- rbind(df, data.frame(x=x, y=y, group="B"))

#calculating ellipses
library(ellipse)
df_ell <- data.frame()
for(g in levels(df$group)){
df_ell <- rbind(df_ell, cbind(as.data.frame(with(df[df$group==g,], ellipse(cor(x, y), 
                                         scale=c(sd(x),sd(y)), 
                                         centre=c(mean(x),mean(y))))),group=g))
}
#drawing
library(ggplot2)
p <- ggplot(data=df, aes(x=x, y=y,colour=group)) + geom_point(size=1.5, alpha=.6) +
  geom_path(data=df_ell, aes(x=x, y=y,colour=group), size=1, linetype=2)

введите описание изображения здесь

Парень л
источник