PFilter 周期をふるいにかける

先週、OQ Gemの二回目の極小を観測するまでに幾つもの食外の観測がありました。（Mdyさんも何度か手伝って下さいました。）

そこで、PFinderが出した良さそうな周期を使って何十も位相図を作り、今まで観測した食外のデータが極小付近に来ないかどうかで周期をふるい分けました。かなり大変な仕事なので、それをスクリプトでできないか考えて、PFilterというのを作りました。

これは今までの食外の観測の開始時刻と終了時刻をリストにしたマトリックスと候補の周期の集合をインプットして、その周期集合をふるいにかけるものです。

細かいことはまた書きます。

結果を載せますと、まず周期集合

これをふるいにかけました。

ほとんど全滅で困り、相当のストレスでした。

そうこうしていたら１８日夜の観測で極小が引っかかったのです。

これでいろんな基礎データが大きく変わりました。まず極小データが一つ増え、周期の可能性が、僕の二回の極小の間隔28.8688日の1/1、1/1.5、1/2、1/2.5、、、1/14.5という２８個だけのの数列になりました。（周期が2日未満は不可とわかっています。）

PFinderと同じような考えでその中から良さそうな周期を選択するスクリプトPSelectorを作り、周期を選ばせたら上位から３つが

となりました。これは周期集合がたった28個なので瞬時に答えが出ます。

このスクリプトは手作業で周期7.2160日を見つけた後で作ったのですが、この結果をみてだいぶ安心しました。

その計算を載せておきます。

> div <- seq(1,14.5,0.5)

> P <- 28.8688/div
>
> n <- length(P)
> pow <- vector("numeric",n)
>
> for (i in 1:n) {
+ N <- round(d / P[i])
+ res <- lm(d ~ N)
+ P[i] <- coef(res)[2]
+ dev <- deviance(res)
+ pow[i] <- (P[i]^2)* m / dev
+ }
>
> max(pow)
[1] 281.268739107
> P0 <- P[which.max(pow)]
> N <- round(d /P0)
> res <- lm(d ~ N)
> coef(res) # epoch & Period
(Intercept) N
0.203144855046 7.216070079084 　選ばれた周期
> residuals(res) # residuals for each minimum
1 2 3 4
-0.2031448550418 0.0144779481667 0.9367946040411 -0.1907731464772
5 6
-0.2809371170861 -0.2764174336027 　残差です。
>
> ord <- order(pow, decreasing=T)
> P[ord[1:3]] # the best 3 periods
[1] 7.21607007908 4.12424539395 3.20775292472

やはり３番目の極小のずれが0.9367..と大きいです。

３番目の極小データを取り除いて同じことをさせるとPが少し減り7.2160056..となって手仕事で見つけた7.2160と同じになります。

> coef(res) # epoch & Period
(Intercept) N
0.11360649727 7.21600564715 　　選ばれた周期。
> residuals(res) # residuals for each minimum
1 2 3 4
-0.1136064972650 0.1098796119930 0.0358763680185 -0.0184634469843
5
-0.0136860357622 　　残差です。
>
> ord <- order(pow, decreasing=T)
> P[ord[1:3]] # the best 3 periods
[1] 7.21600564715 28.86955769417 3.60871702552