Sunday, April 25, 2010

【Ubuntu】AucTexをインストール

UbuntuノートPCでTeXを使いたいので、ここを参考にしてAucTexというのをインストールする。


$ sudo apt-get install auctex



次の環境設定をしておく。

"~/.emacs"(".emacs.el"がある場合はそちらに)の最後に以下を書き加える:

;===================================
; AUCTeX
;===================================
(autoload 'japanese-latex-mode "tex-site" "tex-site" t)
(require 'tex-site)
(add-to-list 'auto-mode-alist '("\\.tex\\'" . japanese-latex-mode))
(setq TeX-default-mode ' japanese-latex-mode)
(add-hock 'TeX-mode-hook (function (lambda () (outline-minor-mode t))))


【音楽】ハーモニカ特訓その1

自分が練習しないといけないハーモニカの曲のキーはA。ということで、

音階は A(ラ)、B(シ)、C#(ド#)、D(レ)、E(ミ)、F#(ファ#)、G#(ソ#)、A(ラ)

となる。

このとき、キーがAのハーモニカの音の配列は、
12345678910
吹くA(ラ)C#(ド#)E(ミ)A(ラ)C#(ド#)E(ミ)A(ラ)C#(ド#)E(ミ)A(ラ)
吸うB(シ)E(ミ)G#(ソ#)B(シ)D(レ)F#(ファ#)G#(ソ#)B(シ)D(レ)F#(ファ#)
第1オクターブ第2オクターブ第3オクターブ

となる。

楽譜上での位置関係は

Saturday, April 24, 2010

【CSS】枠を作る


<div style="border:1px dashed gray; background-color:#FFE3E3">

</div>



これまでに何度かソースコードを枠で囲って記事を書いたことがある。その場合、機能によって文字を色分けしていた(例えば、''call''ってサブルーチンの場合は青色にするとか)。これは、ソースコードをSciTEで開いて、html出力して、ソースコードを貼り付けただけ。

【MPI】CとFortranとのインターフェイスの違い

 C言語とFortranではインターフェイスに違いが見られる。

サブルーチン
○C言語の場合

ierr = MPI_Xxxx(...);

○Fortranの場合

call mpi_xxxx(....,ierr)

システム用の配列
○C言語の場合

MPI_Status istatus;

○Fortranの場合

integer istatus(mpi_status_size)

データ型の指定
○C言語の場合
文字:MPI_CHAR
整数:MPI_INT
実数:MPI_FLOAT
倍精度実数:MPI_DOUBLE


○Fortranの場合
文字:MPI_CHARACTER
整数:MPI_INTEGER
実数:MPI_REAL
倍精度実数:MPI_DOUBLE_PRECISION
複素数:MPI_COMPLEX

Friday, April 23, 2010

ユークリッド空間をcountableな形に分解

 ユークリッド空間R^nがあるとする。このとき、

R^nはseparableである。

ちなみに、metric space Xがseparableであるとは、Xがcountable dense subset Bを含むということ。さらに、Xの部分集合Bがdenseであるとは、topological closure(Bを含む全ての閉集合のintersection)がXに等しくなること。

□DEFINITION (Berge1997 p.93)
A fundamental family for a topological space X is a family of open sets A={A_i:i∈I} such that, given any non-empty set G there exists a subset J of I for which

  G = ∪_{i∈J}A_i.

□THEOREM (Berge1997 p.93)
A metric space is separable if and only if it possesses a countable fundamental famility.

以上のことから、次のようなことが言える:

X⊆R^nとする。さらに、h:R^2→Rという関数があると、

 {y:sup_{x∈X}h(x,y) > b} = ∪_{i∈J}{y:sup_{x∈X_i}h(x,y) > b}.

ここで、Jはcountable set.

つまり、coutableな形に分解できるというメリットがある。テクニカルなところでこういう性質を使う場合があるかもしれない。

Thursday, April 22, 2010

ベイズ統計学のメモ2

 標本数が無限大のとき、すなわち大標本のとき、最尤法とベイズの事後分布の関係について。

 最尤法は、サンプルを所与として、likelihoodが最大になるようなパラメータを求める。そのとき、推定量の分布(母数=真のパラメータは定数であって、推定量が確率変数なだけ)は正規分布で近似できる。

 一方で、ベイズの事後分布はというと、まだ理解してない。ただ、最尤法との関係性についていうと、最尤法はlikehoodの頂点だけに着目しているのに対して、ベイズは全体を扱っている。つまり、最尤法は頂点(と頂点=推定量の分布)を見ているのに対して、ベイズは分布(の変化)を見ている。

 また、古典論では推定量の分布を扱っている。大標本になれば、一致性から真のパラメータに収束する。

 一方、ベイズ統計では、大標本になると母数の分布がある分布に収束する(はず)。この場合、パラメータは確率変数であって、ベイズは大標本だとその分布が収束すると言っているに過ぎない。つまり、ある特定の定数に収束するとは言ってない。

 ここに、両者の本質的な違いがある。

 ここから一つステップを踏んで、「では、ベイズ統計で得られたパラメータの事後分布の平均が、最尤法における推定量に一致するか?」というquestionを立ててみる。

 おそらく、一般的にはNOだろう。事後分布が正規分布となっていれば、頂点は平均に等しいのだから一致する。でも、そうじゃない場合、必ずしも頂点と平均が一致するとは限らない。

 だとしたら、より一般的なベイズの方が妥当じゃないかって思ってしまうのだが、まだ確信を持てないのでひとまず保留にしておこう。

 てか、混乱していて、全く理解できてないな(笑)

 一つの評価方法として考えられるのは、まず、仮定として「真のパラメータ(定数)が存在する」を設ける。そのとき、最尤法を使えば、必ず(適当な仮定の下で)真のパラメータに収束する。このとき、「じゃあベイズさんはどうなりますか?事後分布の期待値が真のパラメータにちゃんと収束しますか?」という問を立てればOK。

 これは古典論をまず基準として評価する方法。でも、そもそも、古典論の「母数は定数」であるって仮定がおかしい!という哲学的な立場があるとすれば、この評価方法はあまりのもバイアスがかかっている。じゃあ、これとは逆方向から評価する方法を考えるとすればどうなるか。それは、一般にはない。そもそも確率変数と定数を比較することができない。では、「古典論の推定量がベイズの事後分布の期待値に一致するか?」という問を立てるのはどうか。でもこの場合、だから何なの?ってなる。

 ようは、まずは母数が定数であるか確率変数であるか、という立場を明確にする必要がある。仮に、定数であるとする。この場合、ベイズの事後分布は分析者の予想を表しているに過ぎない。


【追記】
仮に母数が定数であるというのが真だとする。この場合、明らかに古典論は頑健。一致性などを示せているならば。

Sunday, April 18, 2010

Intel Fortran Compiler for Linuxを利用する

●利用するときは、

<コンパイラがインストールされているフォルダ> ./bin/ifortvars.sh ia32

と端末で入力。''ia32''の部分はアーキテクチャによって変更する。ただし、これによると、ログアウトするまで有効になるだけで、ログイン時は毎回この設定を行う必要がある。それが面倒な場合は、ホームディレクトリにある''.bashrc''ファイルの一番最後の行に

. /opt/intel/Compiler/11.1/069/bin/ifortvars.sh ia32

を書き加えればOK。


●ソース・コードのコンパイル方法は

$ ifort filename.f90

となる。いろいろとオプションをつけることができるが、それは勉強しよう。

 上のようにコンパイルすると、''a.out''という実行ファイルが作成される。これを実行するには端末で

$ ./a.out

と入力すればよい。

Saturday, April 17, 2010

Dell Inspiron mini 10v: Ubuntuを8.04から9.04にアップグレード


 Dell Inspiron mini 10vが届いた。注文してから約1週間。

 Intel Fortran Compiler 11.1をインストールしようとしたら、Ubuntu 8.04では対応してないとエラーが出た。サポートしているのはUbuntu 9.04。ということで、Ubuntu 9.04をインストール必要が出てきた。

 そこで、ここを参考にしてUSBブートのインストールをしようと思う。

●Intel Fortran Compilerをダウンロード。これを参考にしてインストールを進めていたが、次のようなエラーメッセージが表示される:

ステップ: 4 / 7 | インストール設定 - 重要な必要条件の不足
--------------------------------------------------------------------------------
1 つまたは複数の重大な未解決問題が原因でインストールを続行できません。
インストールを終了せずに、問題を解決して再度チェックすることができます。
インストールを終了して問題を解決し、再度インストールを実行することもできます。
--------------------------------------------------------------------------------
重要な必要条件の不足
-- 不明なシステムコマンド
--------------------------------------------------------------------------------
1. 問題の詳細情報を表示する [デフォルト]
2. 必要条件を再度チェックする

h. ヘルプ
b. 前のメニューに戻る
q. 中止
--------------------------------------------------------------------------------
オプションを選択するか、Enter を押してデフォルトを選択してください。 [1]:


●そこで、こんなのを発見。

$ sudo apt-get install libstdc++5 g++ binutils sun-java6-jre

で''libstdc++5'',''g++'',''binutils'',''sun-java6-jre''というパッケージをインストールしたら、上記の問題は解決。

非線形方程式の解法:特殊なケース

 今、次のような2変数の非線形方程式を数値計算で解くことを考えているとする:

 f_1(x_1, x_2) = 0
 f_2(x_1, x_2) = 0.

さらに、変数がとりうる値の範囲が次のようになっているとする:

 l_1 ≦ x_1 ≦ u_1
 l_2 ≦ x_2 ≦ u_2.

 内点解が保証され、かつf_1とf_2が連続微分可能であれば、Newton-Raphson法を修正して問題を解くことができるはず。

 では、もし、微分可能性が満たされない場合はどうだろうか。この場合はSecant methodを使うことが考えられる。さらに、特殊な場合として、区間[l_2,u_2]が非常に小さいとする。
 このとき、2変数のシステムにsecant methodを適用するのではなく、x_1についてだけSecant methodを適用し、x_2については単純に加重平均でアップデートするという方法が考えられる。f_1とf_2が連続であれば、根の付近においてはそのような近似が利用できると考えられる。

非線形方程式を解く場合は次元を下げて

 ''f(x) = 0''を満たすようなn次元ユークリッド空間内の点xを求める方法について、Newton-Raphson methodやMultivariate Secant methodなどの方法がある。
 次元nが大きくなればなるほど、計算時間が掛かるだけでなく、iterationが収束しにくいということが予想される。後者について、システムの均衡条件を使えば実は次元が2になるのに、それに気が付かずに次元が4のままで計算してしまうということもあったりする。特に、次元が大きいままでは収束しなかったものが、次元を低くすれば収束したというような場合ではかなり大きな違いが出てくる。均衡の条件を使うことでより効率的に数値計算を行うことができる可能性が示唆される。

Thursday, April 15, 2010

Intel Visual FortranでLAPACKを使う方法

 Intel Visual FortranでMKLのFortran95インターフェイス仕様のLAPACKを使う方法についてメモしておく。マニュアルを見ただけではなかなかうまくいかず、結構試行錯誤した(汗)

 次の簡単な線形連立方程式をIntel MKL LAPACK95のdriver routineである"gesv"を使って解く:

[2 1][x_1] = [4]
[3 1][x_2] = [2]
⇔ Ax = B

(※行列による表現です。わざわざGoogle Mathematical Formulasを使うのは面倒なので。)

 次のコードは、この連立方程式をIntel MKL LAPACK95のdriver routineである"gesv"を使って解くためのコード。
program sample
use mkl95_lapack
use mkl95_precision
implicit none
real, dimension(2,2) :: A
real, dimension(2,1) :: B

a(1,:)=(/2.0, 1.0/)
a(2,:)=(/3.0, 1.0/)
b(:,1)=(/4.0, 2.0/)

call gesv( A, B )

print *, 'solution:'
print *, B(1,1)
print *, B(2,1)

end program sample


 なお、driver routineとはLU分解等を行うcomputational routineを組み合わせたもので、直接解を出力できる。driver routineを使わない場合、いくつかのcomputational routineを自前で接続しないといけない。


(注意点)
・プロジェクトの「ソースファイル」フォルダを右クリックし、「追加」->「既存の項目」を選択。"...\Intel\Compiler\11.1\048\mkl\include"から"lapack.f90"を選び、追加する。

・ソースコードの"implicit none"の前に、"use mkl95_lapack"と"mkl95_precision"(もしくは、"use lapack95"と"use f95_precision")を書き加える。

・「プロジェクト」タブの「プロパティ」を選択。「リンカー」->「入力」と進み、「追加の依存ファイル」に、アーキテクチャがia32の場合、"mkl_lapack95.lib"を書き加える(アーキテクチャがia64の場合、"mkl_lapack95_ilp64"と"mkl_lapack95_lp64"。ただ、こちらの場合は試していないのでこれは予想)。"mkl_lapack95.lib"にIntelがLAPACKを修正したものが格納されていて、企業秘密ということになっている。コンパイラー利用者はコンパイル時にそのブラックボックスにリンクすることでIntel MKL版のLAPACKを利用する。なお、"mkl_lapack95.lib"は"...\Intel\Compiler\11.1\048\mkl\ia32\lib"にある。

・同じく「プロジェクト」->「プロパティ」->「Fortran」->「ライブラリー」->「インテル(R)マス・カーネル・ライブラリーを使用する」の部分を、「並列」か「シーケンシャル」を選ぶ(「クラスター」はおそらくPCクラスター向け)。

連立方程式を解く

 連立方程式を数値計算で解く方法にはいくつかある。写像が(連続)微分可能か否かで分類することができる。

 ○微分可能な場合
 ・Newton-Raphson method

 ○微分可能でない場合(※もちろん、微分可能な場合にも適用できる)
 ・Bisection search
 ・Secant method

いずれも一変数である場合は利用は単純。ただ、変数が複数ある場合、Bisection searchはよくわからん。Secant methodは複数あるときでも何らかの仮定の下では比較的簡単に実践できると思われる。Newton-Raphson methodについてとにかく偏微分してしまえばいいので、特に複雑なことをするわけでもない。

ベイズ統計学のメモ

 ほとんど勉強したことないから、ベイズ統計学の考え方がよくわかってない。覚書をいくつか。

(1)ベイズ統計学では母数(パラメータ)が定数である必要はない。確率変数でもOK。

(2)事前確率密度と尤度から事後確率密度を計算する。解釈としては、パラメータの分布をアップデートする。事前・事後確率密度の両者はパラメータの分布についての分析者の主観的な予想と解釈することができる。

(3)小標本においても柔軟に統計的推論ができるという点でベイズ統計学にはメリットがある。例えば、パラメータについて、どんな値をとりそうか全くわからないし、偏見も何もないという状況にあるとしたら、十分大きな閉区間上の一様分布を事前確率密度として与えるとか。

(4)重要な論点として、標本が大きくなるにつれて、ベイズ的な推論の下で真の分布に収束するのかということがある。ちゃんと真の分布に収束するなら、ベイズ統計学は古典的統計学を一般化させたものと考えることができる。一方で、収束しないなら、思想的選択を迫られる?パラメータが定数だと考えるなら古典的な統計学を使えばよい。じゃあ、その場合のベイズ的な推論って何なの?この辺が理解できていない。

Thursday, April 8, 2010

Rをインストール

 講義に使うかもしれないので、統計解析向けのプログラミング言語のRをインストール。versionは2.10.1。windows版だからexeファイルをクリックするだけでらくちんインストール。いちおう、Matlabと同じような使い方ができるフリーソフト。インタプリタ言語で遅いのでは?と思われるが、C/C++やFortranと連携することができるらしい。詳しい方法はわからないが。

Wednesday, April 7, 2010

Fortran: 変数のデフォルト値

 Fortranでは、実数(スカラーで、配列ではないもの)のデフォルト値はわけのわからん数値になってしまうので、注意する必要がある。
 その他の実数の配列とスカラーとしての整数、整数の配列についてはデフォルト値は0になっている。

Tuesday, April 6, 2010

Ubuntu搭載ノートPC

 だいぶ前からDellがLinux OS搭載のノートPCを販売していた。現時点では、Inspiron Mini 10vというバージョンで、これにはUbuntu 8.04がインストールされているらしい。とっても安いし、UNIXの練習用としてこれを買ってみようかな。

Thursday, April 1, 2010

MPIを使った並列化の実践(その1)

 某所でMPI(Message Passing Interface)の講習を受けてきた。話を聞いていたら、自分の目的を達成するには、ほんとに基本的な構文だけ知っていればOKのような気がした。

 少なくとも今の俺には、次の問題が解ければいい:


ここで、gは非負の関数とする。さらに、Θはコンパクト。具体例としては、非線形最小二乗法。

 この場合、最上層において、つまりθについてMPIを使って仕事を振り分けることが望ましい(もちろん、g(θ)の計算時間が莫大であることを前提とする)。完全に並列化可能で、ほぼ分散処理となる。
 具体的なコードを次のように書いてみた:

module INCMPI
implicit none
include 'mpif.h'
integer istatus(mpi_status_size)
integer nprocs, myrank, ista, iend, ierror
end
-----------------------------------------------------------------
program MPI
use INCMPI
implicit none
integer, parameter :: N=900
real, dimension(N) :: Theta
real, dimension(:), allocatable :: Objv
real, parameter :: theta_l=0.5, theta_u=1.0
real c, d, obj, int_r
integer i, arg, dummy, int_i

! the space of parameter
Theta(1)=theta_l
Theta(N)=theta_u
c=(theta_u-theta_l)/(N-1)
do i=2,N-1
    Theta(i)=Theta(i-1)+c
enddo

! multi-processing through MPI
call mpi_init(ierror)
call mpi_comm_size(mpi_comm_world,nprocs,ierror)
call mpi_comm_rank(mpi_comm_world,myrank,ierror)

allocate(Objv(nprocs))
int_r=N/nprocs
int_i=floor(int_r)
if (myrank+1<nprocs) then
    ista=int_i*myrank+1
    iend=int_i*(myrank+1)
else
    ista=int_i*myrank+1
    iend=int_i*(myrank+1)+mod(N,nprocs)
endif

! minimization
obj=99999
do i=ista, iend
    call f(Theta(i),d)
    if (d<=min) then
        obj=d
        arg=i
    endif
enddo
Objv(myrank)=obj

! communication
call mpi_allgather(mpi_in_place, dummy, dummy, Objv, 1, mpi_real, mpi_comm_world, ierror)

! summarizing the result
if (obj==min(Objv)) then
    print *, 'The optimal parameter is', Theta(arg),'.'
    print *, 'The minimized value of the objective function is', obj,'.'
endif

deallocate(Objv)
call mpi_finalize(ierror)

end program MPI

-----------------------------------------------------------------
subroutine f(theta,d)
real theta, d
d = g(theta)
end subroutine f
-----------------------------------------------------------------
real, function g(theta)
real theta
g = ... ! nonnegative real number
end function g


ただし、このコードはまだデバッグしてない。まだMPIを使える環境にないため、こんなコードでよかろうに、という想像の世界でございます。

なお、''mpi_in_place''を使っているので、MPI-2に対応している必要がある。