PDM on R
PDM on R
<加藤太一さん[vsolj 9967] PDM on R 2010.03.09より>
しばらく書いていなかったのでご存じでない方もおられるかも知れません。 周期解析(特に変光星の)ならば、Pernaso なんかより (^^; R を使って PDM がおすすめ。京都のメンバーが使い込んでいます(^^; PDM は古いアルゴリズムなので、と新しいのを好まれる方もあると思いますが、 性質が非常に自然で多くの方に好まれています(最近出た "Asteroseismology" にもそう書いてあります)。含まれている周期成分の探索よりも精密周期決定 向けの方法です。 R は最初はちょっととっつきにくいかも知れませんが、データ処理は大変便利です。 以下にコードを挙げておきます。 === ごとに分割してください。
=== (pdmr.c)
#include <stdio.h> #include <math.h> #include <stdlib.h> static void initsum(int nb, int *ns, double *sum, double *ssum); static double phi(const double *jd, const double *mag, int n, double per, double sigma, int nbin, int cover, int *ns, double *sum, double *ssum); void pdmr(const double *jd, const double *mag, const int *n, const double *per, double *theta, const int *nper, const int *nbin, const int *cover, int *ns, double *sum, double *ssum) { int i; double s = 0; double ss = 0; double sigma; if (*n <= 0 || *nbin <= 0 || *cover <= 0) return; for (i=0; i < *n; i++) { double y = mag[i]; s += y; ss += y*y; } sigma = (ss-s*s/(*n))/(*n-1); for (i=0; i < *nper; i++) { theta[i] = phi(jd,mag,*n,per[i],sigma,*nbin,*cover,ns,sum,ssum); } } static double phi(const double *jd, const double *mag, int n, double per, double sigma, int nbin, int cover, int *ns, double *sum, double *ssum) { int i; double sigmasum = 0; int ndata = 0; initsum(nbin*cover,ns,sum,ssum); for (i=0; i<n; i++) { double phase; double val; int base; int j; phase = jd[i]/per; phase -= floor(phase); base = (int)floor(phase*nbin*cover); val = mag[i]; for (j=0; j < cover; j++) { int k = (base+j) % (nbin*cover); ns[k]++; sum[k] += val; ssum[k] += val*val; } } for (i=0; i<nbin*cover; i++) { int pn = ns[i]; double ps = sum[i]; double pss = ssum[i]; if (pn > 0) { sigmasum += (pss - ps*ps/pn); ndata += pn; } } return sigmasum/(ndata-nbin*cover)/sigma; } static void initsum(int nb, int *ns, double *sum, double *ssum) { int i; for (i=0; i<nb; i++) { ns[i] = 0; sum[i] = 0; ssum[i] = 0; } }
=== (pdminit.R)
dyn.load("pdmr.so") source("PDM.R") # ここは必要に応じて絶対パスを記述ください
=== (PDM.R)
source("lfit.R") # ここは必要に応じて絶対パスを記述ください PDM <- function(dat, minp, maxp, fineness, freq=FALSE, nbin=25, cover=5) { jd <- dat$V1 mag <- dat$V2 fineness <- fineness + 1 ndata <- length(jd) if (freq) { per <- seq(1/maxp, 1/minp, length=fineness) per <- 1/per } else { per <- seq(minp, maxp, length=fineness) } theta <- double(fineness) pdm <- .C("pdmr", jd = as.double(jd), mag = as.double(mag), n = as.integer(ndata), per = as.double(per), theta = as.double(theta), nper = as.integer(fineness), nbin = as.integer(nbin), cover = as.integer(cover), ns = integer(nbin*cover), sum = double(nbin*cover), ssum = double(nbin*cover) ) if (freq) r <- list(pdm=data.frame(V1=1/per, V2=pdm$theta),ndata=ndata) else r <- list(pdm=data.frame(V1=per, V2=pdm$theta),ndata=ndata) class(r) <- "PDM" return(r) } list.minima <- function(pdm, num=1, dalias=1,prec=7) { minper <- rep(0,num) fit <- rep(0,num) pdmpdm <- pdm$pdm fmt <- sprintf("%%.%df ± %%.%df",prec,prec) for (i in 1:num) { pos <- which.min(pdmpdm$V2) minper[i] <- p <- pdmpdm$V1[pos] avoid <- (1/(1/p-1) - p)*0.7/dalias pdmpdm$V2[abs(pdmpdm$V1-p) < avoid] <- 2 e <- minerr(pdm,p,pdm$ndata) fit[i] <- e$minper print(sprintf(fmt,e$minper,e$err)) } return(list(minper=minper,fit=fit)) } minerr <- function(pdm, per, ndata) { pdmpdm <- pdm$pdm pos <- which.min(abs(pdmpdm$V1-per)) aicmin <- Inf xx <- NULL for (i in 3:30) { aa <- pdmpdm[(pos-i):(pos+i),] ma <- mean(aa$V1) aa$V1 <- aa$V1 - ma xxold <- xx xx <- lm(aa$V2 ~ poly(aa$V1,2,raw=TRUE)) aic <- AIC(xx) if (aic > aicmin) break aicmin <- aic } ce <- as.vector(coef(xxold)) a <- ce[3] b <- ce[2] c <- ce[1] m <- ma - b/(2*a) e <- sqrt(abs((4*a*(c-0)-b*b)/(4*a*a*(ndata-1)))) return(list(minper=m,err=e,range=i-1)) } PDM.pf <- function(phi,ndata,nbin=25,cover=5) { return(pf(phi,ndata-1,ndata*cover-nbin*cover)) } PDM.qf <- function(phi,ndata,nbin=25,cover=5) { return(qf(phi,ndata-1,ndata*cover-nbin*cover)) } plot.PDM <- function(pdm,ylim=NULL,xlab="Period",ylab=expression(theta),...) { yl <- ylim pdmpdm <- pdm$pdm if (is.null(ylim)) { ymax <- max(pdmpdm$V2) ymin <- min(pdmpdm$V2) y1 <- ymin - (ymax-ymin)*0.15 y2 <- max(1,ymax + (ymax-ymin)*0.15) yl <- c(y1,y2) } plot(x=pdmpdm$V1, y=pdmpdm$V2, ylim = yl, typ="l", xlab = xlab, ylab = ylab, ...) }
=== (lfit.R)
lfit <- function(dat) { xm <- mean(dat$V1) x <- dat$V1 - xm res <- lm(dat$V2 ~ x) dat$V3 <- dat$V2 dat$V2 <- res$residuals return(dat) } lfit0 <- function(dat) { m <- mean(dat$V2) dat$V3 <- dat$V2 dat$V2 <- dat$V2-m return(dat) } lfitn <- function(dat,n) { xm <- mean(dat$V1) x <- dat$V1 - xm res <- lm(dat$V2 ~ poly(x,n)) dat$V3 <- dat$V2 dat$V2 <- res$residuals return(dat) } lfit.param <- function(dat) { m <- lm(dat$V2 ~ dat$V1) c <- summary(m)$coeff print(sprintf("%.6f %.6f",c[1,1],c[1,2])) print(sprintf("%.10f %.10f",c[2,1],c[2,2])) }
===
使い方 (1) install R は最近の Linux だと R-core のパッケージを選んでインストールすれば 簡単に入るようです。 OS のプロンプトから R CMD SHLIB pdmr.c エラーが出ればライブラリが足りないなど、適宜インストールしてください。 (コンパイルするような場合 R-devel が入っていないことがある) (2) 実行 R のプロンプトから > source("pdminit.R") > d <- read.table("data") (JD, mag の並んだデータ) > pdm <- PDM(d,0.05,0.08,1000) > plot(pdm) などすれば図が出ます。 > list.minima(pdm) で極小値を求めます。(このあたりはあまり修正していなくて古いままです) データのトレンドを除くには d <- lfit(d) # 1次関数 d <- lditn(d,2) # 2次関数 など データの一部を使うには t <- subset(d,V1 > 55200) など
Keyword(s):[PDM] [R言語] [周期解析]
References: