用 r 在双峰分布中找到局部最小值
Posted
技术标签:
【中文标题】用 r 在双峰分布中找到局部最小值【英文标题】:Find local minimum in bimodal distribution with r 【发布时间】:2014-10-05 12:45:59 【问题描述】:我的数据是经过预处理的图像数据,我想将两个类分开。理论上(并希望在实践中)最佳阈值是双峰分布数据中两个峰值之间的局部最小值。
我的测试数据是:http://www.file-upload.net/download-9365389/data.txt.html
我尝试关注this thread: 我绘制了直方图并计算了核密度函数:
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)
中两个峰值之间的坑的局部最小值?
【问题讨论】:
您能添加一些示例数据并演示您尝试过的内容吗?这应该可以更轻松地为您提供帮助。 【参考方案1】:有几种方法可以做到这一点。
首先,在您的问题中使用d
作为密度,d$x
和d$y
包含密度图的 x 和 y 值。当导数 dy/dx = 0 时出现最小值。由于 x 值是等距的,我们可以使用 diff(d$y)
估计 dy,并在 abs(diff(d$y))
最小化的位置寻找 d$x
:
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")
【讨论】:
感谢@jlhoward 的详细解答!以上是关于用 r 在双峰分布中找到局部最小值的主要内容,如果未能解决你的问题,请参考以下文章