2

我有一些光谱数据，在绘制时看起来像多重峰之一：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 是峰值最大值（长度（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]) 进行评估。

4

1 回答 1

0

我找到了'as.formula' 命令，它允许将任何字符串转换为公式。因此，我设法创建了一个创建洛伦兹曲线总和的 for 循环。在我的示例中，参数 a[1:5] 现在定义为 ae，但使用 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))

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

于 2014-01-29T16:36:15.430 回答