无法在 R 中拟合函数

Question

我想做一个快速测试：

1. 创建一个有很多点的函数。向数据中添加一些随机噪声。

2. 对于每 n 个点，计算平均值和 sd。

3. 使用误差线（类似于 sd*1.96/sqrt(n)）显示新的（“阻塞”）数据

4. 使用 nls 拟合和权重（即 (sd/sqrt(n))^(2)）拟合“阻塞”数据

然后我得到一个适合的错误 .swts * attr(rhs, "gradient") 中的错误：....

5. 没有重量，这很好用。

par(mfrow=c(1,1))
# sampling frequency
interval<-0.05
# max frequency
max_F<-200
Freq<-seq(1,max_F,by=interval)
# define Lorentzian Function parameters
a1<-2
a2<-3
a3<-60
# Lorentzian function
y<- a1/((Freq-a2)^2 + a3^2)
# add some noise
c.norm <- rnorm(length(y))
amp <- 0.2
y_noise<-y+y*c.norm*amp
# block data (every n points)
n<-200
# The next few lines are to
# calculate the mean and sd values over n points
reminder_F<-length(Freq)%%n
length_N<-length(Freq)-reminder_F
Freq<-Freq[1:length_N]
y_noise<-y_noise[1:length_N]
temp_group<-(rep(1:(length_N/n), each=n))
Y <- data.frame(value =y_noise, group=temp_group )
# I use tapply over groups (each group is n points)
# for the mean value
mean_val<-tapply(Y$value, Y$group, FUN = mean)
# and the sd
sd_val<-tapply(Y$value, Y$group, FUN = sd)
# I then create a new x axis
Freq_N<-seq((n*interval)/2,max_F,by=(n*interval))
Freq_N<-Freq_N[1:length(mean_val)]
# a quick plot to make sure everything is ok
par(mfrow=c(1,1))
plot(Freq_N,mean_val,xlim=c(50,100),ylim=c(1E-4,3E-4))
par(new=TRUE)
plot(Freq,y_noise,xlim=c(50,100),ylim=c(1E-4,3E-4),  cex = .1,col="red")
grid()
par(new=FALSE)
# Now use ggplot to show error bars
# I use a log log plot (not mandatory here)
IC<-1.96*sd_val/sqrt(n)
Y_N=data.frame(Freq_N,mean_val,IC)

# Plot
p1<-ggplot(Y_N, aes(x=Freq_N,y=mean_val))+geom_point()+
  geom_errorbar(aes(ymin=mean_val-IC, ymax=mean_val+IC), width=.02,
                )+scale_x_log10()+scale_y_log10()
p1
# FiT
my_weights<-(sqrt(n)/sd_val)^2
myfit<-nls(formula= mean_val ~  a1/((Freq_N-a2)^2 + a3^2)
           , data=Y_N, start=list(a1=1,a2=1,a3=1),weights = my_weights )
myfit

Answer 1

你的重量太大了。标准化它们。这在不使用如此大的权重的情况下是等效的。

nls(formula= mean_val ~  a1/((Freq_N-a2)^2 + a3^2)
           , data=Y_N, start=list(a1=1,a2=1,a3=1),
           weights = c(prop.table(my_weights)) )

给予：

Nonlinear regression model
  model: mean_val ~ a1/((Freq_N - a2)^2 + a3^2)
   data: Y_N
     a1      a2      a3 
  2.021   2.029 -60.137 
 weighted residual sum-of-squares: 1.163e-12

Number of iterations to convergence: 7 
Achieved convergence tolerance: 1.219e-06