球面上の点を平面上にステレオ投影する

先週のさらに続き。先週は、球面上の点を平面上に投影するのに X, Y 座標をそのまま使って水平投影したけど、今度はステレオ投影というのをやってみる。
ステレオ投影は Wikipedia に詳しく載っている。
先週と同じように、下半球の点 300 個を XY 平面(赤道面)上に投影する。やってみたのが下の図。

コードは次の通り。

やっつけ仕事なので、ヘルプもないけど、先週の gss_gen コマンド(--cartesian オプション付き)の結果を入力として、投影点の座標を CSV で出力する。

多数の点を球面上に一様分布させる(2)

一昨日の続き。最初と最後の点を極から移動(位置調整)するのを実装した。ついでに、gem 化して RubyGems.org にアップしておいた。gem install gss_generator でインストールできる。
使い方は次の通り。

^o^ > gss_gen --cartesian --relocate 1.0 600 > c600r.csv

--cartesian はデカルト座標系で出力、--relocate は位置調整のオプション。結果は CSV 形式で標準出力に出力するので c600r.csv ファイルにリダイレクトしている。

結果を一昨日と同様の図にしたものを載せる。

最初の点が極(中心)からずれて、一様性が増しているのがわかる。

多数の点を球面上に一様分布させる

ちょっと面白いものを見つけた。

 cf. LEDドームのLEDの並びを決めるのに使用した計算式 – jakaladaのブログ

球面上に任意個の点を一様分布させる、座標を求めるもの。一般化螺旋集合(generalized spiral set)というものを使っている。元の資料(論文)はここ。

 cf. 多数の点を球面上に一様分布させるソフトウェアGSS Generator

上のリンク先のブログでは Python を使っているので、Ruby でやってみた。

最初のがコマンドで、あとの2つがライブラリ。引数に球の半径と配置したい点の個数を指定すると、各点の極座標を CSV 形式で出力する。

^o^ > ruby -Ilib exe/gss_gen 1.0 600

極座標じゃなく、デカルト座標(XYZ座標)がほしい時には --cartesian オプション。

^o^ > ruby -Ilib exe/gss_gen --cartesian 1.0 600

下の図は、発生させた600個の点のうち下半分の300個を XY 平面上に投影した図。

各点を順に線でつないでみると、螺旋になっているのがよくわかる。

……なんか螺旋が逆になってるな。元論文じゃ球面を下から見たと書いてあるからそのせいか?

ともあれ、それらしいのはできた。ただし、最初と最後の点の座標を調整するのはやってない。時間があったらやってみよう。

文字列間のレーベンシュタイン距離を求める(3)Haskell版ふたたび

去年の12月に2つの文字列のレーベンシュタイン距離を求めるっていうのを、JavaScriptHaskell でやった。もともとは Python での実装を参考に2次元配列を使ってやってみたもので、JavaScript版はともかく Haskell版は2次元配列の更新に苦労した。
それが、今日ふと1次元配列でやってみたらどうだろうを思いついた。
つまりこうだ。(m + 1)行×(n + 1)列(m、n は比較する文字列2つの長さ)の各行をつなげた1次元配列を考えて、2次元配列の時の座標を表すタプル (i, j) で初期化しておく。0 < i, 0 < j のとき、(i, j) の値はひとつ左(i-1, j)、ひとつ上(i, j-1)、左上(i-1, j-1)から決まるから、これを1次元配列のインデックスに直すと次のようになる:

  • ひとつ左: i * (n + 1) + j – 1
  • ひとつ上: (i – 1) * (n + 1) + j
  • 左上: (i – 1) * (n + 1) + j – 1

これをコードに落としこんでやるとこうなった:

結果はこうだ。

takatoh@apostrophe $ runhaskell ld2.hs apple play
4
takatoh@apostrophe $ runhaskell ld2.hs perl pearl
1

OK、うまくいった。

だけど、上のコードはまだ2次元配列の意識を引きずっている。もっと単純にできるはずだ。1次元配列のインデックスを x とすると:

  • ひとつ左: x – 1
  • ひとつ上: x – (n + 1)
  • 左上: x – (n + 1) – 1

となる。これで一般部については2次元配列を気にしなくて良くなった。ただし問題がある。いちばん上の行(第0行)といちばん左の列(第0列)だ。少し考えて、x を (n + 1) で割った商と余りを使えばいいと気がついた。コードにするとこう:

1行だけだけど長くなってしまった。だけど考え方はシンプルのように思う。実行してみよう。

takatoh@apostrophe $ runhaskell ld3.hs apple play
4
takatoh@apostrophe $ runhaskell ld3.hs perl pearl
1

OK。

逆関数法で指数分布する乱数を生成する

[0,1)区間の一様乱数から、指数分布にならう乱数を生成するには、逆関数法というのが使える。
指数分布の密度関数は、パラメータをτとすると:
f(\tau)=\lambda e^{-\lambda\tau}
であり、分布関数 g(τ) は:
g(\tau)=\int^\tau_{-\infty}{\lambda e^{-\lambda\tau}}d\tau=1-e^{-\lambda\tau}
となる。g(τ)は 0~1 の値をとるので、この逆関数:
\tau=-\frac{1}{\lambda}log(1-g(\tau))
の g(τ) の代わりに一様乱数を入力してやれば、τ は指数分布する乱数になる。

じゃあ Ruby でやってみよう。

λ=0.5とし、10000個の乱数を発生させている。
これを Excel でグラフ化したのがこれ。
expon
「指数分布」の曲線は、上に書いた密度関数の曲線を、スケールを合わせるために8000倍して描いている。乱数はちゃんと指数分布になっているようだ。

参考にしたページ:
 cf. http://www.ishikawa-lab.com/montecarlo/4shou.html

Box-Muller法で正規分布する乱数を生成する

一様分布する乱数から、正規分布に従う乱数を生成する方法に、Box-Muller法というのがある。
Wikipediaによれば、(0,1) 区間の一様分布乱数2つ(X,Y)から、下の式で2つの正規分布乱数 Z_1Z_2 が得られる。

Z_1=\sqrt{-2log{X}}\cos{2\pi{Y}}
Z_2=\sqrt{-2log{X}}\sin{2\pi{Y}}

Z_1Z_2 は標準正規分布となるので、これらに標準偏差 σ をかけて平均 μ を足してやれば、任意の正規分布に従う乱数が得られる。

Ruby で 10000個の乱数を発生させるスクリプトを書いてみた。ここでは平均 μ=1.0、標準偏差 σ=0.2 とした。

結果を Excel でグラフ化してみた。水色の点が 0.1 単位のヒストグラム。黄緑の線が Excel に用意されている NORM.DIST 関数で描いたもの(スケールを合わせるために NORM.DIST 関数の値は 1000 倍している)。

Box-Muller

こうしてみると、ほぼぴったりと正規分布になっているようだ。

ちなみに Excel で平均値と標準偏差を求めたら、それぞれ μ=0.997、σ=0.201 となった。

ダイクストラ法(2)

昨日のダイクストラ法を C でやってみた。

takatoh@nightschool $ ./dijkstra
10
0 -> 2 -> 4 -> 5

ダイクストラ法

ダイクストラ法とは、グラフ上の最短経路問題をとくアルゴリズム。↓このページに詳しいアルゴリズムの説明がある。

 cf. ダイクストラ法(最短経路問題) – deq notes

Ruby でやってみた。
例題は、上のページにもあるこのグラフ(ただしノードにつけた番号のせいか上下が逆になっている)。

g

円(ノード)に番号がふられ、ノードをつなぐ辺(エッジ)にはそこを通るときのコスト(距離とも時間とも解釈できる)が付されている。この中の任意の 2 つのノードをつなぐ最短経路とコストを求めるのが最短経路問題だ。
今回は 0 番のノードをスタートし 5 番のノードをゴールとする最短経路とそのコストを求めてみる。

実行結果:

takatoh@nightschool $ ruby dijkstra.rb
10
0 -> 2 -> 4 -> 5

というわけで、最短経路は 0 -> 2 -> 4 -> 5、コストは 10 という答えが得られた。これは上のリンク先の答えと同じなので、あっていると思っていいだろう。

g1

文字列間のレーベンシュタイン距離を求める

レーベンシュタイン距離というのを知った。

Wikipedia のレーベンシュタイン距離の項によると:

レーベンシュタイン距離(レーベンシュタインきょり)は、二つの文字列がどの程度異なっているかを示す距離(距離空間を参照)である編集距離(へんしゅうきょり)の一種類である。具体的には、文字の挿入や削除、置換によって、一つの文字列を別の文字列に変形するのに必要な手順の最小回数として与えられる。名称は、1965年にこれを考案したロシアの学者ウラジミール・レーベンシュタインにちなむ。

だそう。上記ページに擬似コードによる実装例が載ってるのと、次のページの Python での実装例が参考になった。なので JavaScript でやってみた。

 cf. 編集距離 (Levenshtein Distance) – naoyaのはてなダイアリー

実行例:

takatoh@nightschool $ ./ld.js apple play
4
takatoh@nightschool $ ./ld.js perl pearl
1

二分木を使ったソート

単純な二分木を使ったソート。
構造体で二分木を表現するのと、malloc() / free() の練習。

takatoh@nightschool $ ./btreesort
unsorted: 80 77 21 42 50 76 88 48 32 19 
sorted:   19 21 32 42 48 50 76 77 80 88