1

我为这个问题经历了许多解决方案，但找不到可靠的自动化解决方案。请在下面找到详细的自给自足的描述。

这些是数据：data.txt

x y
1 1
2 2
3 3
4 2
5 1

绘制为散点图：

t=read.table("data.txt",header=T)
plot(t$x,t$y,"l")

你会看到一个峰值，我现在的问题是：假设我对线性插值很满意，那么“曲线”的半最大值宽度是多少？因此，对于 x 的哪些值 x0，我有 f(x0)=max(y)/2，其中 f 是线性插值。我尝试了 approxfun 和一些内核密度，但我不想平滑我的数据。

非常欢迎任何输入。

4

2 回答 2

1

可能有很多更好的方法可以做到这一点，但这里有一个解决方案。如果我们一开始就知道你是如何进行插值的，那就更容易了。

# Extend your line for testing
y<-c(1,2,3,2,1,2,3,4,5,4,3)

# Split into linear segments.
segments<-c(1,diff(diff(y)))!=0
seg.points<-which(c(segments,TRUE))

# Calculate desired value
val<-max(y)/2

# Loop through and find value
res<-c()
for(i in 1:(length(seg.points)-1)) {
  start<-y[seg.points[i]]
  end<-y[seg.points[i+1]]
  if ((val>=start & val<=end) | (val<=start & val >=end)) {
    slope=(end-start)/(seg.points[i+1] - seg.points[i])
    res<-c(res,seg.points[i] + ((val - start) / slope))
  }
}
res
# [1] 2.5 3.5 6.5
于 2013-09-16T11:55:11.700 回答
1

这是一个简单的函数来做你想做的事，假设你的数据是单峰的，这里使用线性插值：

# FUNCTION TO INFER FWHM USING LINEAR INTERPOLATION
fwhm = function(x, y) { 
  halfheight = max(y)/2
  id.maxy = which.max(y)
  y1 = y[1:id.maxy]
  y2 = y[id.maxy:length(y)]
  x1 = x[1:id.maxy]
  x2 = x[id.maxy:length(y)]
  x.start = approx(x=y1,y=x1, xout=halfheight, method="linear")$y # or without interpolation: x[which.min(abs(y1-halfheight))]
  x.stop = approx(x=y2,y=x2, xout=halfheight, method="linear")$y # or without interpolation: x[which.min(abs(y2-halfheight))+id.maxy]
  fwhm = x.stop-x.start
  width = fwhm/(2*sqrt(2*log(2)))
  return(list(fwhm=fwhm, width=width))
  }

# GAUSSIAN PEAK FUNCTION WITH MODE AT u, WIDTH AT INFLECTION POINT w (FWHM=2.355*w) AND PEAK HEIGHT h
gauspeak=function(x, u, w, h=1) h*exp(((x-u)^2)/(-2*(w^2)))

# EXAMPLE
x = seq(0,200,length.out=1000)
y = gauspeak(x=x, u=100, w=10, h=100)
fwhm(x=x, y=y)
# $`fwhm`
# [1] 23.54934
# 
# $width
# [1] 10.00048

您也可以使用该spline()函数代替approx()三次样条插值而不是线性插值。如果您的数据是单峰的但有噪声，您可以首先使用平滑您的数据smooth.spline()（如果您的数据严格为正，则可能在一定log(y+1E-20)范围内，之后您可以反向转换为原始比例）。

于 2018-08-21T05:08:27.383 回答