2010-09-25 17 views
0

我试图将几个高斯响铃组成的函数拟合成一些实验数据。所使用的方法是R的nls函数。但很难获得足够好的初始猜测,使得该方法可以收敛。初步猜测用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') 

回答

1

你去那里(见订货?)

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) 

如果这没有解决您的问题,想改写这个问题,并补充说,说明了这个问题的一些运行的代码。