ある工場で製造された部品の強度を測定しているとします。強度が基準値を下回る部品は出荷前に取り除かれるため、手元のデータには「基準値以上の強度を持つ部品」しか含まれていません。このデータに通常の正規分布を当てはめると、平均や分散の推定値は実際の母集団のものとは異なってしまいます。
「観測されたデータが、元の分布のある範囲に限定されている場合、どのように分析すべきか？」
この問いに答えるのが切断分布（truncated distribution）です。切断分布は、ある連続分布から特定の範囲のデータだけを取り出したときに得られる分布です。元の分布の一部を「切り取った」ように見えることから、この名前がつけられています。
切断分布を理解すると、以下のような応用が開けます。
- 品質管理: 検査で選別された部品の母集団パラメータ推定
- 経済学: トービットモデル（所得、消費データなどでゼロ以下が観測されない場合）
- 生存分析: 観測開始時点以降に発症した患者のみの分析（左切断）
- 環境科学: 検出限界以下の汚染物質濃度のモデリング
本記事の内容
- 切断分布の一般的な定義と直感的理解
- 切断正規分布の確率密度関数・モーメント
- 切断と打ち切り（censoring）の違い
- トービットモデルへの応用
- Pythonによる実装と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
切断分布とは
日常のアナロジー
大学の入試を想像してください。受験者全体の成績は正規分布に近い形をしていますが、合格者だけを見ると分布の形が変わります。合格ラインより低い成績の人がいないため、分布は合格ラインで「切り取られ」、残った部分の面積が1になるように高さが調整されます。
これが切断分布の本質です。元の分布から特定の範囲を取り出し、確率の合計が1になるように再正規化した分布です。
一般的な定義
確率変数 $X$ が確率密度関数 $f(x)$ と累積分布関数 $F(x)$ を持つとき、区間 $(a, b)$ への切断分布の確率密度関数は
$$ \begin{equation} f_T(x) = \begin{cases} \frac{f(x)}{F(b) – F(a)} & \text{if } a < x < b \\ 0 & \text{otherwise} \end{cases} \end{equation} $$
で定義されます。ここで $F(b) – F(a) = P(a < X < b)$ は元の分布で区間 $(a, b)$ に入る確率であり、これが正規化定数の役割を果たします。
式の構造は非常に直感的です。元の密度 $f(x)$ をそのまま使いますが、切断範囲外の確率をゼロにしたことで失われた確率の分だけ、全体を $1/(F(b) – F(a))$ 倍して確率の合計が1になるようにしています。
切断は範囲によって3つのタイプに分類されます。
- 左切断（left truncation）: $a$ のみ指定（$x > a$ の範囲）。例: 検出限界以上のデータのみ
- 右切断（right truncation）: $b$ のみ指定（$x < b$ の範囲）。例: 寿命が観測期間内のデータのみ
- 両側切断（double truncation）: $a$ と $b$ の両方を指定。例: ある範囲内の測定値のみ
この定義は任意の連続分布に適用できます。最も広く使われるのは正規分布への適用です。次に、切断正規分布を詳しく見ていきましょう。
切断正規分布の定義
直感的な理解
正規分布のベルカーブを思い浮かべてください。左切断正規分布は、このベルカーブの左側を $x = a$ のところでスパッと切り落とし、残った部分の面積が1になるように上に引き伸ばした分布です。
元の正規分布の平均 $\mu$ と分散 $\sigma^2$ は、切断後の分布の平均と分散とは異なります。切断すると左側（平均より小さい部分）が削られるため、切断後の平均は元の平均より大きくなり、分散は小さくなります。
確率密度関数
$X \sim N(\mu, \sigma^2)$ を区間 $(a, b)$ に切断した切断正規分布の確率密度関数は
$$ \begin{equation} f_T(x) = \frac{1}{\sigma} \cdot \frac{\phi\left(\frac{x – \mu}{\sigma}\right)}{\Phi\left(\frac{b – \mu}{\sigma}\right) – \Phi\left(\frac{a – \mu}{\sigma}\right)} \end{equation} $$
ここで $\phi(\cdot)$ は標準正規分布の密度関数、$\Phi(\cdot)$ は標準正規分布の累積分布関数です。
標準化変数 $\alpha = (a – \mu)/\sigma$、$\beta = (b – \mu)/\sigma$ を導入すると、正規化定数は $Z = \Phi(\beta) – \Phi(\alpha)$ と書けます。
左切断のみ（$a$ が有限、$b = +\infty$）の場合は特に単純になります。$\Phi(+\infty) = 1$ なので、
$$ f_T(x) = \frac{\phi((x-\mu)/\sigma)}{\sigma(1 – \Phi(\alpha))}, \quad x > a $$
右切断のみ（$a = -\infty$、$b$ が有限）の場合は $\Phi(-\infty) = 0$ なので、
$$ f_T(x) = \frac{\phi((x-\mu)/\sigma)}{\sigma\Phi(\beta)}, \quad x < b $$
では、この切断によって平均と分散がどのように変化するかを導出しましょう。
切断正規分布のモーメント
平均の導出
切断正規分布の平均を求めるために、$E[X | a < X < b]$ を計算します。
$$ E[X | a < X < b] = \int_a^b x \cdot \frac{\phi((x-\mu)/\sigma)}{\sigma Z} dx $$
$\xi = (x – \mu)/\sigma$ と置換すると、$x = \mu + \sigma\xi$、$dx = \sigma d\xi$ であり、積分範囲は $\alpha$ から $\beta$ になります。
$$ E[X | a < X < b] = \frac{1}{Z}\int_\alpha^\beta (\mu + \sigma\xi)\phi(\xi) d\xi $$
第1項は $\mu \cdot \frac{1}{Z}\int_\alpha^\beta \phi(\xi)d\xi = \mu$（正規化の性質）です。第2項は $\xi\phi(\xi) = -\phi'(\xi)$（標準正規分布の密度の微分の性質）を使うと、
$$ \frac{\sigma}{Z}\int_\alpha^\beta \xi\phi(\xi)d\xi = \frac{\sigma}{Z}[-\phi(\xi)]_\alpha^\beta = \frac{\sigma}{Z}(\phi(\alpha) – \phi(\beta)) $$
したがって、
$$ \begin{equation} E[X | a < X < b] = \mu + \sigma \cdot \frac{\phi(\alpha) - \phi(\beta)}{Z} \end{equation} $$
左切断（$\beta = +\infty$、$\phi(\beta) = 0$）の場合は
$$ E[X | X > a] = \mu + \sigma \cdot \frac{\phi(\alpha)}{1 – \Phi(\alpha)} $$
ここで $\lambda(\alpha) = \phi(\alpha)/(1 – \Phi(\alpha))$ は逆ミルズ比（inverse Mills ratio）と呼ばれ、切断分布の理論で頻繁に登場する重要な量です。逆ミルズ比は常に正なので、左切断すると平均は必ず $\mu$ より大きくなります。これは直感に合います。
分散の導出
分散を求めるには $E[X^2 | a < X < b]$ も必要です。同様の置換を行い、$\xi^2\phi(\xi) = -(\xi\phi(\xi))' + \phi(\xi)$ の関係を使って部分積分すると、
$$ E[X^2 | a < X < b] = \mu^2 + \sigma^2 + \sigma^2 \cdot \frac{\alpha\phi(\alpha) - \beta\phi(\beta)}{Z} + 2\mu\sigma \cdot \frac{\phi(\alpha) - \phi(\beta)}{Z} $$
$V[X | a < X < b] = E[X^2] - (E[X])^2$ を計算して整理すると、
$$ \begin{equation} V[X | a < X < b] = \sigma^2\left[1 + \frac{\alpha\phi(\alpha) - \beta\phi(\beta)}{Z} - \left(\frac{\phi(\alpha) - \phi(\beta)}{Z}\right)^2\right] \end{equation} $$
左切断の場合は $\delta(\alpha) = \lambda(\alpha)(\lambda(\alpha) – \alpha)$ として
$$ V[X | X > a] = \sigma^2(1 – \delta(\alpha)) $$
$0 < \delta(\alpha) < 1$ であることが示せるので、切断正規分布の分散は常に元の正規分布の分散 $\sigma^2$ より小さくなります。これも直感に合っています。分布の片側を切り落とせば、データの散らばりは減少するはずです。
モーメントの計算で逆ミルズ比が鍵を握ることがわかりました。次に、切断と混同されやすい「打ち切り」（censoring）との違いを明確にしておきましょう。
切断と打ち切りの違い
切断（Truncation）
切断では、ある範囲外のデータは存在そのものが観測されません。サンプルサイズ自体が小さくなります。
例: 身長180cm以上の人だけが選ばれるバスケットボール選手のデータ → 180cm未満の人は最初から含まれない
打ち切り（Censoring）
打ち切りでは、ある範囲外のデータは正確な値は不明だが、存在は認識されます。サンプルサイズは変わりません。
例: 測定器の上限が100で、100を超える値はすべて「>100」と記録される → 正確な値は不明だが、データポイントの存在は把握している
統計的な影響
切断と打ち切りでは、尤度関数の構成が異なります。
切断の場合の尤度:
$$ L(\theta) = \prod_{i=1}^{n} \frac{f(x_i | \theta)}{P(a < X < b | \theta)} $$
打ち切り（右打ち切り、$x_i > b$ のとき $b$ と観測される）の場合の尤度:
$$ L(\theta) = \prod_{i: x_i < b} f(x_i | \theta) \cdot \prod_{i: x_i \geq b} (1 - F(b | \theta)) $$
打ち切りでは $x_i \geq b$ のデータ点も $1 – F(b|\theta)$ という確率として尤度に寄与しますが、切断ではそもそも観測されないため、条件付き密度を使います。
この違いを押さえたうえで、切断正規分布の重要な応用であるトービットモデルを見ていきましょう。
トービットモデルへの応用
トービットモデルとは
経済学でよく使われるトービットモデル（Tobit model）は、切断・打ち切りの概念を回帰分析に組み込んだものです。
たとえば、家計の耐久消費財（車や家電）への支出を分析したいとします。多くの家計では支出がゼロですが、支出する家計ではさまざまな金額を支出しています。支出額は負になりえないため、通常の線形回帰は不適切です。
トービットモデルでは、潜在変数 $y_i^*$ を導入します。
$$ y_i^* = \bm{x}_i^T\bm{\beta} + \varepsilon_i, \quad \varepsilon_i \sim N(0, \sigma^2) $$
観測される値は
$$ y_i = \begin{cases} y_i^* & \text{if } y_i^* > 0 \\ 0 & \text{if } y_i^* \leq 0 \end{cases} $$
$y_i > 0$ の観測に対して、$y_i$ の条件付き分布は $\mu_i = \bm{x}_i^T\bm{\beta}$、分散 $\sigma^2$ の正規分布を0で左切断したものになります。
トービットモデルの尤度
対数尤度は
$$ \ell(\bm{\beta}, \sigma^2) = \sum_{i: y_i > 0}\left[\ln\phi\left(\frac{y_i – \bm{x}_i^T\bm{\beta}}{\sigma}\right) – \ln\sigma\right] + \sum_{i: y_i = 0}\ln\Phi\left(-\frac{\bm{x}_i^T\bm{\beta}}{\sigma}\right) $$
第1項は正の観測値からの寄与（通常の正規分布の密度）、第2項はゼロ観測からの寄与（正の観測値が得られない確率）です。
条件付き期待値
$y_i > 0$ が観測されたとき、$y_i$ の条件付き期待値は切断正規分布の平均公式を使って
$$ E[y_i | y_i > 0] = \bm{x}_i^T\bm{\beta} + \sigma \cdot \lambda\left(\frac{\bm{x}_i^T\bm{\beta}}{\sigma}\right) $$
ここで $\lambda(\cdot)$ は逆ミルズ比です。逆ミルズ比の補正項がなければ、OLS回帰の推定値が得られますが、切断の影響を無視するとバイアスが生じます。
ヘックマンの二段階推定法は、この逆ミルズ比を利用して、選択バイアスを補正する手法です。
次に、切断正規分布とトービットモデルをPythonで実装し、その特性を確認しましょう。
Pythonでの実装と可視化
切断正規分布の密度と元の分布との比較
切断正規分布の密度関数が、元の正規分布からどのように変形するかを可視化します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
# 切断正規分布のクラス
class TruncatedNormal:
def __init__(self, mu, sigma, a=-np.inf, b=np.inf):
self.mu = mu
self.sigma = sigma
self.a = a
self.b = b
self.alpha = (a - mu) / sigma if a != -np.inf else -np.inf
self.beta = (b - mu) / sigma if b != np.inf else np.inf
self.Z = stats.norm.cdf(self.beta) - stats.norm.cdf(self.alpha)
def pdf(self, x):
"""確率密度関数"""
xi = (x - self.mu) / self.sigma
density = stats.norm.pdf(xi) / (self.sigma * self.Z)
# 切断範囲外はゼロ
mask = (x >= self.a) & (x <= self.b)
return np.where(mask, density, 0.0)
def mean(self):
"""平均"""
phi_a = stats.norm.pdf(self.alpha) if self.alpha != -np.inf else 0
phi_b = stats.norm.pdf(self.beta) if self.beta != np.inf else 0
return self.mu + self.sigma * (phi_a - phi_b) / self.Z
def var(self):
"""分散"""
phi_a = stats.norm.pdf(self.alpha) if self.alpha != -np.inf else 0
phi_b = stats.norm.pdf(self.beta) if self.beta != np.inf else 0
a_val = self.alpha if self.alpha != -np.inf else 0
b_val = self.beta if self.beta != np.inf else 0
term1 = (a_val * phi_a - b_val * phi_b) / self.Z
term2 = ((phi_a - phi_b) / self.Z) ** 2
return self.sigma**2 * (1 + term1 - term2)
mu, sigma = 0, 1
x = np.linspace(-4, 4, 500)
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# (a) 左切断
ax = axes[0]
normal_pdf = stats.norm.pdf(x, mu, sigma)
ax.plot(x, normal_pdf, 'b--', linewidth=1.5, label='Original $N(0,1)$')
ax.fill_between(x, normal_pdf, alpha=0.1, color='blue')
for a_val in [-1, 0, 1]:
tn = TruncatedNormal(mu, sigma, a=a_val)
ax.plot(x, tn.pdf(x), linewidth=2,
label=f'Left trunc. $a={a_val}$, $E={tn.mean():.2f}$')
ax.set_xlabel('$x$', fontsize=12)
ax.set_ylabel('Density', fontsize=12)
ax.set_title('Left Truncation', fontsize=13)
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
# (b) 右切断
ax = axes[1]
ax.plot(x, normal_pdf, 'b--', linewidth=1.5, label='Original $N(0,1)$')
ax.fill_between(x, normal_pdf, alpha=0.1, color='blue')
for b_val in [-1, 0, 1]:
tn = TruncatedNormal(mu, sigma, b=b_val)
ax.plot(x, tn.pdf(x), linewidth=2,
label=f'Right trunc. $b={b_val}$, $E={tn.mean():.2f}$')
ax.set_xlabel('$x$', fontsize=12)
ax.set_ylabel('Density', fontsize=12)
ax.set_title('Right Truncation', fontsize=13)
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
# (c) 両側切断
ax = axes[2]
ax.plot(x, normal_pdf, 'b--', linewidth=1.5, label='Original $N(0,1)$')
intervals = [(-2, 2), (-1, 1), (-0.5, 0.5)]
for a_val, b_val in intervals:
tn = TruncatedNormal(mu, sigma, a=a_val, b=b_val)
ax.plot(x, tn.pdf(x), linewidth=2,
label=f'$[{a_val},{b_val}]$, $V={tn.var():.3f}$')
ax.set_xlabel('$x$', fontsize=12)
ax.set_ylabel('Density', fontsize=12)
ax.set_title('Double Truncation', fontsize=13)
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('truncated_normal_basic.png', dpi=150, bbox_inches='tight')
plt.show()
上のグラフから、切断の効果が視覚的に理解できます。
- 左図（左切断）: 切断点 $a$ を右に移動させるにつれて、左側が削られて密度関数が右にシフトし、高さが増加します。切断後の平均は $a = -1$ で $E = 0.29$、$a = 0$ で $E = 0.80$、$a = 1$ で $E = 1.53$ と、元の平均0から大きくずれます
- 中央図（右切断）: 右切断は左切断の鏡像です。切断点 $b$ を左に動かすと密度が左にシフトし、平均が負の値になります
- 右図（両側切断）: 区間を狭くするほど密度のピークが高くなり、分散が小さくなります。$[-0.5, 0.5]$ に切断すると分散が0.067まで縮小します
逆ミルズ比の振る舞い
逆ミルズ比は切断正規分布の理論で鍵となる関数です。その振る舞いを可視化します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
alpha_values = np.linspace(-3, 4, 500)
# 逆ミルズ比
phi_a = stats.norm.pdf(alpha_values)
Phi_a = stats.norm.cdf(alpha_values)
mills = phi_a / (1 - Phi_a)
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))
# (a) 逆ミルズ比 λ(α)
ax = axes[0]
ax.plot(alpha_values, mills, 'b-', linewidth=2, label=r'$\lambda(\alpha) = \phi(\alpha)/(1-\Phi(\alpha))$')
ax.plot(alpha_values, alpha_values, 'r--', linewidth=1.5, alpha=0.5, label=r'$y = \alpha$ (asymptote)')
ax.set_xlabel(r'$\alpha = (a - \mu)/\sigma$', fontsize=12)
ax.set_ylabel(r'$\lambda(\alpha)$', fontsize=12)
ax.set_title('Inverse Mills Ratio', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_ylim(-1, 6)
# (b) 切断後の平均と分散の変化
ax2 = axes[1]
# E[X|X>a] = μ + σ λ(α) → 標準化: E = λ(α)
# V[X|X>a] = σ²(1 - δ(α)) where δ = λ(λ-α)
delta = mills * (mills - alpha_values)
mean_shift = mills # (E[X|X>a] - mu) / sigma
ax2.plot(alpha_values, mean_shift, 'b-', linewidth=2, label=r'Mean shift: $\lambda(\alpha)$')
ax2.plot(alpha_values, 1 - delta, 'r-', linewidth=2, label=r'Variance ratio: $1 - \delta(\alpha)$')
ax2.axhline(0, color='gray', linestyle=':', alpha=0.5)
ax2.axhline(1, color='gray', linestyle=':', alpha=0.5)
ax2.set_xlabel(r'$\alpha = (a - \mu)/\sigma$', fontsize=12)
ax2.set_ylabel('Value', fontsize=12)
ax2.set_title('Effect of Left Truncation on Moments', fontsize=13)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.set_ylim(-0.5, 5)
plt.tight_layout()
plt.savefig('mills_ratio.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、逆ミルズ比の重要な性質が読み取れます。
- 左図: 逆ミルズ比 $\lambda(\alpha)$ は $\alpha$ が大きくなる（切断点が平均から右に離れる）と急速に増加し、漸近的に $\alpha$ に近づきます。つまり、切断が厳しくなる（多くの観測が切り落とされる）ほど、平均の補正量が大きくなります
- 右図: 青線は平均のシフト量（$\sigma$ 単位）、赤線は分散比 $V[X|X>a]/\sigma^2$ を示しています。$\alpha$ が増えるほど平均シフトが増大し、分散比は0に向かって減少します。これは切断が厳しくなると、残ったデータが狭い範囲に集中することを反映しています
トービットモデルの推定
最後に、トービットモデルを実装して最尤推定を行い、切断を無視した場合の推定バイアスと比較します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy.optimize import minimize
np.random.seed(42)
# --- データ生成 ---
n = 500
x = np.random.uniform(0, 5, n)
beta_true = np.array([1.0, 0.8]) # 切片, 傾き
sigma_true = 1.5
X = np.column_stack([np.ones(n), x])
y_star = X @ beta_true + np.random.normal(0, sigma_true, n)
# ゼロで切断（観測値）
y = np.maximum(y_star, 0)
is_positive = y > 0
print(f"ゼロ観測の割合: {np.mean(~is_positive)*100:.1f}%")
# --- トービット対数尤度 ---
def tobit_neg_loglik(params, X, y):
beta = params[:-1]
log_sigma = params[-1]
sigma = np.exp(log_sigma)
mu = X @ beta
pos = y > 0
# 正の観測値: 正規分布の対数密度
ll_pos = np.sum(stats.norm.logpdf(y[pos], mu[pos], sigma))
# ゼロ観測: P(y* <= 0) = Φ(-μ/σ)
ll_zero = np.sum(stats.norm.logcdf(-mu[~pos] / sigma))
return -(ll_pos + ll_zero)
# --- 推定 ---
# トービットMLE
init_params = np.array([0.5, 0.5, np.log(1.0)])
result_tobit = minimize(tobit_neg_loglik, init_params, args=(X, y), method='Nelder-Mead')
beta_tobit = result_tobit.x[:-1]
sigma_tobit = np.exp(result_tobit.x[-1])
# ナイーブOLS（正のデータのみ、切断を無視）
from numpy.linalg import lstsq
beta_ols_pos, _, _, _ = lstsq(X[is_positive], y[is_positive], rcond=None)
# 全データOLS（打ち切りを無視）
beta_ols_all, _, _, _ = lstsq(X, y, rcond=None)
print(f"\n推定結果:")
print(f"{'':15s} {'切片':>8s} {'傾き':>8s} {'σ':>8s}")
print(f"{'真の値':15s} {beta_true[0]:8.3f} {beta_true[1]:8.3f} {sigma_true:8.3f}")
print(f"{'Tobit MLE':15s} {beta_tobit[0]:8.3f} {beta_tobit[1]:8.3f} {sigma_tobit:8.3f}")
print(f"{'OLS (y>0 only)':15s} {beta_ols_pos[0]:8.3f} {beta_ols_pos[1]:8.3f} {'---':>8s}")
print(f"{'OLS (all data)':15s} {beta_ols_all[0]:8.3f} {beta_ols_all[1]:8.3f} {'---':>8s}")
# --- 可視化 ---
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# (a) データと回帰直線
ax = axes[0]
ax.scatter(x[is_positive], y[is_positive], s=15, alpha=0.5, color='blue', label='$y > 0$')
ax.scatter(x[~is_positive], y[~is_positive], s=15, alpha=0.5, color='red', label='$y = 0$')
x_grid = np.linspace(0, 5, 100)
ax.plot(x_grid, beta_true[0] + beta_true[1]*x_grid, 'k-', linewidth=2,
label=f'True: $y = {beta_true[0]:.1f} + {beta_true[1]:.1f}x$')
ax.plot(x_grid, beta_tobit[0] + beta_tobit[1]*x_grid, 'g--', linewidth=2,
label=f'Tobit: $y = {beta_tobit[0]:.2f} + {beta_tobit[1]:.2f}x$')
ax.plot(x_grid, beta_ols_pos[0] + beta_ols_pos[1]*x_grid, 'r:', linewidth=2,
label=f'OLS(y>0): $y = {beta_ols_pos[0]:.2f} + {beta_ols_pos[1]:.2f}x$')
ax.axhline(0, color='gray', linestyle='-', alpha=0.3)
ax.set_xlabel('$x$', fontsize=12)
ax.set_ylabel('$y$', fontsize=12)
ax.set_title('Tobit Model vs Naive OLS', fontsize=13)
ax.legend(fontsize=8, loc='upper left')
ax.grid(True, alpha=0.3)
# (b) 条件付き期待値 E[y|y>0]
ax = axes[1]
mu_grid = beta_true[0] + beta_true[1] * x_grid
alpha_grid = -mu_grid / sigma_true
mills = stats.norm.pdf(alpha_grid) / (1 - stats.norm.cdf(alpha_grid))
ey_given_pos = mu_grid + sigma_true * mills
ax.plot(x_grid, mu_grid, 'b-', linewidth=2, label='$E[y^*] = X\\beta$')
ax.plot(x_grid, ey_given_pos, 'r-', linewidth=2, label='$E[y|y>0]$')
ax.fill_between(x_grid, mu_grid, ey_given_pos, alpha=0.15, color='orange',
label='Mills ratio correction')
ax.set_xlabel('$x$', fontsize=12)
ax.set_ylabel('Expected value', fontsize=12)
ax.set_title('Conditional Expectation $E[y|y>0]$', fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# (c) モンテカルロシミュレーション: バイアスの比較
ax = axes[2]
n_sim = 500
n_obs = 200
biases_tobit_beta1 = []
biases_ols_pos_beta1 = []
biases_ols_all_beta1 = []
for _ in range(n_sim):
x_sim = np.random.uniform(0, 5, n_obs)
X_sim = np.column_stack([np.ones(n_obs), x_sim])
y_star_sim = X_sim @ beta_true + np.random.normal(0, sigma_true, n_obs)
y_sim = np.maximum(y_star_sim, 0)
pos_sim = y_sim > 0
# Tobit
try:
res = minimize(tobit_neg_loglik, init_params, args=(X_sim, y_sim),
method='Nelder-Mead', options={'maxiter': 2000})
biases_tobit_beta1.append(res.x[1] - beta_true[1])
except:
pass
# OLS (y>0 only)
if np.sum(pos_sim) > 2:
b, _, _, _ = lstsq(X_sim[pos_sim], y_sim[pos_sim], rcond=None)
biases_ols_pos_beta1.append(b[1] - beta_true[1])
# OLS (all)
b, _, _, _ = lstsq(X_sim, y_sim, rcond=None)
biases_ols_all_beta1.append(b[1] - beta_true[1])
ax.boxplot([biases_tobit_beta1, biases_ols_pos_beta1, biases_ols_all_beta1],
labels=['Tobit MLE', 'OLS (y>0)', 'OLS (all)'],
patch_artist=True,
boxprops=[dict(facecolor='lightgreen'),
dict(facecolor='lightsalmon'),
dict(facecolor='lightblue')][0:1]*3)
# 個別に色を設定
bp = ax.boxplot([biases_tobit_beta1, biases_ols_pos_beta1, biases_ols_all_beta1],
labels=['Tobit MLE', 'OLS (y>0)', 'OLS (all)'],
patch_artist=True)
colors = ['lightgreen', 'lightsalmon', 'lightblue']
for patch, color in zip(bp['boxes'], colors):
patch.set_facecolor(color)
ax.axhline(0, color='red', linestyle='--', linewidth=1.5, alpha=0.7)
ax.set_ylabel(r'Bias in $\hat{\beta}_1$', fontsize=12)
ax.set_title(f'Estimation Bias ({n_sim} simulations)', fontsize=13)
ax.grid(True, alpha=0.3, axis='y')
plt.tight_layout()
plt.savefig('tobit_model.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、トービットモデルの有効性が確認できます。
-
左図（データと回帰直線）: 真の回帰直線（黒実線）に対して、トービットMLE（緑破線）はほぼ一致していますが、ナイーブOLS（赤点線、$y > 0$ のデータのみ）は切片が過大推定され、傾きがバイアスを持っています。切断を無視すると、切断点近傍のデータが過少に代表されるためです
-
中央図（条件付き期待値）: 潜在変数の期待値 $E[y^*] = X\beta$（青）と正の観測の条件付き期待値 $E[y|y>0]$（赤）の間には逆ミルズ比による差（オレンジ領域）があります。$x$ が小さい（$\mu$ が0に近い）ときほど差が大きく、切断の影響が強いことがわかります
-
右図（バイアスのシミュレーション）: 500回のシミュレーションにおける推定バイアスの分布を比較すると、トービットMLEはゼロを中心に対称に分布しており、不偏に近い推定を実現しています。一方、OLS（$y > 0$ のみ）とOLS（全データ）はともに系統的なバイアスを示しています
まとめ
本記事では、切断分布の定義から切断正規分布の性質、トービットモデルへの応用までを解説しました。
- 切断分布は、元の分布を特定の範囲に制限し、確率が1になるように再正規化した分布である
- 切断正規分布の平均と分散は、逆ミルズ比 $\lambda(\alpha) = \phi(\alpha)/(1-\Phi(\alpha))$ を通じて簡潔に表現できる
- 切断は平均をシフトさせ、分散を縮小させる。この効果は切断が厳しいほど顕著になる
- 切断と打ち切りは異なる概念であり、尤度関数の構成が異なる
- トービットモデルは切断の概念を回帰分析に組み込み、選択バイアスを適切に補正する
- 切断を無視したナイーブな推定は系統的なバイアスを生じる
次のステップとして、以下の記事も参考にしてください。