簡単なテストをしたかった：
多くのポイントを持つ関数を作成します。データにランダム ノイズを追加します。
n ポイントごとに、平均と標準偏差を計算します。
エラーバー (sd*1.96/sqrt(n) のようなもの) を使用して、新しい (「ブロックされた」) データを表示します。
nls フィットと重み ((sd/sqrt(n))^(2)) を使用して、「ブロックされた」データをフィットします。
次に 、.swts * attr(rhs, "gradient") の fit Error でエラーが発生します: ....
- 重みがなければ、これはうまく機能します。
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