たつをの ChangeLog : 2014-04-16

2014年4月16日のヲハニュースをお届けします。



▼「スイッチ押したら朝が来る」アニメの世界が現実に。『Sony α7S』の超高感度が、夜を消し去る驚愕レベル | DDN JAPAN
おお！


▼店頭でむき出しになっている食品は不衛生
買うのに躊躇するよなあ……。
（関連：コンビニのおでんは食べません[2007-11-13-2]）

▼電磁波が怖い。wifiが怖い。電波恐怖症の人が次々の避難するアメリカの小さな村 : カラパイア
Wikipedia:電磁波過敏症

▼2013年電子書籍の収支 - CHINGE
ふむふむ。

▼育休中には自宅では仕事はできないと言われているけど、本当のところはどうなの？ | 菅原康之
今日のもやもや話。在宅勤務で育休だといろいろと。「育休中でも仕事できるだろ？ やれよ」みたいな話につながりそう。

▼BiC4色やフリクションボール3に装着できるタッチペン「SMART-TIP」と筆箱まとめ - goodegg.jp
BICの4色ボールペンを使っているので気になるなあ(ref. [2013-11-13-1])。それゆえタッチペン買うとボールペンついてくる、ってのは個人的に余計。
Amazon.co.jp： SMART-TIP タッチペン （Bic 4色ボールペンセット） シルバー for iPhone iPad tablet PC product by UNUS PRODUCT SERVICE.: パソコン・周辺機器


▼斜めに書く人のための「まっすぐノート」が斬新かつ実用的！(Excite Bit コネタ) - エキサイトニュース(1/2)
おもしろーい。
Amazon.co.jp： 斜めに書く人のための【まっすぐノート】 SL3070: 文房具・オフィス用品


▼【SF旅行2014】 スマートフォンは現代のタバコである | IDEA*IDEA
なるほどー！

▼「STAP細胞の有無」をアンケートで調べて、何がしたいのか - 俺のメモ
「信仰」というよりも、「宇宙人はいる？いない？」みたいなネタレベルになってしまったということかも。

▼カレーのココイチが、タイで人気沸騰のワケ （東洋経済オンライン） - Yahoo!ニュース BUSINESS
高級路線。なるほど。

メモメモ。この前の週刊アスキー[2014-03-25-2]の特集「Amazonだけで新生活」より、レンジ調理器具いろいろ。


レンジでゆで卵を作る道具。10分くらいでできる。3個用、4個用もある。会社でゆで卵作れるなあと思うも、10分も電子レンジを占拠するのは迷惑だなあ。

■レンジでらくチン ゆでたまご2ケ用 RE-277



レンジでラーメン。金属製でもなく陶器製でもなくプラスチック製だからか安い。でも普通のラーメンどんぶりにフタすれば良いかも？

■パール金属 電子レンジで作る新潟産のラーメンどんぶり C-249



レンジで炊飯。いざというとき便利か。うちは電気炊飯器買ったばかりだからいらないけど。

■電子レンジ専用炊飯器 備長炭入り ちびくろちゃん 2合炊き

取りうる値の個数が有限個の任意の離散分布に従う乱数を発生させる「別名法 (alias method)」を Perl で実装してみました。ロジックは下記参考文献に載っていたのそのままで、ソース中のコメントは引用となっています。

■東京大学教養学部統計学教室 (編集), "自然科学の統計学", 東京大学出版会, 1992.


別名法は、例えば「大吉15%、中吉30%、吉30%、凶20%、大凶5%」の割合でランダムにおみくじを出すプログラムを書くときの効率的なアルゴリズムです。

深く考えない実装だと、

$r = rand(1)
if ($r < 0.15) {
  return "大吉";
} elsif ($r < 0.45) {
  return "中吉";
} elsif ($r < 0.75) {
  return "吉";
} elsif ($r < 0.95) {
  return "凶";
} else {
  return "大凶";
}

みたいになるかと思います。このおみくじは5種類だからそれほどではないけど、これが100種類とか1000000種類とかだとどうでしょう。種類数に比例した数の比較演算が必要で時間がかかります。O(N)ですね。（ちょっと工夫して二分探索を使うと O(logN) になります。）

ところが、別名法では準備にほんの少し時間がかかりますが、それ以降は種類数によらずたった1回の比較演算で済みます。O(1)です。すごいですねえ。

どういう仕組みなのかというと……。

って、余裕があるときに「図を用いた分かりやすい説明」を書こうと思います。のんびりお待ちください。

以下、コードと実行例です。

■動作確認スクリプト (aliasmethod.pl)：

#!/usr/bin/perl
use strict;
use warnings;

###### 別名法 (alias method)
### 取りうる値の個数が有限個の任意の離散分布
###   P(X=x_{k}) = p_{k} (1 <= k <= K)
### に従う乱数を発生させるためにウォーカー(Walker, A.J.)が考案した
### 巧妙な方法である。逆関数法と違って、K がいくら大きくても
### 1個の乱数の発生に必要な比較演算が1回だけであるという特徴を持つ。
### ただし、そのために最初に次のように表を用意しておく必要がある。
### (参考文献："自然科学の統計学", 東京大学出版会, 1992)

my @ns = @ARGV;
my $K = @ns; # K
my $sum = 0;
foreach my $i (@ns) {
    $sum += $i;
}
my @p = map {$_ / $sum} @ns; # p_{k}
print "p_{k} :@p\n\n";

### [初期設定] 次の手順に従って a_{k}, v_{k} (1 <= k <= K) を求める

### 1 : v_{k} ← K*p_{k} (1 <= k <= K)
my @v = map {$K * $_} @p; # v_{k}

### 2 : v_{k} >= 1 を満たす k の集合を G,
###     v_{k} < 1 を満たす k の集合を S　とする。
my @G;
my @S;
foreach my $k (0..$K-1) {
    push @G, $k if $v[$k] >= 1;
    push @S, $k if $v[$k] < 1;
}

### 3 : S が空集合でない限り、3.1 から 3.4 までを繰り返し実行する。
my @a; # a_{k}
while (@S and @G) {
    ### 3.1 : S, G から1つずつ任意の要素を選ぶ。それを k, l とする。
    my $k = $S[0];
    my $l = $G[0];

    print "v_{k}: @v\n";
    print "S: @S\n";
    print "G: @G\n";
    print "k=$k, l=$l\n";
    print "\n";

    ### 3.2 : a_{k} ← x_{l}, v_{l} ← v_{l} - (1 - v_{k})
    $a[$k] = $l;
    $v[$l] = $v[$l] - (1 - $v[$k]);

    ### 3.3 : S から k を除く。
    shift @S;

    ### 3.4 : v_{l} < 1 なら l を G から S に移す。
    if ($v[$l] < 1) {
        push @S, (shift @G);
    }
}

print "k p_{k} Kp_{k} v_{k} a_{k}\n";
for (my $i = 0; $i < $K; $i++) {
    print "$i $p[$i] ".($K*$p[$i])." $v[$i] ".(defined $a[$i]?$a[$i]:"*")."\n";
}
print "\n";

### [乱数発生]
sub gen {
    my ($a_r, $v_r) = @_;
    ### 1 : (0, K) 上の一様乱数 V = KU を発生する。
    my $V = rand(@$v_r);
    ### 2 : k ← L V 」+ 1, u ← k - V
    my $k = int($V);
    my $u = $V - $k;
    ### 3 : u <= v_{k} ならば X ← x_{k}, そうでなければ X ← a_{k}
    if ($u <= $v_r->[$k]) {
        return $k;
    } else {
        return $a_r->[$k];
    }
}

### 実験
my $N = 1000000;
my @count;
for (my $i = 0; $i < $N; $i++) {
    my $X = gen(\@a, \@v);
    $count[$X]++;
}
for (my $i = 0; $i < $K; $i++) {
    print "$i ".("*"x($count[$i]/10000))." $count[$i]\n";
}

■実行例：

% aliasmethod.pl 0.12 0.24 0.38 0.16 0.10
p_{k} :0.12 0.24 0.38 0.16 0.1

v_{k}: 0.6 1.2 1.9 0.8 0.5
S: 0 3 4
G: 1 2
k=0, l=1

v_{k}: 0.6 0.8 1.9 0.8 0.5
S: 3 4 1
G: 2
k=3, l=2

v_{k}: 0.6 0.8 1.7 0.8 0.5
S: 4 1
G: 2
k=4, l=2

v_{k}: 0.6 0.8 1.2 0.8 0.5
S: 1
G: 2
k=1, l=2

k p_{k} Kp_{k} v_{k} a_{k}
0 0.12 0.6 0.6 1
1 0.24 1.2 0.8 2
2 0.38 1.9 1 *
3 0.16 0.8 0.8 2
4 0.1 0.5 0.5 2

0 *********** 119632
1 ************************ 240086
2 ************************************* 379950
3 *************** 159991
4 ********** 100341

■短いソース：お作法的にはアレ（特にmap）だけど短い版。

#!/usr/bin/perl
use strict;
use warnings;

my @ns = @ARGV;
my ($a_r, $v_r) = pre(@ns+0, \@ns);

my $N = 1000000;
my @c = do {my @t; map {$t[gen($a_r, $v_r)]++} 1..$N; @t};
print join("", map {"$_ ".("*"x($c[$_]/($N/100)))." $c[$_]\n"} 0..$#ns);

sub pre {
    my ($K, $ns_r) = @_;
    my $sum = do {my $s; map {$s += $_} @$ns_r; $s};
    my (@v, @a, @p, @G, @S);
    map {$_ /= $sum; push @p, $_; push @v, $K * $_} @$ns_r;
    map {$v[$_] >= 1 ? push @G, $_ : push @S, $_} 0..$K-1;
    for (my ($k, $l); $k = shift @S, $l = $a[$k] = $G[0]; ) {
        $v[$l] = $v[$l] - (1 - $v[$k]);
        push @S, shift @G if $v[$l] < 1;
    }
    return \@a, \@v;
}

sub gen {
    my ($a_r, $v_r) = @_;
    my $V = rand(@$v_r);
    my $k = int($V);
    return $V - $k <= $v_r->[$k] ? $k : $a_r->[$k];
}

関連

- unnonouno: 高速な復元抽出の直感的な説明

- discrete_distribution(C++11)

- 非一様乱数 - ORWiki

別名法(alias method)。「取りうる値の個数(n), が有限個の離散分布に従う乱数を, nに無関係なO(1)の速度で生成できる, 簡単で効率的な方法である. ただし, 初期設定にO(n)の時間とO(n)のメモリを必要とする.」

- Alias method - Wikipedia, the free encyclopedia

In computing, the alias method is a family of efficient algorithms for sampling from a discrete probability distribution, due to A. J. Walker. The algorithms typically use O(n log n) or O(n) preprocessing time, after which random values can be drawn from the distribution in O(1) time.

追記140605:
- RSS を読み込んでランダムな文章を生成する[2008-09-07-3]
（別名法を知る前の私の非効率的な実装。確率分布によりランダムにざっくり選ぶ。）