gnuplot

2カラムのデータが入ったファイルのデータをプロットするとき。

plot "ファイル名" 1:2

"< UNIXコマンド名"の標準出力も解釈できる。

batchモードのR

とりあえずコマンドラインで入力するコマンドを書き込んだファイルをつくって

R --slave --vanilla <**.r

で動くのでそのままやっていたが引数を渡したくなった。こうすればいいらしい。Rスクリプトファイルで

args <- commandArgs()

とし

R --slave --vanilla --args a b c < **.r

commandArgsはどうも全引数を取得するので、第一要素がRのパス、二つ目が--slave、三つめが--vanilla、四つ目が--argsで、五つめ以降にa,b,cというベクトルがargsに代入されることになる。一瞬--argsは必要ないのかと思ったが、--argsを使わないとコマンドに対して無意味な引数を渡せないのでこれは必要。

library(lattice)

層別ヒストグラムを書きたいとき

library(lattice)

histogram(~変数名|層名1+層名2,データ名)

texファイルのレンダリング

いつも忘れてしまうので。hoge.texを作成したら

platex hoge.tex
dvipdfmx hoge

図を入れるときは、プリアンブルに

\usepackage{graphicx}

として、本文中で

\includegraphics[scale=1.5]{sin.ps}

とかする。

Perk/Tkでtableウィジェット

ここを参考にした。新たなモジュールが必要。

Use Tk::Table;

ウィジェットの定義は

my $table_frame = $mw->Frame()->pack();
my $table = $table_frame->Table(-columns => 8,
-rows => 4,
-fixedrows => 1,
-scrollbars => 'oe',
-relief => 'raised');

という感じ。オプションの意味はまだ全部は分からず。scrollbarsはnが上、sが下、wが左、eが右にスクロールバーをつける。reliefは形状、flat (平坦)、raised (出っぱり)、sunken (引っ込み)、groove (溝)、 ridge (土手)となっている。

で、タイトル行をまず入力するには

foreach my $col (1 .. 8)
{
my $tmp_label = $table->Label(-text => "COL " . $col, -width => 8, -relief =>'raised');
$table->put(0, $col, $tmp_label);
}

テキストデータの場合、labelウィジェットを作成しputで入れていく感じになるようだ。次に内容のデータは

foreach my $row (1 .. 8)
{
foreach my $col (1 .. 8)
{
my $tmp_label = $table->Label(-text => $row . "," . $col,
-padx => 2,
-anchor => 'w',
-background => 'white',
-relief => "groove");
$table->put($row, $col, $tmp_label);
}
}
$table->pack();

最後のpack()で反映、はいつも通り。あとputしてpackしていけば追加していけるが、最初に指定した行数をこえて追加したい時はどうすればよいのか、まだわからない。

Tableを空っぽにするには、

$table->clear

Tableの各cellの選択、コピペはできないようだ。

いろいろ

prop.tableは、count tableをfrequencyにするとき役立つ。このときprop.table(..,margin=x)だけど、x=0なら全部わる。x=1なら一つ目の変数？ごとにわる。x=c(1:2)なら、1つめ、2つめでわる。という感じでやってくれるようだ。助かるなあ。

xtabsはcrossのtableを作ってくれる。xtabs(~A+B+C,data=x)的な書き方で。

library(lattice)のbarchartをこれらと組み合わせると超いい。

ロジスティックの説明変数の検定とか

ポイント！full modelの全係数=0の検定が強い関連を示しているのに、各係数の検定があまり対した関連を示さない場合、multicolinearityの存在が示唆される！のだそうな

Waldの検定

t or z=β/SEがt分布または標準正規分布にしたがうことを利用するか、またはz²がdf=1のカイ二乗分布にしたがうことを利用する。

尤度比検定

モデルの対数尤度は、以下のようにして求める。

>logLik(glm(Y~x1,data=x,family=binomial(link="logit")))
'log Lik.' -97.22633 (df=2)

変数なしモデルの対数尤度は次のようにすればよい。

>logLik(glm(Y~1,data=x,family=binomial(link="logit")))
'log Lik.' -112.8793 (df=1)

検定はこんな感じで。

>1-pchisq(-2*(-112.8793+97.23),1)
[1] 2.212391e-08

スコア検定

度数分布表

次のようにする。

table(cut(x,seq(0,100,10)))

誤差つき棒グラフ

>library(gregmisc)
>barplot2(x,plot.ci=T,ci.l=x-se,ci.u=x+se

という感じで

Cygwinでscreen

で、あるが、どうもデタッチができないように思っていた。一度teratermを終了すると、デタッチしていたscreenがdeadになる。某サイトにもそう書いてあった。しかし、「コンピュータ」のプロパティから環境変数「CYGWIN」を作成し、値を「tty」とすると・・・とりあえずteratermを終了し再起動しても、アタッチできた。これでPC再起動後もアタッチできればよいのだけど。

作図tips

軸の数値なしのグラフをplotで書きたいときは、xaxt="n"もしくはyaxt="n"と指定する。

その後、axis(side,at=..)で好きなように軸をかける。sideは1=下、2=左、3=上、4=右。atは数値でも文字ベクトルでもよい。文字ラベルにしたいときはlabels=..とする。

たとえば、x軸が28000000から34000000とかだと、Rは自動的に軸注釈を2.8e+7.0とかにしてしまう。これを避けるための方法はみつからなかったので、xaxt="n"とし、

axis(1,at=seq(2.8e+7,3.4e+7,1e+6),labels=paste(28:34,"000000",sep=""))

とすれば指数表記でない軸が作成可能。

右や上軸にしても、xlabは左にしか軸ラベルを描画しない。そこでこの場合は作図領域外に文字を書くmtext(text,side,line,outer=T)を使う。sideは上と同じで、lineは領域外何行目に書くか。x軸ラベルの場合、line=1としてちょうどよい。

軸注釈の大きさを変えるのはcex.axis。

Cygwin文字化け対策の覚書

メッセージの文字化けについて

シェルのメッセージの文字化け

export OUTPUT_CHASET=sjis

これではうまくいかないと書いてあるサイトもあるけど、いちおううまくいっているみたい

vimの文字化け

まだ解決していない

lessの文字化け

lvを使えばいいらしい

モンテカルロ法

正規化定数の不明な分布P(x)が与えられている。統計量A(x)の分布P(x)のもとでの期待値を求めたい。

定数Zについて

P(x)=\frac{\tilde{P}(x)}{Z}

となる~P(x)が書け、これが容易に計算できるものであるとする。

また、P(x)によく似た分布Q(x)があるとする。

1. Q(x)からサンプルx_αを生成。
   
2. A(x_α)と

   \omega^{(\alpha)}=\frac{\tilde{P}(x^{(\alpha)})}{Q(x^{(\alpha)})}

   を計算する。

期待値は、

\sum_x A(x)P(x)=\frac{\sum_{\alpha}A(x^{(\alpha)})\omega^{(\alpha)}}{\sum_{\alpha} \omega^{(\alpha)}}

具体例

円の面積を求めるモンテカルロ法では、P(x)が単位円の一様密度、Q(x)が外接する一辺2の正方形の中の一様密度、円の面積がZ。半径1の円の中にxがあるなら~P(x)=1、そうでなければ~P(x)=0、Q(x)も同様。

1. 正方形の中の一点x_αをサンプリングする
   
2. ω_α=~P(x)/Q(x)
   
3. 連続変数なので期待値は

わっかんないな～

Box-Cox変換

なかなか正規化しない分布を無理やり正規化する方法。
λを用いてこんな変換をする

T(x)=\frac{x^{\lambda}-1}{\lambda}

MASSライブラリをロードし、

library(MASS)
box_par <- boxcox(y~x,lambda = seq(-0.25, 1, length = 100))

などとするとλの尤度が計算される。

p <- box_par$x[which.max(box_par$y)]

で、最尤推定量を取り出す。

TeX

これもすっかり忘れていたので。

platex **.tex
dvipdfmx **.dvi

って感じ

bloggerで数式

下の投稿で試したけど、ダッシュボード->レイアウト->HTMLの編集、で

<script type="text/javascript" src="http://tex.yourequations.com/"></script>

を、</body>の直前に挿入する。

そうすると、bloggerの投稿としてはHTMLの編集のところで

<pre lang="eq.latex">
f(x|\mu,\sigma)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^2}{\sigma^2}}
</pre>

とかやると、

f(x|\mu,\sigma)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^2}{2\sigma^2}}

となる。素晴らしい。

ちなみにこれに気付く前はクリボウのページにあったスクリプトを使おうとしたんだけど、mimetexのサーバ使用制限でこれはすでに使えない。

混合正規分布

いまのところまったくわかっていません。とりあえず混合比

p_r

、成分分布（、というのは一般化した場合で、この場合正規分布）

h(x,\omega_r)

として

p_{mix}(x,\theta)=\sum_{r=1}^{R}p_r h(x,\omega_r)

と表わされるんだそうだ。そんで、ダミー変数Zを使って

Prob(Z=r)=p_r

X|Z=r \sim h(x,\omega_r)

(X,Z=r) \sim P_r h(x,\omega_r)

とすると、xが与えられた時のZ=rの条件付き密度が

p(Z=r|x,\theta)=\frac{p_r h(x,\omega_r)}{p_{mix}(x,\theta)}

となって、これがEMアルゴリズムに利用されるんだって。

RでCを使うとき

完全に忘れていたので覚書。ヘッダは

#include<R.h>
#include<Rdefines.h>

で、あと必要なもの。Rのベクトル変数はSEXPで指定。関数宣言は、

SEXP likelihood(SEXP InitComplex,SEXP scn,SEXP ssn1,SEXP ssn2,SEXP theta);

とすれば返り値もベクトルになる。このRのベクトル変数を引数でもらったときはメモリ確保が必要。

SEXP p;

PROTECT(ssn1=AS_INTEGER(ssn1));
PROTECT(ssn2=AS_INTEGER(ssn2));
PROTECT(theta=AS_NUMERIC(theta));
PROTECT(p=allocVector(REALSXP,muallele*2));

最後のは関数中で宣言したR形式ベクトル変数の領域確保。

これはreturnする前にかならず解放する。

UNPROTECT(4);
return p;

ソース作成したらコンパイルは、

R CMD SHLIB **.c

で、.soファイルと.oファイルができる。
Rからの読み込みは、

dyn.load("**.so")

で、関数呼び出しは、

.Call("関数名",引数)

とする。この場合の引数はRらしからぬ厳密な型（モード？）宣言が必要。as.numericなどを使用する。

Metropolisアルゴリズム

まだキモはわかっていないが、A-Rアルゴリズムからの発展系みたいなもんのよう。また、Importance samplingはMetropolis Hastingsに対応するようだ。

確率密度関数π(x)からサンプリングしたい。ここで対称サンプラーq(x,y)=q(y,x)を利用する。

1. 初期値X₀を決める。
2. q(Y|X₀)からサンプルYを取りだす。
3. u~U[0,1]を発生。
4. u < π(Y)/π(X₀)ならX₁=Y、そうでなければX₁=X₀

この操作を繰り返して得たX₁,..,X_nはπ(x)からのサンプルである。

u < \frac{\pi(Y)}{\pi(X_0)}

を基準としているが、これが

1 < \frac{\pi(Y)}{\pi(X_0)}

なら逐次改良法となる。メトロポリス法では、かならずしも新しいモデルが古いモデルより確率が高くなくても、確率uで採用するという方法である。

local maximaから逃れられる方法であることは理解できるが・・・

A-R サンプリングの例

たとえばπ(x)=N(1,3)からのサンプリングはこんな感じ。h(x)=1/60、c=10としてる。なんで？ch(x)=0.167となる。これがN(1,3)の最大値0.133より大きいから。c=8くらいでもよさそうだけど。とりあえずh(x)は均一分布の確率密度関数(-30~30の範囲について）になる。

u~U(0,1)

Z~h(x)=U(-30,30)

について、

u<=f(Z)/ch(Z)

にてサンプリングされたものが赤丸で表わされている。棄却されたものは黒。10000のZに対し、998のサンプルがacceptされた。

たとえばh(x)=1/40、c=7とするなら1412とかのサンプルがacceptされたりする。

いずれにせよA-Rサンプリングでは棄却が多すぎて非効率的なこと、またch(x)が求められない場合があることから、MCMCを行うことになるということらしい。

acceptance-rejection (A-R) サンプリング

確率密度関数π(x)=f(x)/Kからサンプリングしたい。f(x)はunnormalized (非標準？)分布で、Kはnormalizing constant。f(x)は離散分布でも連続分布でもよい。

ここでf(x)<=ch(x)なるcとh(x)をとる。

アルゴリズム

1. Z~h(・)、u~U(0,1)をとる。
2. u<=f(Z)/ch(Z)ならZ=Xを返す。
3. そうでなければ元に戻る。

2がacceptance、3がrejection。意味的には、横にZ、縦にuの座標を取り各Zにおいてランダムにとったu<=f(Z)/ch(Z)となる点を返す(f(x)/ch(x)<=1なので、点が返されないZはない)。

意味は？

P(X≦x)=P(Z≦x|u≦f(Z)/ch(Z))

ところでP(Z≦x)=∫-∞~x h(y)dyで、

P(Z≦x, u≦f(Z)/ch(Z))=∫ f(y)/ch(y) ・ h(y)dy=(K/c) ∫-∞~x π(y)dy

ここでx→∞とすると

P(u≦f(Z)/ch(Z))=K/c

なので

P(X≦x)=P(Z≦x|u≦f(Z)/ch(Z))=P(Z≦x, u≦f(Z)/ch(Z))/P(u≦f(Z)/ch(Z))=∫-∞~x π(y)dy

というわけでZは確率密度関数f(Z)/cからのランダムなサンプリングとなる。