5

我的数据是经过预处理的图像数据，我想分开两个类。理论上（并希望在实践中）最佳阈值是双峰分布数据中两个峰值之间的局部最小值。

我的测试数据是：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)？

4

1 回答 1

10

有几种方法可以做到这一点。

首先，使用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-means 聚类。

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

于 2014-08-13T01:25:08.110 回答