Pyhonで確率分布を解説・実装してみた

統計学における確率についてまとめました。数式も技術的なことも避けずに解説・実装しています。

バージョンは以下の通りです。

❯ python -V                                                                                 
Python 3.7.4

❯ conda -V
conda 4.9.1

❯ jupyter --version
jupyter core     : 4.5.0
jupyter-notebook : 6.0.1
qtconsole        : 4.5.5
ipython          : 7.8.0
ipykernel        : 5.1.2
jupyter client   : 5.3.3
jupyter lab      : 1.1.4
nbconvert        : 5.6.0
ipywidgets       : 7.5.1
nbformat         : 4.4.0
traitlets        : 4.3.3

Jupyter Labをつかっています。以下のライブラリをインポートしておきます。

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()

確率の基礎

まずは確率の基礎。

サイコロを一回投げて出た目を観測する時、

1,2,3,4,5,6の目が出る結果のことを標本点と言い、

標本点をまとめた集合$\Omega = \{1,2,3,4,5,6\}$のことを、標本空間といいます。

標本空間$\Omega$の部分集合を事象といいます。

偶数の目が出る事象は、$\{2,4,6\} \subset \Omega$です。

特になにも起きないという事象を空事象と言い、$\phi$で表します。

標本空間がk個の要素からなるとき、事象は$2^k$個存在します。

コイン投げをして表を1、裏を0とするとき、標本空間は$\Omega = \{0,1\}$であり、事象（部分集合）を全て列挙すると、

$$\phi, \{0\}, \{1\}, \{0,1\}$$

となります。

事象Aの補集合を余事象といい、$A^c$と書き、$A^c = 1 - A$となります。

事象Aと事象Bが同時に起きるという事象をAとBの積事象といい、$A \cap B$と書き、

事象Aと事象Bのいずれかが起きるという事象をAとBの和事象といい、$A \cup B$と書きます。

条件つき確率

条件付き確率は
$$ P(A|B) = \frac{P(A \cap B)}{P(B)} $$
で表せます。これはBという事象が起こることが確定している時にAが起こる確率です。

$ \frac{P(A \cup B)}{\Omega} $ から、$ \frac{P(A \cup B)}{P(B)} $になるので、分母が小さくなった分だけ値が大きくなります。

サイコロの例をPythonで実装してみます。

fig = plt.figure(figsize=(8,4))
axes = fig.add_subplot(111)

np.random.seed(1)
s = np.random.randint(1,7,10000)

axes.hist(s,bins=6)
axes.set_title('サイコロの確率分布')
axes.set_xlabel('サイコロの目')
axes.set_ylabel('度数')

{idx:len(s[s==idx]) for idx in range(1,7)}

{1: 1648, 2: 1609, 3: 1693, 4: 1679, 5: 1635, 6: 1736}

それぞれの出目の出方は一様であることがわかりました。

set_a = s[s<=5]
prob_a = len(set_a) / len(s)
prob_a

0.8373

サイコロを振って5以下の目が出る確率$P(A)$は0.84でした。

set_b = s[s%2==0]
prob_b = len(set_b) / len(s)
prob_b

0.5024

サイコロを振って偶数が出る確率$P(B)$は$\frac{1}{2}$でした。

set_a_and_b = s[(s<=5) & (s%2==0)]
prob_a_and_b = len(set_a_and_b) / len(s)
prob_a_and_b

0.3311

サイコロを振って5以下かつ偶数が出る確率$P(A \cap B)$は0.33でした。

では、$ P(A|B) = \frac{P(A \cap B)}{P(B)} $を計算してみます

prob_a_and_b / prob_b

0.670514378290806

事象Bが起こった場合に事象Aも起こる確率$P(A|B)$は0.67でした。

このように、事象Bが起こった場合に限定すると確率が高くなることがわかります。

では確率分布を実装します。

def f(n): # 条件付き確率関数
  p_arr = []
  for _ in range(10000):
    s = np.random.randint(1,7,10000)
    set_a = s[s<=n]
    prob_a = len(set_a) / len(s)
    set_b = s[s%2==0]
    prob_b = len(set_b) / len(s)
    set_a_and_b = s[(s<=n) & (s%2==0)]
    prob_a_and_b = len(set_a_and_b) / len(s)
    p_arr.append(prob_a_and_b / prob_b)
    return np.mean(p_arr)
  
g_data = {n:f(n) for n in range(1,7)}
  
fig = plt.figure(figsize=(8,4))
axes = fig.add_subplot(111)

x = g_data.keys()
y = g_data.values()

axes.bar(x,y)
axes.set_title('条件付き確率分布')
axes.set_xlabel('サイコロの目')
axes.set_ylabel('確率')

確率変数と確率分布

たとえば、サイコロの目1,2,3,4,5,6があって、重心がサイコロのちょうど中心にある場合、それぞれの出目が出る確率は$\frac{1}{6}$ですよね。

このように、それぞれの結果とその確率の組み合わせ両方が存在するもののことを、確率変数といいます。

そしてその値と確率の分布のことを確率分布と言います。

離散確率変数

サイコロの出目のように、出る結果とその確率が存在するものは確率変数と考えることができます。Xをサイコロの出目とすると、

$$ P(X = x) = \frac{1}{6}, x = 1,2,…,6 $$

となります。

サイコロのようにとびとびの値（離散値）を持つ確率変数を離散確率変数といいます。

離散確率変数の期待値は、

$$ E(X) = \sum_{x} x P(X = x) $$

分散は

$$ V(X) = \sum_{x} (x - \mu)^2 P(X = x) $$

となります。

離散一様分布

n個の相異なる値$a_1,a_2,… , a_n$を等しい確率$\frac{1}{n}$でとる確率変数が従う確率分布を、$\{a_1,a_2, … , a_n\}$を台に持つ離散一様分布といいます。

ではサイコロの確率分布を見てみましょう。

fig = plt.figure(figsize=(8,4))
axes = fig.add_subplot(111)

sample_points = np.arange(1,7) #標本点
prob = np.ones(len(sample_points))/6 #標本点にインデックスで対応した確率
random_variable = [sample_points, prob] #確率変数

axes.bar(random_variable[0], random_variable[1])
axes.set_ylim(0,1)
axes.set_title('サイコロの確率分布')
axes.set_xlabel('度数')
axes.set_ylabel('サイコロの目')

このような確率分布になります。どの目が出る確率も$\frac{1}{6}$になっています。

離散確率変数の標本とその確率を対応づける関数を、確率質量関数、あるいは単に確率関数といい、

離散一様分布の場合、確率質量関数は、

$$P(X = k) = \frac{1}{N} (k = 1,2,3, ... N)$$

あるいは

$$P(X = k) = \frac{1}{b - a + 1} (k = a, a+1,a+2, ... b)$$

で表すことができます。サイコロの場合、xの取りうる値は{1, 2, 3, 4, 5, 6}であり、aは1、bは6です。

一様分布の期待値は$\frac{a + b}{2}$、分散は$\frac{(b - a)^2}{12}$になります。

幾何分布

離散一様分布では標本点は有限個しかありません。たとえば確率p > 0で起きる事象Aが、n回目に初めて起きる確率は、Aが起きるまでの回数をXとしたとき、

$$P(X = n) = p(1 - p)^{n - 1}$$

で表すことができます。

この確率分布は、幾何分布といいます。たとえばサイコロの目3が10回サイコロを振って初めて出る確率は、

$$ P(X = 10) = \mathrm{\frac{1}{6}(1 - \frac{1}{6})}^{10 - 1} = 0.0323011...$$

となります。

ではこの幾何分布を実装してみます。

fig = plt.figure(figsize=(8,4))
axes = fig.add_subplot(111)

def f(num, prob): # 幾何分布の確率関数
  return (prob)*(1 - prob)**(num-1)

nums = np.arange(1,101) # 1~100（サイコロの目xが初めて出るまでの回数の配列）
probs = [f(num,1/6) for num in nums]

axes.bar(nums, probs)
axes.set_xlim(0,30)
axes.set_ylim(0,0.2)
axes.set_title('サイコロの幾何分布')
axes.set_xlabel('サイコロの目が初めて出るまでの回数')
axes.set_ylabel('確率')

幾何分布の期待値は$\frac{1}{p}$、分散は$\frac{1 - p}{p^2}$です。

幾何分布は離散確率変数でも連続確率変数でも同じように記述できます。

二項分布

一回の試行で確率$p$で起きる事象が$n$回の試行で$k$回起こる確率は、

$$ {}_n C_k p^k(1 - p)^k $$

で表すことができます。これを二項分布と呼びます。

二項分布の期待値は$np$、分散は$np(1 - p)$になります。

また、後述する正規分布は二項分布の極限として捉えることができます。つまり、二項分布$Bi(n,p)$のサンプル数nが十分に大きい時、正規分布$N(np, np(1-p))$によく近似するということです。これを、ド・モアブル＝ラプラスの定理と呼びます。

ポアソン分布

一定の期間に平均$\lambda$回起きる事象（独立かつランダムに起きる）が$k$回起きる確率はポアソン分布に従います。ポアソン分布の確率質量関数は

$$P(X=k)=\dfrac{\lambda^ k e^{-\lambda}}{k!}$$

と表すことができます。

二項分布の$n$が十分に大きく$p$が非常に小さい時に$np=一定$として考え、$np=\lambda$と置いたとき、ポアソン分布になります。

例えば、ある地域における一年間の交通事故の回数や、一定時間に届くメールの数などがポアソン分布に従います。

左に山がある非対称な確率分布です。

期待値は$\lambda$、分散も$\lambda$になります。

負の二項分布（またはパスカル分布）

一回の試行で確率$p$で起こる事象を$r$回行って失敗回数$X$が$k$である確率

$$P(X=k)={}_{r+k-1}C_k p^ r (1-p)^k$$

という確率質量関数で表すことができ、$NB(r,p)$と略記します。

$r = 1$の時は幾何分布に一致します。また負の二項分布は一般化線形モデルで用いられます。

期待値は$\frac{k}{p}$、分散は$\frac{k(1-p)}{p^2}$になります。

連続確率変数

たとえばネジの直径が$\sqrt{2}$mmである確率は$\frac{1}{\infty}$なので0です。そこで、$1.41 \leq X \leq 1.45$のような幅を持たせて考えます。多くの連続確率変数は、適当な関数$f(x)$を用いて、

$$ P(a \leq X \leq b) = \int_{a}^{b}f(x)dx  $$

と表すことができます。

連続確率変数の期待値は

$$ E(X) = \int_{a}^{b}x f(x)dx$$

分散は

$$ V(X) = \int_{a}^{b}(x - \mu)^2 f(x)dx$$

となります。

この時の関数$f(x)$を確率密度関数といい、分布の種類によって異なります。

たとえば、正規分布の確率密度関数は、

$$ f(x) = \frac{1}{\sqrt{2\pi}\sigma}e^{- \frac{(x - \mu)^2}{2\sigma^2}} $$

確率変数Xの対数をとると正規分布する、対数正規分布の確率密度関数は

$$ f(x) = \frac{1}{\sqrt{2\pi}\sigma x}e^{- \frac{(\log x - \mu)^2}{2\sigma^2}} $$

連続一様分布の確率密度関数は

$$ f(x) = \begin{cases} \frac{1}{b - a} (a \leq x \leq b) \\ 0　(その他) \end{cases} $$

などのように表すことができます。

累積分布関数

Xの累積分布関数を
$$ F(x) = P(X \leq x) $$
で表します。これは確率変数Xがx以下の値をとる確率です。

確率密度をもつ確率分布の累積分布関数は、
$$ F(x) = \int_{-\infty}^{x}f(s)ds $$
で表します。$F(x)$が微分可能であれば、$F'(x) = f(x)$が成り立ちます。つまり累積分布関数を微分すると、確率密度関数になるということです。

正規分布の、160cm以下の割合をPythonで求めます。

from scipy import integrate

def f(s, x): #　累積分布関数
  return integrate.quad(s,-np.inf,x)

y = lambda x: norm.pdf(x,170,10)
ie = f(y,160)

print(ie)

(0.15865525393145707, 1.0511793462042972e-11)

タプルの左の値が積分値、右の値が推定誤差です。

正規分布

身長のような連続的な確率変数を取る場合、ちょうど170.12345cmの人の割合はものすごく小さなものになります。ほぼ0になるでしょう。そこで、170cmから171cmの間の人の割合を出すようにすれば、その割合は使える指標になります。

まず、度数分布を見てみましょう。

fig = plt.figure(figsize=(8,4))
axes = fig.add_subplot(111)

np.random.seed(1)
x = np.random.normal(170, 10, 1000)
axes.hist(x,bins=60)
axes.set_xlim(130,210)
axes.set_title('身長の度数分布')
axes.set_xlabel('身長')
axes.set_ylabel('度数')

身長は正規分布に従う確率変数なので、np.random.normalで平均170、標準偏差10で1000個の正規分布に従うデータをランダムに取得します。

この一つ一つの棒の横幅が身長の区間、高さがその区間に属するデータの数になります。

ではnp.random.normalでは、どんなことを実際にしているのかというと、確率密度関数を使っています。正規分布に従う確率密度関数は以下のような式で表されます。

$$ f(x) = \frac{1}{\sqrt{2\pi}\sigma}e^{- \frac{(x - \mu)^2}{2\sigma^2}} $$

では、Pythonで実装してみましょう。

from scipy.stats import norm

fig = plt.figure(figsize=(8,4))
axes = fig.add_subplot(111)
x = np.linspace(100,300,1000)
#正規分布の確率密度関数,平均μ: 170、標準偏差σ：10
y = norm.pdf(x,170,10)

axes.plot(x,y)
axes.set_xlim(130,210)
axes.set_title('確率密度関数の例')
axes.set_xlabel('$x$')
axes.set_ylabel("$f(x)$")

xは100から300までの等差数列1000個です。そのxの数字それぞれにたいする、

$$ f(x) = \frac{1}{\sqrt{2\pi}\sigma}e^{- \frac{(x - \mu)^2}{2\sigma^2}} $$

がy軸です。scipy.statsモジュールのnorm.pdfでこれを実装しています。

上述しましたが、正規分布は二項分布の極限として捉えることができ、二項分布$Bi(n,p)$のサンプル数nが十分に大きい時、正規分布$N(np, np(1-p))$はよく近似します。これを、ド・モアブル＝ラプラスの定理と呼びます。

連続一様分布

閉区間[a,b]を台に持つ連続一様分布は、以下のような確率密度関数に従います。

$$ f(x) = \begin{cases} \frac{1}{b - a} (a \leq x \leq b) \\ 0　(その他) \end{cases} $$

Pyhonで実装してみます。

# np.random.uniform: 下限 low、上限 high、要素数 size の一様分布に従う配列を生成する関数

def f(x,a,b): # 連続一様分布の確率密度関数
  if (x > a) & (x < b):
    return 1/(b-a)
  else:
    return 0

fig = plt.figure(figsize=(8,4))
axes = fig.add_subplot(111)

x = np.random.uniform(0, 300, 10000)
y = [f(x_i,160,180) for x_i in x]

axes.scatter(x,y,marker='.')
axes.set_title('連続一様分布の確率密度関数')
axes.set_xlabel('$x$')
axes.set_ylabel("$f(x)$")

指数分布

一定時間に届くメールの数など、独立な事象が一定期間の間に起きる回数の分布はポアソン分布に従います。メールが届いてから次にメールが届くまでの期間は指数分布にしたがいます。ある事象がある期間内に起きる回数の平均が$\lambda$ですると、指数分布の確率密度関数は

$$f(x) = \lambda e^{-\lambda x} (x \geq 0)$$

と表せます。

期待値は$E(X) = \frac{1}{\lambda}$、分散は$V(X) = \frac{1}{\lambda^2}$です。

コーシー分布

期待値が存在しないため分散も存在せず、また中心極限定理にも当てはまらないという珍しい分布です。以下のような確率密度関数に従います。

$$ f(x) = \frac{1}{\pi} \frac{\gamma}{\gamma^2 + (x - x_0)^2} $$

$x_0$は最頻値を与える位置母数、$\gamma > 0$は尺度母数です。

コーシー分布は一見すると正規分布と似ています。

正規分布は外れ値を取る確率が低く（裾が軽い）、コーシー分布は「外れ値を取る確率が高い（裾が重い）ため、

外れ値が重要な意味を持つような状況、例えば「証券価格の大暴落が起こる確率はそこまで低くなく，正規分布のような裾の軽い分布でモデル化するのは難しいためコーシー分布を使うことがあります。

ワイブル分布

生命表の10万人あたり死亡数や、機械部品の寿命などはワイブル分布に従います。

ワイブル分布の確率密度関数は

$$ f(x) = \frac{k}{\lambda} (\frac{x}{\lambda})^{k - 1} e ^{- (\frac{x}{\lambda})^k} (x > 0)$$

と表すことができます。

$k > 0$は形状母数（母斑）、$\lambda > 0$は位置母数（母斑）と呼ばれます。確率分布の特性を表す定数のことを母数（母斑）と言います。例えば正規分布の母数は平均と分散、二項分布の母数はサンプルサイズと確率です。形状母数は形状を特定するための定数、位置母数は位置を特定するための定数、後述のガンマ分布の確率密度関数に現れる尺度母数は尺度を特定するための定数です。

独立な確率変数の和の分布

離散確率変数の場合

独立な二つの離散確率変数$X$、$Y$の和$S = X + Y$の分布は以下のように表せます。

$$ P(S = s) = P(X + Y = s) = \sum_{x_i + y_j = s} P(X = x_i)P(Y = y_j) $$

この$X$が二項分布$Bi(m,p)$、$Y$が$Bi(n,p)$として計算すると

$$ P(S = k) = {}_{m+n} C_k p^k (1 - p)^{m+ n -k} $$

が得られます。これは$Bi(m+n, p)$であることがわかります。これを二項分布は再生性があるといいます。

先ほどの計算は高度な計算技術が必要であり、かつ大変です。これを積率母関数を使うと、シンプルに求めることができます。たとえば二項分布の積率母関数は

$$ M_X(t) = E(e^{tX}) = (pe^t + (1 - p))^m$$

です。$S = X + Y$の積率母関数を計算すると、

\begin{eqnarray}
E(e^{t(X+Y)}) &=& E(e^{tX})E(e^{tY}) \\
&=& (pe^t + (1 - p))^m (pe^t + (1 - p))^n \\
&=& (pe^t + (1 - p))^{m+n} \\
\end{eqnarray}

となり、これは$Bi(m+n, p)$であることがわかります。

連続確率変数の場合

独立な二つの連続確率変数$X$、$Y$のそれぞれの確率密度関数を$f(x)$、$g(x)$として、その和$S = X + Y$の累積分布関数は以下のように計算できます。

\begin{eqnarray}
H(S = s) &=& P(X + Y = s) \\
&=& \sum_{x_i + y_j = s} P(X = x_i)P(Y = y_j) \\
&=& \int \int_{x + y \leq s} f(x)g(y)dxdy \\
&=& \int_{-\infty}^{\infty} (\int_{-\infty}^{s - y} f(x) d(x))g(y)dy \\
&=& \int_{-\infty}^{\infty} F(s - y)g(y)dy \\
\end{eqnarray}

確率密度関数は、この累積分布関数を微分したものなので

\begin{eqnarray}
h(s) = H'(s) &=& \int_{-\infty}^{\infty} \frac{d}{ds} F(s - y)g(y)dy \\
&=& \int_{-\infty}^{\infty}f(s - y)g(y)dy
\end{eqnarray}

となります。

正規分布の積率母関数は

$$ e^{\mu t + \frac{\sigma^2 t^2}{2}} $$

であり、$N(\mu_1, \sigma_1^2)$）である確率変数$X$、$N(\mu_2, \sigma_2^2)$である確率変数$Y$に対し、$X + Y$の積率母関数は、

\begin{eqnarray}
E(e^{t(X+Y)}) &=& E(e^{tX})E(e^{tY}) \\
&=& e^{\mu_1 t + \frac{\sigma_1^2 t^2}{2}} e^{\mu_2 t + \frac{\sigma_2^2 t^2}{2}} \\
&=& e^{(\mu_1 + \mu_2) t \frac{1}{2}(\sigma_1^2 + \sigma_2^2) t^2}
\end{eqnarray}

となり、これは$N(\mu_1 + \mu_2, \sigma_1^2 + \sigma_2^2)$の積率母関数です。これを正規分布が再生性を持つと言います。

ガンマ分布

一定期間に1回起きると期待されるランダムな事象が複数回起きるまでの時間の分布です。確率密度関数は

$$ f(x) = \frac{1}{\Gamma(k) \theta^k} x^{k-1} e^{- \frac{x}{\theta}} (x > 0)$$

です。この$\Gamma(k)$はガンマ関数であり、以下のように表せます。

$$\Gamma(k) = \int_{0}^{\infty} x^{k - 1} e^{-x}dx$$

また$ k > 0 $は形状母数、 $\theta > 0 $は尺度母数です。

ガンマ分布の期待値は$ E(X) = k\theta $、分散は$ V(X) = k\theta^2 $です。

後述のアーラン分布やカイ二乗分布は、このガンマ分布の特別な形です。

アーラン分布

アーラン分布はガンマ分布で形状母数を正の整数に限定したものです。確率密度関数は、

$$ f(x) = \begin{cases} \frac{1}{(k - 1)! \theta^k} x^{k - 1} e^{- \frac{x}{\theta}} (x > 0) \\ 0　(x \leq 0) \end{cases} $$

です。ガンマ分布と同様アーラン分布の期待値は$ E(X) = k\theta $、分散は$ V(X) = k\theta^2 $です

カイ二乗分布

カイ二乗分布は「母分散の区間推定」や「適合度の検定」、「独立性の検定」を行う際に使われます。自由度$k$のカイ二乗分布の確率密度関数は

$$ f_k(x) = \begin{cases} \frac{1}{2^{\frac{k}{2}}\Gamma(\frac{k}{2})} x^{\frac{k}{2} - 1} e^{- \frac{x}{2}} (x > 0) \\ 0　(x \leq 0) \end{cases}$$

に従います。

カイ二乗分布に関する定理をいくつか紹介します。

定理：（カイ二乗分布の再生性）カイ二乗分布に従う確率変数$X$、$Y$がそれぞれ自由度$k_1$、$k_2$である時、確率変数$X+Y$は自由度$k_1 + k_2$のカイ二乗分布に従う

定理：標準正規分布$N(0,1)$に従う確率変数$X_1, X_2, ... X_k$に対し、$ Z = \sum_{j =1}^{k} X_j^2 $は自由度$k$のカイ二乗分布に従う

多変量確率分布

離散確率変数の場合

離散確率変数$X_1, X_2, … X_n$を並べて作られる確率変数$X = (X_1, …, X_n)$の確率分布
$$ P^X(B) = P(X \in B) $$
により、Xの確率分布を定義します。$X \in B$は、「XはBに含まれる」という意味です。

これを$X_1, X_2,…X_n$の結合分布といいます。BはXの取りうる値全体からなる集合の部分集合です。つまり、複数の確率変数をまとめてひとつの確率変数と見た時の確率分布を結合分布といいます。

$X_1=x_1,X_2=x_2,…,X_n=x_n$となる確率を、$x_1,x_2,…,x_n$の関数と見たもの
$$ P(x_1,x_2,…,x_n) = P(X_1=x_1,X_2=x_2,…,X_n=x_n) $$
を結合確率関数といいます。

たとえば、2枚のコイン1、コイン2を投げて表なら1、裏なら0と書き、コイン1の値を確率変数と見たものを$X_1$、コイン2の値を確率変数と見たものを$X_2$とします。

その時の$(X_1, X_2)$の取る値は$(0,0)(0,1)(1,0)(1,1)$の4通りで各々の確率は$\frac{1}{4}$です。この確率分布が$X_1$と$X_2$の結合分布です。

二つの離散的確率変数X,Yに対して、
$$ P^X(x) = \sum_{y}^{}P(x,y) $$
を、Xの周辺分布といいます。$X_1$,$X_2$,…$X_n$に対し、$X_j$の周辺分布は、結合確率関数の$X_j$以外の変数について和をとったものです

連続確率変数の場合

離散確率変数と同様に、確率ベクトル$X = (X_1,X_2,…,X_n)$に対し、集合Bに対して、
$$ P^X(B) = P(X \in B) $$
でXの確率分布を定義し、$X_1, X_2,…X_n$の結合分布といいます。

連続確率変数の場合、$ P^X(B) $の確率密度関数を結合確率密度関数といいます。つまり

$$ P^X(B) = \int…\int_B f(x_1,x_2,…,x_n)dx_1,dx_2,…,dx_n $$
となる$f(x_1,x_2,…,x_n)dx_1,dx_2,…,dx_n $を結合確率密度関数といいます。

$\int…\int_B$の部分は、$B$全体にわたる積分を表します。

連続確率変数の周辺分布は、離散確率変数における和を積分にして定義します。つまり、$f(x,y)$を$X$、$Y$の結合確率密度関数とした時、
$$ f^X(x) = \int_{-\infty}^{\infty}f(x,y)dy $$
が、Xの周辺確率密度関数です。