カメラが25%しか見えなくても全体の数がわかる：N-mixtureモデルで個体数を推定するPythonシミュレータ

カメラに映らない個体も含めて「全体の数」を推定できるか？

AIカメラによる人数カウント・動物の生態調査・工場の生産ライン管理。こうしたシステムでは、「カメラの視野に映っている数」は測れても、「空間全体にいる数」はわからない という問題が常につきまといます。

• カメラの画角外にいる個体は？
• AIが検知し損ねた個体は？

この問題に対して、生態学が長年研究してきた強力なツールが N-mixture モデル です。

「繰り返し観測したカウントデータの平均と分散だけから、真の総個体数 N と検知確率 p を同時推定する」

ここで重要なのは、カメラの視野面積比（何%をカバーしているか）は推定の入力として使わないということです。面積比を事前に知らなくても、カウントデータだけから N と p を推定できます。推定後に p が面積比と一致していれば、それはモデルが正しく動いていることの検証になります。

今回は、このN-mixtureモデルをPythonで実装し、「本当に推定できるか」を視覚的に検証するシミュレータを作りました。

システム全体の構成

シミュレータは3つのモジュールで構成されています。

① 世界モジュール：ランダムウォークする個体群

設定

各個体はランダムな方向に一定速度で移動し、壁に当たると反射します（完全弾性反射）。壁での反射は速度の該当成分を反転させるだけです。

def step(self):
    self.positions += self.velocities
    for dim in range(2):
        low  = self.positions[:, dim] < 0.0
        high = self.positions[:, dim] > self.space_size
        # 壁を超えた分だけ折り返す
        self.positions[low,  dim] = -self.positions[low,  dim]
        self.positions[high, dim] = 2.0 * self.space_size - self.positions[high, dim]
        # 速度の向きを反転
        self.velocities[low,  dim] =  abs(self.velocities[low,  dim])
        self.velocities[high, dim] = -abs(self.velocities[high, dim])

撮影間の独立性の確保

N-mixtureモデルは、各撮影が互いに独立であることを前提とします。つまり「1回目に75匹映った」という結果が「2回目の結果に影響を与えない」必要があります。

しかしランダムウォークでは、直前のフレームで視野内にいた個体はほぼその場にとどまっているため、連続した2枚の写真はほぼ同じカウントになってしまいます。これでは「500回の独立した観測」ではなく「ほぼ同じ写真を500枚撮っただけ」になります。

そこでシャッターを切るたびに全個体の進行方向だけをランダムにリセットします。位置はそのままなのでアニメーションは自然に見えますが、統計的には「毎回バラバラな状態から出発した独立な観測」として扱えるようになります。

def randomize_directions(self):
    """速度方向をランダムリセット（位置は維持）"""
    angles = np.random.uniform(0.0, 2.0 * np.pi, self.n)
    self.velocities = self.speed * np.column_stack([np.cos(angles), np.sin(angles)])

100試行のシミュレーションで、この処置により推定精度の中央値誤差が 6.5% → 3.2% に改善されることを確認しています。

② 観測モジュール：AIカメラ

空間中央に 50 × 50 の視野（FOV: Field of View）を設定します。視野面積を全体面積で割ったものが面積比率 p_{\text{area}} （p sub area）です：

視野内にいる全個体を100%検知します（第一段階の理想カメラ仮定）。

def observe(self, positions):
    in_fov = (
        (positions[:, 0] >= self.x0) & (positions[:, 0] <= self.x1) &
        (positions[:, 1] >= self.y0) & (positions[:, 1] <= self.y1)
    )
    return int(np.sum(in_fov)), in_fov

これを30フレームおきに500回繰り返し、カウントデータ y = [y_1, y_2, \ldots, y_{500}] を収集します。

③ 推論モジュール：N-mixture モデル（最尤推定）

モデルの前提

「N 匹いて、そのうち p の割合が視野に入る」という状況は、コインを N 回投げて表が出た回数を数えるのと同じ構造です。これを二項分布と呼びます。

• N ：空間全体の真の総個体数（未知・推定したい値）
• p ：1個体が視野内に入る確率（未知・推定したい値）
• y_i ：第 i 回撮影で実際に映った数

期待値と分散から N と p を取り出す

二項分布には便利な性質があります。「平均（ E ）」と「分散（ Var ）」を計算すると、それぞれ N と p を含む式になります：

つまり「平均 E」と「分散 Var」という2つの数字から、2つの未知数 N と p を連立方程式で解けます。これがこのモデルの核心です。

ステップ1：まず大まかに計算する（モーメント法）

500回の撮影データから標本平均 \bar{y} と標本分散 s^2 を計算し、上の2式に当てはめます。

ここで「分散 ÷ 平均」を計算すると、両方に含まれる N が約分されて消え、 p だけの式が残ります：

分子の Np(1-p) を分母の Np で割ると、 N と p が約分されて (1-p) だけ残ります。これを変形すると p が求まります。

次に求めた p を平均の式 E[y] = Np に代入すれば N も求まります：

具体的な数値で確認してみます。 今回のシミュレーション結果（ \bar{y} = 75.26 、 s^2 = 56.27 ）を当てはめると：

真の値（ p = 0.25 、 N = 300 ）に近い値が出ていることがわかります。

なぜこれで推定できるのかを直感的に：

• カウントが毎回ほぼ同じ数（分散が小さい）→ 撮るたびに安定して映っている → 視野が全体をよくカバーしている（ p が大きい）→ 視野外にいる個体は少ない → N は小さい
• カウントが毎回バラバラ（分散が大きい）→ 映ったり映らなかったり不安定 → 個体が広い空間をランダムに動き回っている → 視野の外に大量にいる → N は大きい

つまり「毎回のカウントのばらつき具合」が、視野外の個体数を推測する唯一の手がかりになっています。

ステップ2：より正確に絞り込む（最尤推定）

ステップ1の値はあくまで粗い初期値です。「この N と p なら、実際のカウントデータが得られる確率はどのくらいか」を計算し、その確率を最も高くする N を探します。これが最尤推定（MLE） です。

「この N と p のとき、観測されたカウントデータが得られる確率（尤度）」を全撮影回分だけ掛け合わせ、それを最大化する N を探します。掛け算のままだと数値が極端に小さくなるため、対数をとって足し算に変換したものが対数尤度 \ell です：

式の各項の意味：

• \log \binom{N}{y_i} ：N 匹の中から y_i 匹を選ぶ組み合わせの数（何通りの出方があるか）
• y_i \log p ：実際に映った y_i 匹が「視野内に入った」確率
• (N - y_i) \log(1-p) ：残りの N - y_i 匹が「視野外だった」確率

しかし未知数が N と p の2つあるため、このままでは探索が2次元になってしまいます。 そこで「 N が決まれば p の最適値も自動的に決まる」という性質を使います。対数尤度を p で微分してゼロとおくと：

つまり「カウントの平均 ÷ 候補の N 」が、その N に対して最も尤もらしい p です。これを代入して p を消去した式がプロファイル対数尤度です：

あとは グリッド探索 で \hat{N} を求めます。

from scipy.special import gammaln

def nmixture_mle(counts):
    y = np.array(counts, dtype=float)
    y_mean = np.mean(y)
    y_var  = np.var(y, ddof=1)

    # モーメント法による初期推定（グリッド探索の中心に使用）
    if y_mean > y_var > 0:
        p_mom = 1.0 - y_var / y_mean
        N_mom = y_mean / p_mom
    else:
        # 分散≥平均の場合（過分散）はフォールバック
        N_mom = y_mean / AREA_RATIO

    # プロファイル対数尤度のグリッド探索
    best_ll = -np.inf
    best_N  = int(round(N_mom))
    for N_cand in range(int(N_mom * 0.4), int(N_mom * 3.0) + 100):
        p_cand = y_mean / N_cand
        if not (0.0 < p_cand < 1.0):
            continue
        log_binom = (gammaln(N_cand + 1)
                     - gammaln(y + 1)
                     - gammaln(N_cand - y + 1))
        ll = float(np.sum(log_binom
                          + y * np.log(p_cand)
                          + (N_cand - y) * np.log(1.0 - p_cand)))
        if ll > best_ll:
            best_ll = ll
            best_N  = N_cand

    return float(best_N), y_mean / best_N

gammaln（対数ガンマ関数）を使って \log \binom{N}{y_i} を計算しています。直接 N! を計算するとオーバーフローするため、この形が標準的です。

シミュレーション結果

300個の個体が100×100の空間を動き回り、中央のカメラが500回シャッターを切る様子です。右のグラフにカウントデータが蓄積され、最後に N と p が表示されます。

実行結果：

====================================================
  【答え合わせ / N-mixture 推定結果】
====================================================
  真の総数         : 300 個
  推定総数 (N_hat) : 298 個  （誤差 0.7%）
  カメラ視野面積比 : 0.2500
  推定発見率 (p)   : 0.2526  （誤差 1.0%）
----------------------------------------------------
  観測回数         : 500
  カウント平均     : 75.26
  カウント分散     : 56.27
====================================================

「面積比を入力していないのに p ≈ 0.25 が出る」のはなぜか

出力をよく見ると、推定発見率 \hat{p} がカメラ視野面積比 0.25 にほぼ一致しています。

  カメラ視野面積比 : 0.2500
  推定発見率 (p)   : 0.2526

「結局 50 \times 50 / 100 \times 100 = 0.25 を使ってるだけでは？」と思うかもしれませんが、推論モジュールはこの面積比を一切入力として使っていません。nmixture_mle() が受け取るのはカウントデータ y のみです。

モデルが行っているのは、カウントデータの平均と分散の比だけから p を導出することです：

\hat{p} \approx 0.25 が面積比と一致するのは、モデルが「個体が空間に一様分布しているとき、視野内に入る確率は面積比になる」という物理的事実を、データから独立に復元できていることの検証です。

これは実運用で重要な意味を持ちます。N-mixtureモデルが事前に必要とする情報はカウントデータだけです。

統計的な注意点と考察

注意点1：この推定器は「重尾分布」を持つ

ステップ1で求めた推定式（ \hat{N} = 推定総個体数）を思い出してください：

分母の \bar{y} - s^2 （平均 − 分散）がゼロに近づくと、 \hat{N} が無限大に発散します。標本分散 s^2 はサンプルごとにランダムにばらつくため、まれに分母がゼロ付近になり推定値が極端に大きくなることがあります。

注意点2：観測数 n が重要

推定の核心は「カウントの分散と平均の比」を精度よく測ることです。この比の推定誤差（SE）は \propto 1/\sqrt{n} で減少するため、観測回数 n が多いほど推定が安定します。これを信号（推定したい量 p ）とノイズ（分散推定の誤差）の比として整理すると、推定精度のシグナル対ノイズ比（SNR）は：

p=0.25 （25%の視野面積比）で N と p を同時推定するには、 n ≥ 200 が実用上の目安で、本実装では n=500 を採用しています。

注意点3：過分散（Overdispersion）問題

二項分布の理論では \text{Var} < \text{Mean} （ \text{Var/Mean} = 1 - p < 1 ）が成り立つ必要があります。

しかし連続したランダムウォークでは、隣接する観測の位置が相関しているため、標本分散が理論値より大きくなることがあります（過分散）。この場合：

• s^2 > \bar{y} → \hat{p} = 1 - s^2/\bar{y} がマイナスになる（確率がマイナスは不正な値）
• MLEの尤度が単調増加 → 推定値が上限まで飛ぶ

対処法として実装していること：

1. 撮影ごとに速度方向をランダムリセット（相関を低減）
2. 過分散検出時は面積比をフォールバックとして使用
3. 観測数を500回確保

実際の生態学調査での対応： シミュレータでは「速度ランダムリセット」というハックで独立性を確保していますが、現実の野外調査では「時間を十分に空けて繰り返し観測する」か、「複数のカメラを異なる場所に設置して空間的な繰り返しとする」アプローチが一般的です。後者は個体の移動に依存しないため、より堅牢に独立性を担保できます。

注意点4：モデルの仮定（第一段階の理想化）

今回の実装は 第一段階として以下を仮定しています：

課題・改善の方向性

課題1：不完全検知の導入（第二段階）

現在は「視野内に入れば必ず検知できる」理想的なカメラを仮定しています。実際のAIモデルには見逃し（検知漏れ）があります。

これを取り込むと、カウント y_i （第 i 回の撮影で映った数）のモデルは次のようになります：

• N_{\text{FOV},i} （N sub FOV）：第 i 回撮影時に視野内にいた真の個体数
• p_{\text{detect}} （p sub detect）：AIが視野内の個体を実際に検知できる確率

この場合、推定すべきパラメータは 3つ に増えます：

• N ：空間全体の総個体数
• p_{\text{area}} ：カメラが空間をカバーしている面積の割合（現在は既知の0.25）
• p_{\text{detect}} ：AIの検知確率

パラメータが増える分、より多くの観測回数が必要になります。

課題2：空間的不均一性

現在の実装では個体が空間に一様に分布すると仮定しています。しかし実際には個体が集まる「クラスタ」を形成したり、特定のエリアを好むことがあります。

こうした「ばらつきのばらつき」（過分散）に対応するには、Negative Binomial分布（負の二項分布）や Beta-Binomial分布 を使ったN-mixtureの拡張版が有効です。

• \mu （ミュー）：個体数の平均
• \theta （シータ）：ばらつきの大きさを表すパラメータ（小さいほど過分散が大きい）

課題3：開放個体群への拡張

今回は「調査中に個体が出入りしない」閉鎖個体群を仮定しています。実際のシステムでは個体が領域に出入りします。この場合は 動的N-mixtureモデル（Dail-Madsen model） への拡張が必要です。

課題4：MLEの計算効率

現在のグリッド探索は候補の N を1つずつ全部試しているため、探索範囲に比例した計算量（ O(N_{\max}) ）がかかります。 N_{\max} が大きくなると処理が重くなります。

大規模な問題では MCMC（マルコフ連鎖モンテカルロ法）や変分推論を使ったベイズ推定が現実的です。特に pymc や numpyro を使えばN-mixtureを数行で記述できます：

import pymc as pm

with pm.Model():
    N = pm.DiscreteUniform("N", lower=y_max, upper=2000)
    p = pm.Uniform("p", 0, 1)
    y_obs = pm.Binomial("y_obs", n=N, p=p, observed=counts)
    trace = pm.sample(2000)

実装全体

シミュレータは以下の構成です：

nmixture_simulator.py
├── class World          # ランダムウォーク物理シミュレーション
├── class Camera         # FOV内カウント観測
├── def nmixture_mle()   # N-mixture最尤推定
└── def run_simulation() # アニメーション・UI・結果表示

ソースコード：

GitHub - ramtuc/ram_prog / N-mixture

実行環境：

pip install numpy matplotlib scipy
python nmixture_simulator.py

まとめ

「カメラが空間の25%しか見えない」状況でも、繰り返し観測のカウントデータから統計モデルが真の個体数を推定できることを確認できました。

N-mixtureモデルは生態学では定番の手法ですが、AIカメラによる混雑状況モニタリング・工場の生産管理・セキュリティシステムなど、工学的な応用にも十分使える考え方です。「見えていない部分を推論する」という発想が、センサーデータの分析に新しい視点をもたらしてくれます。

コードは GitHub で公開しています。pip install numpy matplotlib scipy だけで動くので、ぜひ手元で試してみてください。

さらに学ぶには

• Royle, J.A. (2004) - N-Mixture Models for Estimating Population Size from Spatially Replicated Counts. Biometrics, 60(1), 108-115. — N-mixtureモデルの原論文
• Kéry & Schaub (2012) - Bayesian Population Analysis using WinBUGS — 実践的なベイズ個体数推定の教科書
• pymc ドキュメント / numpyro ドキュメント — PythonでのベイズN-mixture実装

関連記事

• Arduinoと2次元セル・オートマトン：64ドットに宿る生命シミュレーション
• Azure Functions + ESP32 でIoTデータを送信する