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) # １次関数
d <- lditn(d,2) # ２次関数　など

　データの一部を使うには

t <- subset(d,V1 > 55200) など