views:

31

answers:

1

I'm trying to fit a function consisting of several gauss bells to some experimental data. The method used is the nls function from R. But it is difficult to get the initial guess good enough, such that the method can converge.

Is it possible to visualize the initial guess BEFORE the optimization routine is called?

The code I'm working on is shown below (I cannot provide access to the datafile).

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  A: 

There you go: (see ?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)

If this doesn't solve your problem, think about rephrasing the question and adding some runnable code that illustrates the problem.

Joris Meys