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

そこに行きます:(？注文を参照）

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 に答える