0

我正在尝试将由几个高斯钟组成的函数拟合到一些实验数据中。使用的方法是来自 R 的 nls 函数。但是很难得到足够好的初始猜测，使得该方法可以收敛。

是否可以在调用优化例程之前可视化初始猜测？

我正在处理的代码如下所示（我无法提供对数据文件的访问权限）。

library(signal)
# Load data from file
spectre <- read.table("LIA159.UXD")

# Extract variables and perform median filtering of the signal count
scatterangle <- spectre$V1
signal <- medfilt1(spectre$V2, n = 5)

#Perform a non linear fit of several gauss bells to the signal peaks
res <- nls( signal ~ bg + a*scatterangle 
    + h1*exp(-((scatterangle - m1)/s1)^2) 
    + h2*exp(-((scatterangle - m2)/s2)^2) 
    + h3*exp(-((scatterangle - m3)/s3)^2)
    + h4*exp(-((scatterangle - m4)/s4)^2)
    + h5*exp(-((scatterangle - m5)/s5)^2)
    + h6*exp(-((scatterangle - m6)/s6)^2)
    + h7*exp(-((scatterangle - m7)/s7)^2)
    , 
    start=list( 
        h1 =  2300, m1 = 23.42, s1 = 0.3, 
        h2 =  900,  m2 = 11.64, s2 = 0.2, 
        h3 =   100,     m3 = 34.80, s3 = 0.6, 
        h4 =   6,   m4 = 39.43, s4 = 1.3, 
        h5 =   3,   m5 = 46.83, s5 = 1.6, 
        h6 =  10,   m6 = 60.23, s6 = 0.3, 
        h7 =  10,   m7 = 61.46, s7 = 0.3, 
        bg=2, a = -0.1))

# Show the values of the fit
print(summary(res))

plot(signal ~ scatterangle, t='l', axes=F, xlab=expression(2*theta), 
ylab="")

# Draw the fitted function on top of the original data.
lines(scatterangle, predict(res, data.frame(scatterangle)), col='red')
4

1 回答 1

1

你去：（见？order）

set.seed(10)
bg <- rnorm(10000,2,0.1)

scatterangle <- runif(10000,5,35)
signal <- bg + -0.4*scatterangle +
     2000*exp(-((scatterangle - 24)/0.4)^2) +
     1000*exp(-((scatterangle - 12)/0.14)^2)+
     rnorm(10000,sd=100)

sv <- list(
    h1 =  2300, m1 = 23.42, s1 = 0.3,
    h2 =  900,  m2 = 11.64, s2 = 0.2,
    bg=2, a = -0.1)


res <- nls( signal ~ bg + a*scatterangle
    + h1*exp(-((scatterangle - m1)/s1)^2)
    + h2*exp(-((scatterangle - m2)/s2)^2)
    ,
    start=sv)

signal2 <- with(sv,{
    bg + a*scatterangle
    + h1*exp(-((scatterangle - m1)/s1)^2)
    + h2*exp(-((scatterangle - m2)/s2)^2)
    }
)

id <- order(scatterangle)
plot(signal[id]~scatterangle[id],
     t='l', axes=F, xlab=expression(2*theta),
    ylab="",col="grey")
lines(scatterangle[id],signal2[id],
    col='blue',lwd=2)
lines(scatterangle[id],
    predict(res, data.frame(scatterangle))[id],
    col='red',lwd=2)

如果这不能解决您的问题，请考虑重新表述问题并添加一些说明问题的可运行代码。

于 2010-09-25T22:22:20.580 回答