muttenz's blog

スイス星空だより

R上でC言語のPDMを動かすには。

以前にも書きましたが、基本的にはIakさんの記事にしたがってやってみました。

 

参照記事

[vsolj 9967] PDM on Rの記事 2010年3月9日

Iakさんの記事

http://blog.goo.ne.jp/imako009/e/5d80ee68888c948637c96566d4ba73d4

http://blog.goo.ne.jp/imako009/e/a200a686eba03b31ae13032889d2f31c

Meinekoさんの記事

http://d.hatena.ne.jp/meineko/20110816

 

R環境はもうコンピュータにインストールされているとします。僕のはR x64 2.15.0のヴァージョンです。

 

Iakさんの記事のようにたどっても良いのですが、そこに述べられている3つのプログラムをインストールするには以下のようにします。

ただし、MinGWを僕はどうもインストールしていないようです。それでもPDMは動いているからPDMを動かすだけならこのプログラムは必要ないかもしれません

1.Rtoolsと2.Perlだけインストールして、うまく動かなかったら3.MinGWも入れてみてください。

 

1. Rtoolsをダウンロードするためにはこれをクリック。

http://cran.md.tsukuba.ac.jp/bin/windows/Rtools/

ここでDownloadの表の一番上の Rtools30.exeをダウンロードします。

これをインストールする時にいろいろな指示が出ますが、その途中でパスを書き込ませますか?と言うのが出たら、書き込ませるようにすると、自分でやる必要がなく安心です。

2.Active Perlをダウンロードするためにはこれをクリック。

http://www.activestate.com/activeperl

するとDownload now”のボタンがありますから、それをクリック。二つの選択肢がでます。僕の場合は64ビット用を選びました。

これもインストールする時にいろいろな指示が出ますが、その途中でパスを書き込ませますか?と言うのが出たら、書き込ませるようにすると、自分でやる必要がなく安心です。

(3.MinGWをダウンロードするためにはここをクリック 

http://sourceforge.net/projects/mingw/files/

画面中央左にLooking for the latest version?という文章の後ろに

Download mingw-get-inst-20120426.exe

が書いてありますから、これをクリックしてダウンロード。

右の方にDownloadと言うボタンがあったりしますが、これは関係ないソフトのダウンロードです。これをクリックしないように!これを僕はクリックしてひどい目にあいました)

これでプログラムのインストールとパスは準備完了。

 

さらにRにもパスを通しておきます。

(僕の場合は)「C:\Program Files\R\R-2.15.0\bin\x64 」をPath 変数の中に書き込む。(パスの通し方はGoogleで引くといろいろ出てきます。Windows7の場合は次のサイトに説明があります。

http://blog.cnu.jp/2009/11/06/windows-7-path/

 

これから、Iakさんが書いておられるようにMLにあるプログラムを分けて次の4つのファイルを作る。(ここからはIakさんの手順をほとんどそのまま採っています)

pdmr.c

PDM.R

pdminit.R

lfit.R

(僕はこの4つのファイルをRのWorkフォルダに入れました。)

それから(アクセサリーから)コマンドプロンプトを開き、pdmr.cのファイルが置いてあるディレクトリに行って、そこで次のコマンドを実行する。

R CMD SHLIB pdmr.c

すると

pdmr.dll

pdmr.o

という二つのファイルがそのディレクトリの中にできる。

つぎにpdminit.Rの内容を少し変更する。

dyn.load("pdmr.dll")

source("PDM.R")

 

これで準備万端整いました。

 

Rをスタートして、Rでのコマンドプロンプトで以下のようにして動けばばんざーい。(ここからはMLの丸写しです)

> 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) など

 

注ですが

pdm <- PDM(d,0.05,0.08,1000)で、0.05が周期を探す下限、0.08が上限で、1000がその区間をいくつに分割して探すかの値です。

また、フェーズの図表示 "plot(phase())"やフェーズをいくつかの区間に分けてその区間の平均値でフェーズの図を表示させる"plot(phave())"にはPergrmのスクリプトをさらに動かしてからやっています。

pergrmがどこにあったかこれから探してみます。どこでしたっけ?)

 pergrmなどに関して、meinekoさんが詳しく述べておられるので、そちらを参照してください。

http://d.hatena.ne.jp/meineko/20110818