任意の確率密度分布に従う乱数の生成（von Neumannの棄却法）

投稿日時: 2013年10月13日投稿者: pacocat

任意の分布に関して乱数を作りたかったのでメモ。
分布の逆関数を容易に求められる場合は、逆関数法を使ってしまうのが簡単なのだけど、

どんな（一価の）関数にも適用出来る
アルゴリズムが作りやすい

といった理由でここでは棄却法を使います。

以下のページを参考にしています。
数値計算法 2011/6/29 (pdf)
15 乱数とモンテカルロ法 (pdf)

von Neumannの棄却法

考え方はとてもシンプルで、

分布関数を含む範囲内で一様乱数のペア ( $x_{i}, y_{i}$ )を生成

具体的には分布関数の範囲を $[x_{min}, x_{max}], [0, y_{max}]$ 、
確率変数 $x_{rand}, y_{rand}$ を $0 \leq x_{rand}, y_{rand} \leq 1$ に従う一様分布と仮定した時に、

$x_{i} = x_{min} + (x_{max} - x_{min})x_{rand}\\ y_{i} = y_{max} y_{rand}\\$

を作る。

生成した乱数ペア( $x_{i}, y_{i}$ ) が分布関数内に入っていれば、 $x_{i}$ を分布に従う乱数として採用（下図でいうと赤点が採用、×が不採用）

アルゴリズム的にはこんな感じ？（コード中にあるテストファイルはこちらから→”random_dist.txt“）

結果のヒストグラム：

正規分布（normal distribution）と指数関数（exponential distribution）の畳み込み積分

投稿日時: 2013年10月12日投稿者: pacocat

手計算でチェックしたかったので過程をメモ。

以下の二つの関数を畳み込み積分した場合、

$f(x) = \begin{cases} \lambda \exp(-\lambda x) & (x\ge 0) \\ 0 & (x<0) \end{cases}$
$\displaystyle g(x) = \frac{1}{\sqrt{2\pi}\sigma}\exp\left\{-\frac{(x-\mu)^2}{2\sigma^2}\right\}$

Exponentially modified Gaussian distribution
にあるようにex-Gaussian distributionと呼ばれる分布が出来る。
$\displaystyle f(x)\ast g(x) = \frac{\lambda}{2}e^{\lambda \mu + \frac{1}{2}\lambda^2\sigma^2}e^{-\lambda x}\mbox{erfc}\left(\frac{\mu+\lambda \sigma^2 -x}{\sqrt{2}\sigma}\right)$
ここでerfcはガウスの相補誤差関数（complementary error function）で以下のように定義されている。erfはガウスの誤差関数。
$\displaystyle \mbox{erfc}(x) = 1-\mbox{erf}(x) = \frac{2}{\sqrt{\pi}} \int ^{+\infty}_{x} e^{-t^2}dt$

Convolution of the normal and exponential probability density functions

以下計算、
$\displaystyle f(x)\ast g(x) = \int^{+\infty}_{-\infty} dx'f(x')g(x-x')\\[1.0ex] =\int^{+\infty}_{0} dx' \lambda \exp(-\lambda x') \frac{1}{\sqrt{2\pi}\sigma}\exp\left\{-\frac{(x-x'-\mu)^2}{2\sigma^2}\right\}\\[1.0ex] = \frac{\lambda}{\sqrt{2\pi}\sigma} \int^{+\infty}_{0} dx' \exp\left[-\frac{1}{2\sigma^2} \left\{x'^2+2x'(\mu-x+\lambda \sigma^2)+(x-\mu)^2 \right\} \right]\\[1.0ex] = \frac{\lambda}{\sqrt{2\pi}\sigma} \int^{+\infty}_{0} dx'\exp\left\{-\frac{1}{2\sigma^2}(x'-x+\mu+\lambda \sigma^2)^2 + \frac{1}{2}\lambda^2 \sigma^2 -\lambda x + \lambda \mu \right\}\\[1.0ex] = \frac{\lambda}{\sqrt{2\pi}\sigma} e^{\lambda \mu + \frac{1}{2}\lambda^2\sigma^2}e^{-\lambda x}\int^{+\infty}_{0} dx'\exp\left\{-\frac{1}{2\sigma^2}(x'-x+\mu+\lambda \sigma^2)^2 \right\}$

ここで、変数変換
$\displaystyle z = \frac{1}{\sqrt{2}\sigma}(x'-x+\mu+\lambda \sigma^2), \; dz' = \frac{dx'}{\sqrt{2}\sigma}$
とガウスの相補誤差関数の定義を使えば、
$\displaystyle f(x)\ast g(x) = \frac{\lambda}{\sqrt{\pi}} e^{\lambda \mu + \frac{1}{2}\lambda^2\sigma^2}e^{-\lambda x}\int^{+\infty}_{\frac{\mu+\lambda \sigma^2-x}{\sqrt{2}\sigma}} dz e^{-z^2}\\[1.0ex] = \frac{\lambda}{2}e^{\lambda \mu + \frac{1}{2}\lambda^2\sigma^2}e^{-\lambda x}\mbox{erfc}\left(\frac{\mu+\lambda \sigma^2 -x}{\sqrt{2}\sigma}\right)\\$
が得られる。

Python code

Pythonに用意されてるモジュールについてはnumpy.convolveも参照

#!/opt/local/bin/python import numpy as np import scipy.special # ex-Gaussian distribution def convolve_exp_norm(alpha, mu, sigma, x): co = alpha/2.0 * np.exp( alpha*mu+ alpha*alpha*sigma*sigma/2.0) x_erf = (mu + alpha*sigma*sigma - x)/(np.sqrt(2.0)*sigma) y = co * np.exp(-alpha*x) * (1.0 - scipy.special.erf(x_erf)) return y # input parameters x = np.arange(-5.0,5.0,0.1) alpha = 1.0 # index of the exponential function sigma = 0.5 # dispersion of the normal distribution mu = -0.5 # center of the normal distribution exponential = [] for i in range(len(x)): value = alpha * np.exp(-alpha*x[i]) if (x[i]>=0.0) else 0.0 exponential.append(value) np.array(exponential) normal = 1.0/(sqrt(2.0*np.pi)*sigma) * np.exp(-(x-mu)*(x-mu)/(2.0*sigma*sigma)) # convolve convolution = convolve_exp_norm(alpha, mu, sigma, x) # plot result plot(x,convolution,color='red',linewidth=2.0) plot(x,exponential,color='blue') plot(x,normal,color='green') legend(('Exponentially modified Gaussian','Exponential Distribution (input)','Normal Distribution (input)'),frameon=False) xlim(-2,5) text(3.0,0.7,r'$\lambda=1.5$',fontsize=20,verticalalignment='center') text(3.0,0.6,r'$\mu=-0.5$',fontsize=20,verticalalignment='center') text(3.0,0.5,r'$\sigma=0.5$',fontsize=20,verticalalignment='center')

結果はこんな感じ。

4月25日

投稿日時: 2013年4月26日投稿者: pacocat

[スペイン語] Colombia extradita a EE.UU. a poderoso ‘narco’

Spanish	English (Japanese)
delincuencial	criminal
fritanga	cheap, greasy
extraditar	extradite（引き渡す）
tonelada	ton
ilicita	ilicit（不法の）
testigo	witness
operativo	operative（刑事）
delinquir	to commit a crime
adicionalmente	additionally
alianza	alliance

マルコフ連鎖モンテカルロ法入門-1 by @teramonagi

マルコフ連鎖モンテカルロ法入門-2 by @teramonagi

未来を予想する最善の方法はそれを発明することではない:動画で学ぶGUIの歴史と未来

MBAとはどういう教育なのか – 愛の日記
これはホントにいいエントリーだった。

4月24日

投稿日時: 2013年4月24日投稿者: pacocat

[スペイン語] Canadá: Supuestos terroristas dan la cara en corte
カナダでテロの容疑者二人逮捕。国籍は明らかにしてないが、関係者の話では一人はチュニジアの出身。先日のアメリカといい、この種の事件では国とか国籍とかもはや何でもありになってきてますよね。

Spanish	English (Japanese)
capturó	to capture, to seize
sujeto	subject, suspicious indivisdual
sospechar	to suspect
vínculo	link, tie, bond
relacionado	related, connected
sacudir	to shake
acusar	to accuse, to charge
culpable	guilty
juez	judge
divulgar	to spread, to disclose
prueba	test, proof, trial
corte	court
ciudadano	citizen
tunecino	Tunisian

島森路子さんについて – CHIAKI KASAHARA BLOG
“「彼女のインタビューの基本は、「話しているその人のパーソナリティと、話している内容を同時に浮かび上がらせる」というものだろう。「人が話す」ということは、「そこに明確なパーソナリティが存在していることだ」という、「考えてみれば当たり前のことを、彼女はいち早く理解していた。だから彼女は「普通の言葉」や「話す言葉」を重要視していた。」（「島森路子インタビュー集」の橋本治氏の寄稿より）”

「盗んだダイコン並べだす」　大根100本盗んだ中3女子を逮捕　路上に並べ車がダイコンをひいていく音と形の変化などを見て楽しむ

IPSJ Magazine 巻頭コラム – 情報と哲学（東浩紀）

哲学は、知覚可能な集合体と、それら世界の背後にあるものを区別する表層-深層対立図式でできている
カントの「現象-物自体」、ハイデガーの「存在者-存在」、深層を否定したのはニーチェ
1930年代からの哲学と科学の分裂（カルナップのハイデガー批判）
科学から自由に世界の背後を追求し続けたのが「大陸哲学」（オカルト）、認識形式や言語形式を地味に分析したのが「分析哲学」（オタク）
情報は抽象的かつ具体的な存在であり、哲学的対象として扱うには構造的問題が存在する

シリコンバレーによろしく: 米Yahooが17歳の少年に社運を賭ける、新体制マリッサ・メイヤーの一手

米Yahooはモバイルへのシフトに失敗しイケてない企業に。2012年7月にマリッサ・メイヤーがCEOに就任。
企業が新しいことを生み出すためには「小さくいる」ことが重要。大きな企業は保守化する。

【SYNODOS】なぜ彼女たちは出会い系に引き寄せられるのか／「彼女たちの売春（ワリキリ）」著者、荻上チキ氏インタビュー

ねとぽよ「朝まで生ソシャゲFINAL」観客の反応

自主規制でトレード禁止になってから無課金ユーザーの逆転が難しくなった（e.g.ドラコレ）
カードゲームの飽き？
北米では時短などデメリットを解消する後ろ向き課金から欲しいからアイテムを買う前向き課金へ
パズドラはゲーム＋ソーシャルであってソーシャルゲームではない？（ソーシャル性が隠されてる）
パズドラの勝因は、対人勝負のソーシャル疲れに誠実に対応したこと
ガラケーの5押しの快感をスマホのなぞる操作に完全に置き換えた
北米ではパズドラにPvPを追加して欲しいとの声も
コンシューマーゲームはコンテンツを楽しむが、ソシャゲはコンテンツを楽しむものではないので課金ポイントが変わる
海外ではゲームをスクロールしない
面白ゲームの探索コストを減らすためにSNSベースのソーシャルグラフを利用したソシャゲが出てきた？

4月22日

投稿日時: 2013年4月24日投稿者: pacocat

[スペイン語] Turquía desoye pedido de EE.UU. sobre Gaza
トルコのパレスチナ訪問が延期になった話。あの辺りはいつも緊張してるよなぁ。

Spanish	English (Japanese)
apropiado	appropriate
hacerse eco de algo	to report something
desde el punto de vista	from the point of view
rueda de prensa	Press Conference
solicitar	to request
citar	to quote, to mention, to summons
juicio	trial
de acuerdo con	in accordance with
acceder	to accept, to access
aplazar	to postpone, to delay

GiGiインタビュー　36歳無職・ゲーム/IT業界未経験者が、人気タイトルの企画屋に転身するまで – ねとぽよ

LINE POPで500万点だと…！？　その”奥義”を廃プレイヤーに聞いてみた – ねとぽよ

新経済サミット2013：「エンジニア1人で世界を変えられる時代、日本からイノベーションを起こすには」　伊藤穰一氏、Rubyまつもと氏など議論 (1/3) – ITmedia ニュース

4月10日

投稿日時: 2013年4月13日投稿者: pacocat

國分功一郎が語る「浪費としてのファッションはありえるのか？」- STUDIO VOICE
“ここで最初の問いに戻ると、最初にぼくが言った意味でファッションが消費的である可能性はもちろんあります。流行っているからという理由で買うのは消費的なファッションです。でも、そもそもファッションは本質的に装飾であって実質や機能性というものはそのあとでくっつけたられたものだとしたら、ファッションというのは本質的に消費というわけではないということになりますよね。”

統計解析ツール「R」、8年半ぶりのメジャーバージョンアップ版「R 3.0.0」リリース | SourceForge.JP Magazine

「経済入門　～金融緩和で円安の仕組み～」：Mプラス：テレビ東京

広まっている「インド人の宇宙観」は誤解がある！

世界で初めてアマゾン違法伐採の監視ツールを開発した72歳技術者が、今も土日にコードを書く理由 – エンジニアtype

4月9日

投稿日時: 2013年4月9日投稿者: pacocat

[スペイン語] Obama: “No olvidamos la promesa que hicimos”
米憲法修正第2条を守ったままどうやって火器の使用を押さえるかという問題。

Spanish	English (Japanese)
cobrar	charge
mandatario	representative, attorney, mandatory
aprobación	approval
atentar	to assault, threaten
asunto	issue
sentido común	common sense
a la vez	at the same time
enmienda	amendment（法改正）
receso	break
medida	measurement, dimensions
propuesta	proposal
antecedente	antecedent（先立つ、前提、仮定的な）
feria	fair
estricta	strict
propósito	purpose
vecindarios	neighborhood

[フランス語] Margaret Thatcher est morte
マーガレット・サッチャー元英首相死去。

French	English (Japanese)
succomber	succumb（屈する、死ぬ）, die
accident vasculaire cérébral
à cause de	because of
fermeté	firmness
paisiblement	peacefully
porte-parole	spokesperson

はてなに入った技術者の皆さんへ – jkondoのはてなブログ
アウトプットを出す事、疑う事、続ける事

4月8日

投稿日時: 2013年4月8日投稿者: pacocat

小浜温泉でバイナリー発電　実証試験スタート　長崎 – MSN産経ニュース
“豊富な源泉を誇る小浜温泉は、湧出した湯の７割を使用せずに廃棄していた。この温泉水（９５度）を使い、沸点が１５度と低いフッ素化合物を気化させ、発電に利用する。”
出力は150kW。

17歳の高校生が弟の赤毛の要因を解明すべくビデオデッキを活用して“DNA増幅装置”を自作

ゲーム内広告の新たな試み。シムシティに日産「リーフ」充電スポットが登場

パスティーシュ – wikipedia

データの鬼、Googleが解析した「よい上司を製造する8つの条件」 – DNA

２７分発の電車が２６分に発車したのですが……

7年以内に中国が世界最大の高級車市場に、日米独は中国人の好みに合わせた開発を―米誌

今後7年以内に中国が米国を追い抜いて世界最大の高級車市場になる。
米国では装備の整ったフォード・エクスプローラーを持っていれば、高級車を持っていると見なされる。だが中国では高級車は基本的にドイツ車と同義語だ。多くの中国人にとって高級車とはベンツ、BMW、アウディのみを指す。
中国の一人っ子世代はすでに結婚し、仕事を持つ年齢に入っている。彼らはより多くのチャンスを享受するだけでなく、親から多額の財産を継承する可能性がある。この世代の若者はまだ一家を養う圧力がないため、消費力が非常に強い。

全米のラジオ局がKE$HA「Die Young」の放送を自粛

インターネットはMtGを変えた　～Magic the Gathering 20周年に寄せて～ – ねとぽよ
メタゲームの活性化とコミュニティのフラット化

4月4日

投稿日時: 2013年4月4日投稿者: pacocat

[スペイン語] Contemplan “tarjeta azul de residencia” para agricultores
移民政策。農業従事者のビザが取りやすくなったという話。調べてみたらH-2Aビザは短期季節的労働者用なんですね。申請基準が緩和されたらしい。

Spanish	English (Japanese)
migratoria	migratory
reclamo	claim, protest
empresarial	business
propuesta	proposal
jornalero	laborer
permanecer	to remain
acceder	to accept, to access
promedio	average
estipular	to stipulate（取り決める、要求する）
diseñar	to draw, to design, to outline
culminar	to conclude, to finish
cosecha	harvest
generar	to generate, to create

[フランス語] Remises en liberté massives de sans-papiers suite aux réductions budgétaires
上と似たような記事。アメリカの移民政策について。

French	English (Japanese)
sans-papiers	undocumented immigrant（正式書類のない移民）
budgétaire	budgetary（予算上の）
initié	initial
disposer	to dispose（配置する、処理する）
en détention	under detention（拘留中）, in costody（保護観察下）
relâche	rest
renvoyé	return
déporter	to deport（追放する）

『LINE』がスペインで会員数1000万人を突破…有名俳優を起用したプロモーションが奏功
“現在、世界1億3000万人・日本国内4500万人以上のユーザーが利用している。今回、ヨーロッパ圏で初めて、スペインの国内ユーザー数が1000万人の大台を突破した。”

芸術係数blog 芸術係数読書会：アートをめぐる1980年代の言葉を読む／ボードリヤール、ブルデュー、オーウェンス、クラウス、フォスター、クルーガー、ピーター・ハリー、ジェイムソン、クーンズ、クリステヴァetc.

「ミッション・インポッシブル」に挑む黒田日銀総裁 – WSJ.com
“エコノミストの中には、デフレが余りに長年続いたため、日銀は真剣であると国民に確信させるためにも過度に大胆な目標を設定しなければならないのだと述べる。”

【WEBマーケター必見】ABテスト最新事例19選

4月3日

投稿日時: 2013年4月4日投稿者: pacocat

[スペイン語] EE.UU. pide se respete libertad de prensa
ホンジュラスの情報通信法（？）改正について。

Spanish	English (Japanese)
actuación	performance
adecuado	appropriate
equilibrio	balance
proteger	to protect
a través de	through, across
revelar	to reveal, to disclose, to show
alentar	to encourage
propietario	owner

[フランス語] Le virus H7N9 de la grippe aviaire fait de nouvelles victimes
中国のH7N9型インフルエンザについて。病気にかかるのは英語でcontractを使うんですね。

French	English (Japanese)
aviaire	avian（鳥類の）
au moins	at least
contracté	contracted
volailles	poultry（家禽）

ドコモのキャリアフリーなO2Oサービス『ショッぷらっと』を生んだ“火星人チーム”とは？（エンジニアtype）
メモ：「iモードが全盛の時代には、NTTドコモが自らビジネスモデルを構築できました。しかし今日、スマートフォン向けサービスの主導権を持っているのは、通信キャリアではなくアップルやグーグルのような企業です。そこでNTTドコモとして、スマートフォンにおいてどのようなサービスを提供すべきか慎重に検討した結果、導き出したのがO2Oサービスだったんです。」

The Economistによる、働く女性として避けたい国で日本はワースト2位。

３Dプリンターの反政府運動と「オープン・ソース」
テクノロジー悲観主義の立場として。確かに”openness”の意味が状況によって変わるのは問題だろうなぁ。

Pacocat's Life

Every day is a new day

作成者アーカイブ: pacocat