nls 当数据是卷积光谱时（变量数量可变）

Posted 2023-02-25

【中文标题】nls 当数据是卷积光谱时（变量数量可变）

【英文标题】：nls when the data is a convoluted spectra (variable number of variables)

【发布时间】：2014-01-27 12:34:47

【问题描述】：

我有一些光谱数据在绘制时看起来像多重峰之一： http://journals.prous.com/journals/dof/19982303/html/df230301/images/keiferf3.gif

如何在图像中看到，所有峰彼此之间非常接近，所以我想使用 nls 函数进行一些反卷积，就像之前发布的 (R: Fitting Gaussian peaks to density plot data using nls) 一样，但使用洛伦兹函数代替：

y <- 1/(pi*a*(1+((x-x0)/a)^2))

在我的例子中，x0 是峰值最大值（而 length(x0) 是峰值的数量），所以我只需要优化 'a'。

但是，我的问题与执行此操作无关，而是编写一个强大的脚本来对任何光谱进行反卷积，将峰数作为输入信息。

我的第一个想法是编写洛伦兹函数并将“a”作为向量（之后应用所有洛伦兹曲线的总和），但 R 无法识别这种结构：

for (i in 1:length(x0)) 
    f[i]<-function(a)  y <- 1/(pi*a[i]*(1+((x-x0[i])/a[i])^2)) 


fit <- nls(sum(f[1:length(x0)]), start=list(a=rep(1, times=length(x0))))

更新：

这是我的示例，采用 .csv 格式 (https://drive.google.com/file/d/0B66EHLI5AufhbjlWcW9rYXl1UFk/edit?usp=sharing)。数据填充为 2 行。第一个具有频率（以 ppm 为单位），第二个具有强度。对于这个数据，我会选择 5 个峰值，所以我会在 formula=f[1]+f[2]+f[3]+f[4]+f[5] 上做 'nls' 并且我会有 5 个参数(a[1:5]) 进行评估。

【问题讨论】：

旁注：f <- 1/pi*a*(1+(x-x0)/a^2) 是循环的矢量化替换。

我知道我的代码写得不好，但我认为发布代码有助于理解我的要求。

当然，不用担心 :) 这只是 R 特定的提示。

查看引用的问题：显然，这仅在峰分离良好时才有效（例如瑞利标准）。鉴于此，我认为您需要向nls“解释”您有一组适合的变量。也许让list 的元素是各种f[i] ？

您能否提供样本数据，比如链接中的光谱数据？

【参考方案1】：

我找到了'as.formula' 命令，它允许将任何字符串转换为公式。因此，我设法创建了一个创建洛伦兹曲线总和的 for 循环。在我的示例中，参数 a[1:5] 现在定义为 a-e，但使用 sprintf 我也可以使用向量命名法。

cac<-"abcdefghijklmnopqrstuvwxyz"
for (i in 1:length(x0)) 
    letter<-substr(cac, i,i)
    formulae[i]<-sprintf("1/(pi*%s*(1+((spectra[1,]-%f)/%s)^2))", letter,x0[i],letter)
    coeff[i]<-sprintf("%s=1", letter)

formula2<-paste(formulae, collapse="+")
formulo<-paste("spectra[2,] ~", formula2, sep="")
coeffs<-paste(coeff, collapse=",")

fit<-as.formula(paste("nls(",formulo,",start=list(",coeffs,"))", sep=""))

因此，所有这些代码都为我提供了发布示例的以下公式：

"nls(spectra[2,] ~1/(pi*a*(1+((spectra[1,]-2.156460)/a)^2))+1/(pi*b*(1+((spectra[1,]-2.170150)/b)^2))+1/(pi*c*(1+((spectra[1,]-2.184820)/c)^2))+1/(pi*d*(1+((spectra[1,]-2.163550)/d)^2))+1/(pi*e*(1+((spectra[1,]-2.142040)/e)^2)), start=list(a=1,b=1,c=1,d=1,e=1))

但是，这个公式似乎不起作用，但这不是这个帖子的目的，所以我可以说现在它已经关闭了（我开了一个新的试图解决它）。