Table of Contents
置換積分
Problem
$a>0$のとき, 以下を求めよ
\[\int^a_0 x\sqrt{a^2 - x^2}dx\]解答
以下のように変数変換を置く:
\[\begin{align*} x &= a\sin\theta\\ dx &= a\cos\theta d\theta \end{align*}\]このとき, $\theta$は$0$から$\pi/2$を動くので
\[\begin{align*} \int^a_0 x\sqrt{a^2 - x^2}dx &= \int^{\pi/2}_0 a\sin\theta\sqrt{a^2 - a^2\sin^2\theta}a\cos\theta d\theta\\[3pt] &= a^3 \int^{\pi/2}_0\sin\theta\cos^2\theta dx\\[3pt] \end{align*}\]あらためて, $u = \cos\theta$とおくと, $du = -\sin\theta d\theta$なので
\[\begin{align*} a^3 \int^{\pi/2}_0\sin\theta\cos^2\theta dx &= a^3\int^0_1 -u^2 du\\[3pt] &= a^3\int^1_0 u^2 du\\[3pt] &= \frac{a^3}{3} \end{align*}\]自然対数と置換積分
Problem
$x\in(0,1)$で確率密度関数が以下のように与えられたとする
\[f(x) = 2\log(1 + \sqrt{x})\](1) $\int^1_0 f(x)dx = 1$になることを確かめよ
(2) $X$について, 期待値と分散を求めよ
確率事象の値全体としての標本空間を$\Omega$としたとき, 連続型の確率変数の密度関数 $f(x)$は次の性質を満たす必要があります
\[\begin{align*} & (a): f(x) \geq 0\\ & (b): \int_\Omega f(x) dx = 1 \end{align*}\]解答
(1)について以下確かめてみる. $1 + \sqrt{x} = t$とおくと
\(\begin{align*} \int_\Omega f(x) dx &= \int^1_0 2\log(1 + \sqrt{x}) dx \\[3pt] &= 2\int^2_1\log(t) 2(t-1)dt\\[3pt] &= 2\int^2_12(t-1)\log(t)dt\\[3pt] &= 2\bigg[(t-1)^2\log(t)\bigg]^2_1 - 2\int^2_1(t-1)^2\frac{1}{t}dt\\[3pt] &= 2\log2 - 2\int^2_1 t - 2 + \frac{1}{t}dt \\[3pt] &= 2\log2 - \bigg[t^2 - 4t + 2\log(|t|)\bigg]^2_1\\[3pt] &= 1 \end{align*}\)
期待値と分散については計算がめんどくさいですが.
\(\begin{align*} \mathbb E[X] &= \int^1_0 2x\log(1 + \sqrt{x}) dx = \frac{7}{12}\\[3pt] \mathbb E[X^2] &= \int^1_0 2x^2\log(1 + \sqrt{x}) dx = \frac{37}{90}\\[3pt] \text{Var}(X) &= \mathbb E[X^2] - \mathbb E[X]^2 = \frac{17}{240} \end{align*}\)
Python simulation: Acceptance-Rejection Sampling
目標密度関数$f$, 参照密度関数$f_0$, 定数$M$を以下のように定義します
\[\begin{align*} f(x) &= \begin{cases}2\log(1 + \sqrt{x}) & x \in(0, 1)\\ 0 & \text{otherwise}\end{cases}\\[3pt] f_0(x) &= \begin{cases}1 & x \in(0, 1)\\ 0 & \text{otherwise}\end{cases}\\[3pt] M & = \sup\left\{\frac{f(x)}{f_0(x)}\right\} = f(1) = 2\log 2 \end{align*}\]$M = 2\log 2$は目標密度関数が定義域内で単調増加関数であることから自明.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
import numpy as np
from scipy import stats
import plotly.graph_objects as go
def cdf(x):
return 2*(x-1)*np.log(1+np.sqrt(x)) - x + 2*np.sqrt(x)
def pdf(x):
return 2*np.log(1 + np.sqrt(x))
## data generating
np.random.seed(42)
N = 20000
u = np.random.uniform(0, 1, N)
v = np.random.uniform(0, 1, N)
M = 2*np.log(2)
thresh = pdf(v)/M
x = v[u <= thresh]
## plot
support = np.linspace(0, 1, 1000)
density = pdf(support)
binsize= 50
fig = go.Figure()
fig.add_trace(go.Histogram(x=x,
nbinsx=binsize,
name='Rejection sampling',
histnorm='probability'))
fig.add_trace(go.Scatter(x=support,
y=density/binsize,
mode='lines',
name='density function'))
シミュレーションで発生させたデータのhistogramとdensity plotを重ねるとそれなりにうまく行っていそうな気配が漂っている.
KS testで統計量を確認してみると, 異なる分布と棄却することはできない結果となった.
1
2
3
4
5
6
## ks test
stats.kstest(x, cdf)
>> KstestResult(statistic=0.005850481461250778,
pvalue=0.7039166342313163,
statistic_location=0.7378783528059887, statistic_sign=1)
一旦は分布が再現できていると判断して期待値と分散を計算すると以下のような結果となり, 上で求めた厳密解と近しい値を取ってることがわかる.
1
2
3
4
5
## mean, variance
print("sample mean:{:.4f}, true-val::{:.4f}".format(np.mean(x), 7/12))
print("sample mean:{:.4f}, true-val::{:.4f}".format(np.var(x), 17/240))
>> sample mean:0.5823, true-val::0.5833
>> sample mean:0.0704, true-val::0.0708
ガンマ関数と置換積分
Problem
$x\in(0,1)$で確率密度関数が以下のように与えられたとする
\[f(x) = -\log(x)\](1) $\int^1_0 f(x)dx = 1$になることを確かめよ
(2) $X$について, 期待値と分散を求めよ
解答
また,
\[\begin{align*} \Pr(X \leq x) &= \bigg[-x\log(x)\bigg]^x_0 + \int^x_0 dx\\[3pt] &= -x\log(x) + x \end{align*}\]期待値は, $-\log(x) = u$と変数変換をおくと
\[\begin{align*} \mathbb E[X] &= \int^1_0 -x\log(x)dx\\[3pt] &= \int^0_\infty u \exp(-u)\exp(-u)\times (-1) du \\[3pt] &= \int^\infty_0 u \exp(-2u) du \\[3pt] &= \int^\infty_0 \frac{1}{2}t\exp(-t)\frac{1}{2}dt\\[3pt] &= \frac{1}{4}\int^\infty_0t^{2-1}\exp(-t)dt\\[3pt] &= \frac{1}{4} \Gamma(2)\\[3pt] &= \frac{1}{4} \end{align*}\]分散については$\mathbb E[X^2]$が求まれば良いので
\[\begin{align*} \mathbb E[X^2] &= \int^1_0 -x^2\log(x)dx\\[3pt] &= \int^\infty_0 u \exp(-3u) du \\[3pt] &= \int^\infty_0 \frac{1}{3}t \exp(-t)\frac{1}{3} du \\[3pt] &= \frac{1}{9} \end{align*}\]従って,
\[\text{Var}(X) = \frac{1}{9} - \bigg(\frac{1}{4}\bigg)^2 = \frac{7}{144}\]三角関数と置換積分
基本問題
Problem
以下の定積分を求めよ
\[\int^{\pi}_0\cos\sqrt{x}dx\]解答
$t = \sqrt{x}$とおいて
\[\begin{align*} \int^{\pi}_0\cos\sqrt{x}dx &= \int^{\sqrt{\pi}}_02t\cos t dt\\[3pt] &= [2t\sin t]^{\sqrt{\pi}}_0 - \int^{\sqrt{\pi}}_0 2\sin t dt\\[3pt] &= [2t\sin t]^{\sqrt{\pi}}_0 + [2\cos t]^{\sqrt{\pi}}_0\\[3pt] &= 2(\sqrt{\pi}\sin \sqrt{\pi} + \cos \sqrt{\pi} - 1) \end{align*}\]$\arcsin$と置換積分
Problem
$x\in(0,1)$で確率密度関数が以下のように与えられたとする
\[f(x) = \frac{2}{\pi\sqrt{1-x^2}}\]- $\int^1_0 f(x)dx = 1$になることを確かめよ
- $X$について, 期待値と分散を求めよ
解答
であることに留意すると
\[\begin{align*} \int^1_0 \frac{2}{\pi\sqrt{1-x^2}}dx &= \frac{2}{\pi}[\arcsin(x)]^1_0\\ &= \frac{2}{\pi} \bigg(\frac{\pi}{2} - 0\bigg)\\ &= 1 \end{align*}\]期待値は$x=\sin\theta$を用いて
\[\begin{align*} \mathbb E[X] &= \int^1_0\frac{2x}{\pi\sqrt{1-x^2}}dx\\[3pt] &= \frac{2}{\pi}\int^{\pi/2}_0 \frac{\sin\theta}{\cos\theta}\cos\theta d\theta\\[3pt] &= \frac{2}{\pi}\int^{\pi/2}_0\sin\theta d\theta\\ &= \frac{2}{\pi} \end{align*}\]$\mathbb E[X^2]$は
\[\begin{align*} \mathbb E[X^2] &= \int^1_0\frac{2x^2}{\pi\sqrt{1-x^2}}dx\\[3pt] &= \frac{2}{\pi}\int^{\pi/2}_0 \frac{\sin^2\theta}{\cos\theta}\cos\theta d\theta\\[3pt] &= \frac{2}{\pi}\int^{\pi/2}_0\sin^2\theta d\theta\\[3pt] &= \frac{1}{\pi}\int^{\pi/2}_0(1 - \cos2\theta)d\theta\\[3pt] &= \frac{1}{2} \end{align*}\]従って,
\[\text{Var}(X) = \frac{1}{2} - \frac{4}{\pi^2}\]シミュレーションは可能なのか？
Acceptance-Rejection samlingでは
\[M = \sup\left\{\frac{f(x)}{f_0(x)}\right\} < \infty\]という条件が必要ですが
\[\lim_{x\to1}f(x) = \frac{2}{\pi\sqrt{1-x^2}}=\infty\]なので, Acceptance-Rejection samlingでsimulationするのは難しそうです.
逆関数法
今回の累積密度関数は
\[\Pr(X\leq x) = \frac{2}{\pi} \arcsin(x)\]なので逆関数が
\[\text{ICDF}(x) = \sin\bigg(\frac{x\pi}{2}\bigg)\]とわかるので, 逆関数法でシミュレーションすることができます.
変数変換法
について$x = \sin\theta$と変数変換し, $\theta$についての密度関数$g$を考えてみると $\theta\in(0, \pi/2)$ について
\[\begin{align*} g(\theta) &= \frac{2}{\pi\sqrt{1-\sin^2\theta}}\cos\theta\\ &= \frac{2}{\pi} \end{align*}\]従って,
\[\theta\sim\text{Unif}(0, \pi/2)\]に従う$\theta$に対して$x = \sin(\theta)$と変換してシミュレーションデータの生成もできそうです.
Column: 変数変換は一意ではない
\[\theta\sim\text{Unif}(0, \pi)\]について$x=\sin(\theta)$と考えたときの$x$の密度関数を考えてみます.
\[\begin{align*} \Pr(X\leq x) &= \Pr(\sin\theta \leq x)\\ &= \Pr(\theta \in (0, \arcsin(x))) + \Pr(\theta \in (\pi-\Pr(\theta \in (0, \arcsin(x))), \pi))\\ &= \frac{2\arcsin(x)}{\pi} \end{align*}\]従って,
\[f(x) = \frac{2}{\pi\sqrt{1-x^2}}\]$\theta\sim\text{Unif}(0, \pi)$からの変数変換でも今回の目的密度関数が再現できる.
Python simulation
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import numpy as np
from scipy import stats
def cdf(x):
return (2 / np.pi) * np.arcsin(x)
def pdf(x):
return 2 / (np.pi * np.sqrt(1 - x**2))
## DGP
np.random.seed(42)
N = 20000
u = np.random.uniform(0, np.pi/2, N)
u2 = np.random.uniform(0, 1, N)
x = np.sin(u)
x2 = np.sin(u2 * np.pi/2)
## ks test
print(stats.kstest(x, cdf))
>>> KstestResult(statistic=0.0038036929545425258, pvalue=0.9333433180624549, statistic_location=0.8256177094009344, statistic_sign=1)
print(stats.kstest(x2, cdf))
>>> KstestResult(statistic=0.004595894932220124, pvalue=0.7902584909796191, statistic_location=0.8961216642195171, statistic_sign=1)
期待値分散についても
1
2
3
4
5
6
7
8
9
print("sample mean:{:.4f}, true-val::{:.4f}".format(np.mean(x), 2/np.pi))
print("sample mean:{:.4f}, true-val::{:.4f}".format(np.mean(x2), 2/np.pi))
print("sample mean:{:.4f}, true-val::{:.4f}".format(np.var(x), 1/2 - 4/np.pi**2))
print("sample mean:{:.4f}, true-val::{:.4f}".format(np.var(x2), 1/2 - 4/np.pi**2))
>>> sample mean:0.6360, true-val::0.6366
>>> sample mean:0.6364, true-val::0.6366
>>> sample mean:0.0946, true-val::0.0947
>>> sample mean:0.0944, true-val::0.0947
$\arccos$と置換積分
Problem
\[\int^1_0 x^2\arccos x dx\]解答
であるので
\[\begin{align*} &\int^1_0 x^2\arccos x dx\\[3pt] &=-\int^{\pi/2}_0t\cos^2t(-\sin t)dt\\[3pt] &= -\bigg[t\frac{\cos^3t}{3}\bigg]^{\pi/2}_0 + \int^{\pi/2}_0\frac{\cos^3t}{3}dt \end{align*}\]このとき, 積和の公式より
\[\cos^3t = \frac{3\cos x + \cos 3x}{4}\]従って
\[\begin{align*} &-\bigg[t\frac{\cos^3t}{3}\bigg]^{\pi/2}_0 + \int^{\pi/2}_0\frac{\cos^3t}{3}dt \\[3pt] &= 0 + \bigg[\frac{\sin x}{4}\bigg]^{\pi/2}_0 + \bigg[\frac{\sin 3x}{36}\bigg]^{\pi/2}_0 \\[3pt] &= \frac{2}{9} \end{align*}\]$\arctan$と置換積分
Problem
$a, b > 0$としたとき, 次の定積分の値を計算せよ
\[I = \int^{\pi/2}_{0}\frac{1}{(a^2\cos^2\theta + b^2\sin^2\theta)^2}d\theta\]解答
ここで,
\[\begin{align*} &t = \bigg(\frac{b}{a}\tan\theta\bigg)\\[3pt] &\frac{d\theta}{dt}= \frac{a}{b}\cos^2\theta \end{align*}\]とおくと
\[\begin{align*} I &= \int^{\pi/2}_{0}\frac{\frac{1}{a^4\cos^2\theta} + \frac{\tan^2\theta}{a^4\cos^2\theta}}{(1 + \frac{b^2}{a^2}\tan^2\theta)^2}d\theta\\[3pt] &= \frac{1}{a^3b^3}\int^{\infty}_0\frac{b^2 + a^2t^2}{(1+t^2)^2}dt \end{align*}\]このとき
\[I_1 = \int^{\infty}_0\frac{b^2 + a^2t^2}{(1+t^2)^2}dt\]に対して, 以下のような変数変換を実施する
\[t = \frac{1}{x}\]すると
\[\begin{align*} \int^{\infty}_0\frac{b^2 + a^2t^2}{(1+t^2)^2}dt &= \int_{\infty}^0\frac{b^2 + a^2x^{-2}}{(1+x^{-2})^2}(-x^{-2})dx\\[3pt] &= \int^{\infty}_0\frac{b^2 + a^2x^{-2}}{(1+x^{-2})^2}\frac{x^2}{x^4}dx\\[3pt] &= \int^{\infty}_0\frac{x^2b^2 + a^2}{(x^2+1)^2}dx \end{align*}\]従って
\[\begin{align*} 2I_1 &= \int^{\infty}_0\frac{b^2 + a^2t^2}{(1+t^2)^2}dt + \int^{\infty}_0\frac{x^2b^2 + a^2}{(x^2+1)^2}dx\\[3pt] &= \int^\infty_0\frac{a^2+b^2}{(1+t^2)}dt\\[3pt] &= (a^2+b^2)[\arctan(t)]^\infty_0\\[3pt] &= (a^2+b^2)\frac{\pi}{2} \end{align*}\]従って,
\[I = \frac{a^2+b^2}{a^3b^3}\frac{\pi}{4}\]Appendix: Acceptance-Rejection Sampling
Acceptance-Rejection Sampling Algorithm
Input: target density $f$, proposal density $f_0$, constant $M$ satisfying the follwoing
\[M = \sup\left\{\frac{f(x)}{f_0(x)}\right\}\]Output: one realization from $f$
Repeat: generate $y$ from $f_0$ and $u$ from $\text{Unif}(0, 1)$, independently
Until
\[u \leq \frac{f(y)}{Mf_0(y)}\]Return the last $y$
- $M$が大きくなると棄却する割合が大きくなり非効率なサンプリング法になる
- $M<\infty$を満たすためには, 提案分布の密度$f_0(x)$は目標分布の密度$f(x)$よりも裾が厚くなる必要がある
- 例: 提案分布の密度がコーシー分布で目標分布が正規分布のときは機能するが, 逆の場合は機能しない
References
- Lecture note 2007 by Karl Sigman > Acceptance-Rejection Method
- Ryo’s Tech Blog > 半円周上に一様分布する点の高さの分布
- KIT数学ナビゲーション > 三角関数の1次化のための公式
(注意：GitHub Accountが必要となります）