Initialisation de la visualisation avec la fonction nls

Question

J'essaie d'ajuster une fonction composée de plusieurs cloches Gauss à des données expérimentales. La méthode utilisée est la fonction nls de R. Mais il est difficile d'obtenir une estimation initiale suffisamment bonne pour que la méthode puisse converger.

Est-il possible de visualiser la supposition initiale AVANT l'appel de la routine d'optimisation?

Le code sur lequel je travaille est illustré ci-dessous (je ne peux pas donner accès au fichier de données).

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') 

Answer 1

Là vous allez: (voir commande?)

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) 

Si cela ne résout pas votre problème, pensez à reformuler la question et en ajoutant un code qui illustre runnable le problème.