Не может вписать функцию в R - PullRequest
Новая
       45

Не может вписать функцию в R

0 голосов
/ 11 октября 2019

Я хотел сделать быстрый тест:

  1. Создать функцию, имеющую много точек. Добавьте к данным случайный шум.

  2. Для каждого n точек вычислите среднее значение и сд.

  3. Отображение новых («заблокированных») данных с использованием полос ошибок (например, sd * 1,96 / sqrt (n))

  4. Подгоните «заблокированные» данные, используя nls fit и weights (которые (sd / sqrt (n)) ^ (2))

Я тогда получаю ошибку при подгонке Ошибка в .swts * attr (rhs, «градиент»): ....

без весов, это просто прекрасно работает. 
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

1 Ответ

0 голосов
/ 11 октября 2019

Ваш вес слишком велик. Нормализуй их. Это эквивалентно без использования таких больших весов. 

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