Найдите локальный минимум в бимодальном распределении с r

Мои данные представляют собой предварительно обработанные данные изображения, и я хочу разделить два класса. Теоретически (и, надеюсь, на практике) лучшим порогом является локальный минимум между двумя пиками в бимодальных распределенных данных.

Мои тестовые данные: http://www.file-upload.net/download-9365389/data.txt.html

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

datafile <- read.table("....txt")
data <- data$V1
hist(data)

d <- density(data) # returns the density data with defaults
hist(data,prob=TRUE)
lines(d) # plots the results

Но как продолжить?

Я бы вычислил первую и вторую производные функции плотности, чтобы найти локальные экстремумы, в частности, локальный минимум. Однако я понятия не имею, как это сделать в R, и density(test), похоже, не является нормальной функцией. Таким образом, пожалуйста, помогите мне: как я могу вычислить производные и найти локальный минимум ямы между двумя пиками в функции плотности density(test)?


r kernel statistics distribution kernel-density
person Iris    schedule 12.08.2014    source источник
comment
Можете ли вы добавить несколько примеров данных и продемонстрировать, что вы пробовали? Это должно облегчить вам помощь.   -  person Andrie    schedule 12.08.2014

Ответы (1)


arrow_upward
10
arrow_downward

Есть несколько способов сделать это.

Во-первых, используя d для плотности, как в вашем вопросе, d$x и d$y содержат значения x и y для графика плотности. Минимум возникает, когда производная dy/dx = 0. Поскольку значения x равномерно распределены, мы можем оценить dy, используя diff(d$y), и искать d$x, где abs(diff(d$y)) минимизируется:

d$x[which.min(abs(diff(d$y)))]
# [1] 2.415785

Проблема в том, что кривая плотности также может быть максимальной, когда dy/dx = 0. В этом случае минимум неглубокий, а максимумы пиковые, так что это работает, но вы не можете рассчитывать на это. .

Таким образом, второй способ использует optimize(...), который ищет локальный минимум в заданном интервале. optimize(...) нужна функция в качестве аргумента, поэтому мы используем approxfun(d$x,d$y) для создания функции интерполяции.

optimize(approxfun(d$x,d$y),interval=c(1,4))$minimum
# [1] 2.415791

Наконец, покажем, что это действительно минимум:

hist(data,prob=TRUE)
lines(d, col="red", lty=2)
v <- optimize(approxfun(d$x,d$y),interval=c(1,4))$minimum
abline(v=v, col="blue")

Другой подход, который на самом деле предпочтительнее, использует кластеризацию k-средних.

df <- read.csv(header=F,"data.txt")
colnames(df) = "X"

# bimodal
km <- kmeans(df,centers=2)
df$clust <- as.factor(km$cluster)
library(ggplot2)
ggplot(df, aes(x=X)) + 
  geom_histogram(aes(fill=clust,y=..count../sum(..count..)),
                     binwidth=0.5, color="grey50")+
  stat_density(geom="line", color="red")

Данные на самом деле выглядят скорее тримодальными, чем бимодальными.

# trimodal
km <- kmeans(df,centers=3)
df$clust <- as.factor(km$cluster)
library(ggplot2)
ggplot(df, aes(x=X)) + 
  geom_histogram(aes(fill=clust,y=..count../sum(..count..)),
                 binwidth=0.5, color="grey50")+
  stat_density(geom="line", color="red")

person jlhoward    schedule 13.08.2014
comment
Спасибо, @jlhoward, за подробный ответ! - person Iris; 15.08.2014