観察データの差分とかAUCとか

• 多変量データの周期性についていじっている
• そこでは、微分や積分を行き来する
• 少し立ち返って、単変量時系列データをRでハンドリングして、その微分(差分)や積分(Area under the curve計算)をすることについてみてみる
• こちらと関連しており、データ発生に借用させていただくこととする
• cumsum(),diff()を使ってみる
• 観察点をばらつかせながら増やしてやって、AUCが収束する様子が掲載図

AUC<- function(xvec,yvec){                      # 台形の面積を計算する。
  int<- 0
  for(i in 1:(length(yvec)-1)){
    int<- int + 
          (yvec[i]+yvec[i+1])*(xvec[i+1]-xvec[i])/2
  }
  return(int)
}
f<- function(x,a){                                       # この関数は適当
  p<- -x*(log(x)-(1+a))
  return(p)
}
T<- c(0,f(1:15,2))                                       # 15h後に服用するくらいの意味で
#T<- T[which(T>0)]
for(i in 1:3){                                           # 継続服用
  T<- append(T,f(1:15,2)+tail(T,1))
}
T<- append(T,f(1:1000,2)+tail(T,1))
T<- T[which(T>=0)]

#plot(T,type="b")
Nt<-length(T)

# 差分は微分みたいなものですが
SABUN<-diff(T)
# 差分から元のデータを作るのは「積分」
# 観察間隔が一定なら
Hukugen.Sekibun<-cumsum(c(T[1],SABUN))
plot(T,Hukugen.Sekibun)

# 観察間隔が一定でないとき
# 部分的な観察データを作る
# 観察時点数を1から全時点まで振る
obsFrac<-2:length(T)
# それごとにAUCを計算して格納する
AUCs<-rep(0,length(obsFrac))
AUCs2<-AUCs
for(i in obsFrac){
	# 観察時点数iを適当に選ぶ
	obsTimes<-sort(sample(1:Nt,i))
	Tobs<-T[obsTimes]

	plot(obsTimes,Tobs,type="b",xlim=c(0,length(T)),ylim=c(0,max(T)))
	# 観察時系列データの時刻の差分をとる
	TimeDiffs<-diff(obsTimes)
	# 観察時系列データの観測値の差分をとる
	TobsDiffs<-diff(Tobs)
	# 隣接観察
	TobsAverage<-Tobs[1:(length(Tobs)-1)]+TobsDiffs/2
	AUCs[i]<-sum(TobsAverage*TimeDiffs)
	AUCs2[i]<-AUC(obsTimes,Tobs)
}

plot(AUCs,AUCs2)
plot(AUCs)