Как определить плавающие параметры при использовании NLS модель R? - PullRequest
0 голосов
/ 26 марта 2020

Я работаю над сценарием R, который будет соответствовать некоторым данным, определенным функцией Лоренца. В общем, это довольно легко сделать, однако я хочу иметь возможность фиксировать указанные c параметры, чтобы они не включались в подбор. Как правило, я бы не включал параметры в команду 'start' в функцию nls, однако я хочу что-то более надежное, которое будет работать с неограниченным числом пиков, т.е. я хочу иметь возможность хранить определенные c параметры в векторе, Например, модель фиксирует [2], но соответствует [1], a [3] et c. Приведенный ниже код должен также помочь объяснить мою проблему:

rm(list = ls()) #clear history

stepsize<-25/500 #Define step size
x=seq(-10, 10, by=stepsize) # Define x-axis

#Generate dummy data
a0=0.5 #width
b0=2 #x-position of peak
c0=1

# Create y-axis based on Lorentzian function with random noise
y <- ((1/pi)*(a0*0.5/((x-(b0-c0/2))^2+(a0/2)^2))+(1/pi)*(a0*0.5/((x-(b0+c0/2))^2+(a0/2)^2)))+rnorm(401,0,0.1) 

plot(x,y) # Plot data to check
f=0
fit<-function(x,a,b){ #Define function which will be used to "Fit" the dummy data.
  for(i in 1:2){
    f<-f+(1/pi)*(a[i]*0.5/((x-b[i])^2+(a[i]/2)^2))  
  }
  f
}

#Define some approximations for the fit
a00<-c(0.2,0.4)
b00<-c(1.3,2.55)

#add the rough guess to the plot
lines(x,fit(x,a00,b00))


# ---- Run the fitting model ---- #
# To hold parameters, give an initial guess and then remove the parameter from the "start" list. 
# The problem is here. I want to be able to selectively remove certain parameters from the fit, 
# e.g. fix a[2] but allow a[1] to float.

# ##option 1: works but does not include any fixed parameters 
# m<-nls(y ~ fit(x,a,b), start = list(a=a00, b=b00),
#        control=nls.control(maxiter=50, tol=1e-5, minFactor = 1/2048, warnOnly=TRUE),
#        trace=TRUE)


##option 2: parameter 'a' is fixed and is not modelled in the fit
# m<-nls(y ~ fit(x,a=c(a00[1], a00[2]),b), start = list(b=c(b00[1],b00[2])),
#        control=nls.control(maxiter=50, tol=1e-5, minFactor = 1/2048, warnOnly=TRUE),
#        trace=TRUE)

# ##option 3: Does not work! But this is what needs to be corrected. 
# I want to be able to fix a[2], but fit a[1]
m<-nls(y ~ fit(x,c(a00[1], a00[2]),b), start = list(a=c(a00[1], a00[2]), b=c(b00[1],b00[2])),
       control=nls.control(maxiter=50, tol=1e-5, minFactor = 1/2048, warnOnly=TRUE),
       trace=TRUE)


plot(x,y)# Plot raw data
lines(x,predict(m),lty=1,col="red",lwd=1) #Add fitted line

...