用r找到双峰分布的局部最小值

Question

我的数据是预处理的图像数据,我想分开两个类.在理论上(并且希望在实践中),最佳阈值是双峰分布数据中两个峰值之间的局部最小值.

我的testdata是: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)？

Answer 1

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

首先,使用问题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")