Cannot fit function in R

Question

I wanted to make a quick test:

1. Create a function having lots of points. Add some random noise to the data.

2. For every n points, calculate the mean and sd.

3. Display the new ("blocked") data using error bars (which are something like sd*1.96/sqrt(n))

4. Fit the "blocked" data using nls fit and weights (which are (sd/sqrt(n))^(2))

I then get an error on the fit Error in .swts * attr(rhs, "gradient") : ....

5. without weights, this just works nicely.

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

Your weights are too large. Normalize them. This is equivalent without using such large weights.

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

giving:

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