【改】第５章「樹木の直径と高さとの関係のモデリング」のベイズモデリングをPyMC Ver.5 で

この記事は、以前投稿した『 第５章「樹木の直径と高さとの関係のモデリング」のPyMC記事 』を改造する記事です。

以前の記事はこちらです。

以前の記事執筆時点では、事後分布を収束させることができず、苦い思いをしました。
しかし、事態は一変します。
生成AIの進化のおかげで、解決策が見つかったのです！

それではテイクツー、スタートです！

リスタート

この記事は、テキスト「たのしいベイズモデリング２」の第５章「樹木の直径と高さとの関係のモデリング」のベイズモデルを用いて、PyMC Ver.5で「実験的」に実装する様子を描いた統計ドキュメンタリーです。

この章では、京都市の天然林の樹木データを用いて、直径と高さの関係をモデリングします。
直径と高さの関係には「拡張相対成長式」が用いられます。
このモデルがあれば、幹の直径の観測値で木の高さを推論できるのです！

しかし今回は、とても複雑なモデル。
前回記事では無事に収束を迎えることができず・・・
今回の再挑戦でテキストのベイズモデルのお気持ちに近づけました！
ではでは一緒に PyMCのベイズモデリングの世界 を楽しみましょう。

テキストの紹介、引用表記、シリーズまえがき、PyMC等のバージョン情報は、このリンクの記事をご参照ください。

テキストで使用するデータは、R・Stan等のサンプルスクリプトとともに、出版社サイトからダウンロードして取得できます。

サマリー

テキストの概要

執筆者　　　：　伊東宏樹 先生
モデル難易度：　★★★・・　（ふつう）

自己評価

評点

$$
\begin{array}{c:c:c}
実装精度 & ★★★★★ & GoooD! \\
結果再現度 & ★★★★★ & 最高！ \\
楽しさ & ★★★★★ & 楽しい！ \\
\end{array}
$$

評価ポイント

• 生成AIのClaudeと一緒に前回モデルの改善点を検討して、以下の２つの対策を実施しました。

  • Divergences 回避：非中心化パラメータ化

  • モデルのバグ修正

工夫・喜び・反省

• ギリギリ収束（まじギリギリ）ですが…
  テキストの推論結果に近づけたと思います！よかったです。

モデルの概要

テキストの調査・実験の概要

■ 樹木の直径と高さのモデリングって？
木の高さを測定するのは、かなり手間がかかるそうです。
そこで、一部の木の直径と高さを測ってモデリングし、残りの木は直径だけを測定、モデルを用いて高さを推定するのです。
直径の測定については、一定の高さで幹の周囲を測って円周率で割って直径値とする方法があり、この方法は測定の手間がかからないというわけです。

■ 相対成長
背の高い人は腕も長い、というような関係を相対成長関係と呼ぶそうです。
直径の大きい木は高い、という関係に通じるのでしょうか？
樹木の場合、直径の大きさと木の高さ「樹高」の関係について、次式の「拡張相対式」がよく使われるそうです。

$${H}$$は樹高、$${H_D}$$は直径$${D}$$を測定した高さ、$${H_{\text{Max}}}$$は潜在的に取り得る最大の樹高、$${a,h}$$は曲線の形状を決めるパラメータです。
ベイズモデリングでは、下式の両辺の対数を取った数式が用いられます。

【お知らせ】
今回は図示を省略いたします。

■ データの概要
京都市の銀閣寺山間国有林の天然林（二次林）で測定した15種類、232本の木の直径と樹高のデータです。
この国有林は五山送り火で知られる大文字山のふもとにあるそうです。
データは次のサイトで公開されています。
「補足電子資料：表S1～S3（エクセル2007：30KB）」

テキストのモデリング

■目的変数と関心のあるパラメータ
目的変数は、樹高の観測値$${H}$$から直径を測定した高さの観測値$${H_D}$$を差し引いた値「$${H_{ij}-H_D}$$」です。
添字$${i}$$は木のID、$${j}$$は樹高測定回です。
また、樹高と関係のある直径の観測値$${D}$$も目的変数に用います。

関心のあるパラメータに関してはテキストに特に明記されていませんが、おそらく、全パラメータに注目します。
直径から樹高を推定するモデルの構築には、全パラメータが必要だからです。

■ テキストのベイズモデル
テキストの数式、文章補足からモデル数式を書きました。
数式の具体的な内容はぜひ、テキストでお確かめください。

主要な変数の内容は次のとおりです。

■分析・分析結果
分析の仕方や分析数値はテキストの記述が正確だと思いますので、テキストの読み込みをおすすめします。
PyMCの自己流モデルの推論値を利用した分析は「PyMC実装」章をご覧ください。

PyMC実装

Let's enjoy PyMC & Python !

準備・データ確認

１．インポート

### インポート

# 数値・確率計算
import pandas as pd
import numpy as np

# PyMC
import pymc as pm
import pytensor.tensor as pt
import arviz as az

# WEBアクセス
import requests

# 描画
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'

# ワーニング表示の抑制
import warnings
warnings.simplefilter('ignore')

２．データの読み込み
WebサイトからEXCELファイルをダウンロードして分析に用います。
まずダウンロードしましょう。
「path」に「Excelファイルの格納フォルダ＋ファイル名」を指定します。

# ### EXCELファイルのダウンロード

## 初期値設定
# EXCEL（.xlsx）ファイルのURL
url = 'https://www.ffpri.go.jp/pubs/bulletin/435/documents/tables-s1-s3.xlsx'
# EXCELファイルの格納フォルダ＋ファイル名
path = 'tables-s1-s3.xlsx'

## ダウンロード
# Web接続
res = requests.get(url)
# ダウンロードファイルの保存
with open(path, mode='wb') as f:
    f.write(res.content)

【実行結果】なし

続いて、EXCELファイルをpandasのデータフレームに読み込みます。

### ダウンロードしたEXCELファイルの読み込み

# 初期値設定：列名
cols = ['和名', '学名', 'No', '胸高直径', '樹高1', '樹高2', '樹高3', '樹高4']
# EXCELファイルの読み込み
data_orgn = pd.read_excel(path, sheet_name=0, header=None, skiprows=2, names=cols)
# データの表示
display(data_orgn)

【実行結果】
Noは木の個体IDです。
胸高直径の単位はｃｍです。
樹高は測定回１～測定回４まで横持ちしています。単位はｍです。
なお、樹高の測定回数は個体によってまちまちです。
例えば１行目のコシアブラNo.1は３回測定、３行目のコシアブラNo.3は２回測定です。

Excelファイルをダウンロードしない（フォルダに保存しない）で、URLから直接データ参照する場合は、次のコードを実行します。

### 参考：データの読み込み（WEBサイトから直接読み込み）

## 初期値設定
# EXCEL（.xlsx）ファイルのURL
url = 'https://www.ffpri.go.jp/pubs/bulletin/435/documents/tables-s1-s3.xlsx'
# 列名
cols = ['和名', '学名', 'No', '胸高直径', '樹高1', '樹高2', '樹高3', '樹高4']

## データの読み込み
# データの読み込み
data_orgn = pd.read_excel(url, sheet_name=0, header=None, skiprows=2, names=cols)
# データの表示
display(data_orgn)

# データをcsvファイル形式で保存する場合のコード
# data_orgn.to_csv('data.csv', encoding='utf_8_sig', index=False)
# 参考：保存したcsvファイルの読み込み
# data_orgn = pd.read_csv('data.csv')

３．データの外観の確認
樹種（木の種類）ごとにデータの要約統計量を表示します。
テキストの表5.1に相当します。

### データに含まれる木の種類とタイプ・平均直径・平均樹高 ※表5.1に相当

## 初期値設定 樹のタイプ
# 樹のタイプ辞書、各樹のタイプリストの作成
type_dic = {1: '落葉広葉樹', 2: '常緑広葉樹', 3: '常緑針葉樹'}
tree_type = [1, 2, 1, 2, 2, 1, 1, 2, 3, 2, 1, 1, 2, 2, 1]

## 平均計算をしながらデータフレームを作成
# 測定本数・平均直径の算出 ※groupbyで和名単位の胸高直径のcountとmeanを計算
summary_df =  (data_orgn.groupby(['和名'])['胸高直径'].agg(['count', 'mean'])
               .set_axis(['測定本数', '平均直径(cm)'], axis=1))
# 平均樹高の算出 ※meltで縦持ち変換した後にgroupbyで和名単位の平均樹高を計算
tmp = (data_orgn[['和名', '樹高1', '樹高2', '樹高3', '樹高4']]
       .melt(id_vars=['和名'], value_vars=['樹高1', '樹高2', '樹高3', '樹高4'],
             var_name='測定回', value_name='平均樹高(m)')
       .dropna()
       .groupby(['和名'])[['平均樹高(m)']].mean())
# 各樹のタイプ（日本語）のデータフレーム化
tree_type_df = pd.DataFrame([type_dic[x] for x in tree_type],
                            index=summary_df.index, columns=['タイプ'])
# ３つのデータフレームを結合してデータフレームを最終化
summary_df = pd.concat([tree_type_df, summary_df, tmp], axis=1)

## データフレームの表示
display(summary_df.round(1))

【実行結果】
各樹種ごとの測定本数は12～17本、平均直径・平均樹高はさまざまです。
平均樹高の値の一部はテキストと異なりますが、原因は不明です。

この平均値データを用いて、木の大きさを可視化しましょう。
棒グラフを描画します。

### 平均直径と平均樹高を棒グラフで描画

# 平均直径と平均樹高を昇順で並び替え
d_sorted = summary_df['平均直径(cm)'].sort_values()
h_sorted = summary_df['平均樹高(m)'].sort_values()

# 描画
fig, ax = plt.subplots(1, 2, figsize=(10, 4), tight_layout=True)
# 平均直径の棒グラフの描画
ax[0].bar(d_sorted.index, d_sorted)
ax[0].xaxis.set_tick_params(rotation=60)
ax[0].set(title='平均直径', ylabel='平均直径(cm)')
# 平均樹高の棒グラフの描画
ax[1].bar(h_sorted.index, h_sorted)
ax[1].xaxis.set_tick_params(rotation=60)
ax[1].set(title='平均樹高', ylabel='平均樹高(m)')
plt.show()

【実行結果】
小さな値から並び替えました。
平均直径・平均樹高の並び順は同じです。
樹種間を比べると「直径が大きくなると樹高が高くなる」といった関係があるように見えます。

４．データの前処理
データ横持ちの測定回数ごとの樹高をモデルで扱いやすいように縦持ちに変換して、樹高データとします。
胸高直径データは１回しか測定されていないので、樹高データから独立したデータにします。
さらに、樹種の和名をラベルにしてインデックスを作成します。

### データの前処理

##  data_orgnのインデックスを付加した一時データの作成
tmp = data_orgn.copy().reset_index(drop=False)

## 樹高データの作成
data_H = (tmp.melt(id_vars=['index', '和名', 'No', '胸高直径'],
                   value_vars=['樹高1', '樹高2', '樹高3', '樹高4'],
                   var_name='測定回', value_name='樹高')
          .dropna().sort_values(['index', '和名', 'No', '測定回'])
          .reset_index(drop=True))
print('【 樹高データ H 】')
display(data_H)

## 直径データの作成
data_D = tmp[['index', '和名', 'No', '胸高直径']].reset_index(drop=True)
print('【 直径データ D 】')
display(data_D)

## 樹高データの樹種インデックスの作成
species_idx, species_labels = pd.factorize(data_H['和名'], sort=False)
print('【 樹種 ラベルとインデックス】')
print(species_labels)
print(species_idx)

【実行結果】
こちらは前処理後の樹高データです。

こちらは前処理後の直径データです。

こちらは樹種のラベルとインデックスです。

モデル構築

モデルの数式表現
目指したいPyMCのモデルの雰囲気を混ぜた「なんちゃって数式」表記です。

$$
\begin{align*}
\sigma_{D},\sigma_{H}, \sigma_{\epsilon}, \sigma_{\epsilon_a},\sigma_{\epsilon_h},\sigma_{\epsilon_{H_{max}}} &\sim \text{HalfNormal}\ (\text{sigma}=10) \\ 
 \\
\log \overline{a} &\sim \text{Normal}\ (\text{mu}=0,\ \text{sigma}=10)\\
\epsilon_a &\sim \text{Normal}\ (\text{mu}=0, \ \text{sigma}=\sigma_{\epsilon_a}) ★ \\
a &\sim \exp(\log \overline{a} + \epsilon_a)\\
 \\
\log \overline{h} &\sim \text{Normal}\ (\text{mu}=0,\ \text{sigma}=10)\\
\epsilon_h &\sim \text{Normal}\ (\text{mu}=0,\ \text{sigma}=\sigma_{\epsilon_h}) ★ \\
h &\sim \exp(\log \overline{h} + \epsilon_h)\\
 \\
\log \overline{H}_{max} &\sim \text{Normal}\ (\text{mu}=4,\ \text{sigma}=4)\\
\epsilon_{H_{max}} &\sim \text{Normal}\ (\text{mu}=0,\ \text{sigma}=\sigma_{\epsilon_{H_{max}}}) ★ \\
H_{max} &\sim \exp(\log \overline{H}_{max} + \epsilon_{H_{max}})\\
 \\
\overline{D}  &\sim \text{HalfNormal}\ (\text{sigma}=100) \\
\epsilon &\sim \text{Normal}\ (\text{mu}=0,\ \text{sigma}=\sigma_{\epsilon},\\ &\quad\quad\quad\quad\quad \textbf{dims} \bm{=} \textbf{data}_{\bm D}) ★ \\
 \\
\log \overline{H} &= \log a + h \times \log \overline{D} + \log H_{max} \\
&\quad - \log (a \times \overline{D}^h + H_{max}) + \epsilon \\
 \\
likelihoodD &\sim \text{Normal}\ (\text{mu}=\overline{D},\ \text{sigma}=\sigma[D])\\
likelihoodH2 &\sim \text{LogNormal}\ (\text{mu}=\log\overline{H},\ \text{sigma}=\sigma[H])\\
\end{align*}
$$

前回記事からのモデル修正点は以下のとおりです。

1. 非中心化パラメータ化
   ★印を非中心化パラメータ化します。Divergences対策になります。

2. バグ修正
   パラメータ $${\epsilon }$$について配列の次元が誤っていました。
   太字のとおり「$${\text{data}_D}$$」（木の個体単位）にします。

１．モデルの定義
数式表現を実直にモデル記述します。

### モデルの定義

# 初期値設定
HD = 1.3  # 直径を測定した高さ

# モデルの定義
with pm.Model() as model:
    
    ### データ関連定義
    ## coordの定義
    # 樹高データのインデックス ※樹の個体別測定別
    model.add_coord('dataH', values=data_H.index)
    # 直径データのインデックス ※樹の個体別
    model.add_coord('dataD', values=data_D.index)
    # 樹種
    model.add_coord('species', values=species_labels)
    # パラメータのσの種類
    model.add_coord('sigmaParam',
                    values=['D', 'H', 'epsilon', 'epsilon_a', 'epsilon_h',
                            'epsilonHmax'])
    
    ## dataの定義
    # 直径データ
    D = pm.Data('D', value=data_D['胸高直径'].values, dims='dataD')
    # H2樹高データ（樹高H－直径を測定した高さHD）
    H2 = pm.Data('H2', value=data_H['樹高'].values - HD, dims='dataH')
    # 樹種インデックス
    speciesIdx = pm.Data('speciesIdx', value=species_idx, dims='dataH')
    # 直径データレファレンス（樹高データに持たせた直径データとのレファレンス）
    refIdx = pm.Data('refIdx', value=data_H['index'], dims='dataH')
  
    ### 事前分布
    ## 6つの標準偏差
    sigmas = pm.HalfNormal('sigmas', sigma=10, dims='sigmaParam')
    
    ## パラメータa
    # 平均
    logabar = pm.Normal('logabar', mu=0, sigma=10)
    # 変量効果 ★非中心化パラメータ化
    offset_a = pm.Normal('offset_a', mu=0, sigma=1, dims='species')
    epsilon_a = pm.Deterministic('epsilon_a', offset_a * sigmas[3], dims='species')
    # パラメータa
    a = pm.Deterministic('a', pt.exp(logabar + epsilon_a), dims='species')
    
    ## パラメータh
    # 平均
    loghbar = pm.Normal('loghbar', mu=0, sigma=10)
    # 変量効果 ★非中心化パラメータ化
    offset_h = pm.Normal('offset_h', mu=0, sigma=1, dims='species')
    epsilon_h = pm.Deterministic('epsilon_h', offset_h * sigmas[4], dims='species')
    # パラメータh
    h = pm.Deterministic('h', pt.exp(loghbar + epsilon_h), dims='species')
    
    ## パラメータHmax
    # 平均
    logHmaxbar = pm.Normal('logHmaxBar', mu=4, sigma=4)
    # 変量効果 ★非中心化パラメータ化
    offset_Hmax = pm.Normal('offset_Hmax', mu=0, sigma=1, dims='species')
    epsilon_Hmax = pm.Deterministic('epsilon_Hmax', offset_Hmax * sigmas[5],
                                    dims='species')
    # パラメータHmax
    Hmax = pm.Deterministic('Hmax', pt.exp(logHmaxbar + epsilon_Hmax),
                            dims='species')
    
    ## 真の直径Dbar
    Dbar = pm.HalfNormal('Dbar', sigma=100, dims='dataD')
    
    ## logHbar
    # 誤差 ★樹の個体単位
    offset_ep = pm.Normal('offset_ep', mu=0, sigma=1, dims='dataD')
    epsilon = pm.Deterministic('epsilon', offset_ep * sigmas[2], dims='dataD')
    # logHbar
    logHbar = pm.Deterministic(
        'logHbar',
        (
            pt.log(a[speciesIdx])
            + h[speciesIdx] * pt.log(Dbar[refIdx])
            + pt.log(Hmax[speciesIdx])
            - pt.log(a[speciesIdx] * pt.pow(Dbar[refIdx], h[speciesIdx])
                    + Hmax[speciesIdx])
            + epsilon[refIdx]
        ),
        dims='dataH'
    ) 

    ### 尤度
    # 直径D
    likeliD = pm.Normal('likeliD', mu=Dbar, sigma=sigmas[0],
                        observed=D, dims='dataD')
    # 樹高H2：H-HD
    likeliH2 = pm.LogNormal('likeliH2', mu=logHbar, sigma=sigmas[1], 
                            observed=H2, dims='dataH')

【モデル注釈】

• coordの定義
  座標に名前を付けたり、その座標が取りうる値を設定できます。
  今回は次の４つを設定しました。

  • 樹高データの行の座標：名前「dataH」、値「行インデックス」

  • 直径データの行の座標：名前「dataD」、値「行インデックス」

  • 樹種の座標：名前「species」、値「樹種の和名」

  • $${\sigma}$$の座標：名前「sigmaParam」、値「$${\sigma}$$の添字含むパラメータの名前」

• dataの定義
  次の４つを設定しました。

  • 直径データ D

  • 樹高から測定高さ1.3を引いたデータ H2

  • 樹種インデックス speciesIdx

  • 直径データレファレンス refIdx
    （樹高データに持たせた直径データとのレファレンス）

• パラメータの事前分布・計算値

  • テキストの記載のとおりです。

• 尤度

  • 直径 D は正規分布、樹高-1.3 H2 は対数正規分布に従うとしています。

２．モデルの外観の確認

# モデルの表示
model

【実行結果】
パラメータが多いのでたくさんの数式が並びます。

# モデルの可視化
pm.model_to_graphviz(model)

【実行結果】
標準偏差 sigmas がヘッドで仕切っているみたいな構図、面白いですねｗ

３．事後分布からのサンプリング
NUTSサンプラーに nutpie を利用します。
処理時間は 7分10秒でした。

### 事後分布からのサンプリング ※NUTSサンプラーにnutpieを使用　7分10秒
# 書籍：iter=6000, warmup=3000, chains=デフォルト（おそらく4）
with model:
    idata = pm.sample(draws=6000, tune=3000, chains=4,
                      target_accept=0.95,
                      nuts_sampler='nutpie', random_seed=1234)

【実行結果】
Divergences は０です！

４．サンプリングデータの確認
$${\hat{R}}$$ と有効サンプル数を確認します。
事後分布の収束確認はテキストにならって$${\hat{R} \leq 1.1}$$とします。

### r_hat>1.1の確認
# 設定
idata_in = idata         # idata名
threshold = 1.1          # しきい値

# しきい値を超えるR_hatの個数を表示
(az.rhat(idata_in) > threshold).sum()

【実行結果】
$${\hat{R} > 1.1}$$のパラメータは「0」件です。
全てのパラメータが $${\hat{R} \leq 1.1}$$ であることを確認できました。

有効サンプル数の確認をします。
ess_bulk、ess_tail の両方が 100 以上であることを確認します。
（100 以上というしきい値はとても緩い基準です…）

### 有効サンプル数の確認：ess_bulk が100以上あればOK
(az.ess(idata, method='bulk') < 100).sum()

【実行結果】
ess_bulk $${< 100}$$ のパラメータは「0」件です。
全てのパラメータが ess_bulk $${\geq 100}$$ であることを確認できました。

### 有効サンプル数の確認：ess_tail が100以上あればOK
(az.ess(idata, method='tail') < 100).sum()

【実行結果】
ess_tail $${< 100}$$ のパラメータは「0」件です。
全てのパラメータが ess_tail $${\geq 100}$$ であることを確認できました。

主だったパラメータのトレースプロットを見ておきます。
sigmas、logabar、loghbar、logHmaxBarを描画します。

### トレースプロットの表示
var_names = ['sigmas', 'logabar', 'loghbar', 'logHmaxBar']
pm.plot_trace(idata, var_names=var_names, compact=False)
plt.tight_layout();

【実行結果】
右のトレースプロットに注目します。
sigma D がちょっと怪しいですが…その他のパラメータは綺麗なゲジゲジを描いていると思います。

ここまでの確認により「収束した」とみなしましょうヽ(=´▽`=)ﾉ

分析

１．注目するパラメータの要約統計量の確認

事後統計量算出の準備をします。
ヘルパー関数の定義とidataからMCMCサンプルを取り出しします。

## 事後統計量算出の準備

# 事後統計量算出関数の定義
def calc_post_stats(samples):
    q = [0.025, 0.25, 0.5, 0.75, 0.975]
    return samples.mean(), *np.quantile(samples, q=q)

# MCMCサンプルデータをパラメータごとに抽出
# 取り出しやすい形式に変換
idata_extract = az.extract(idata)
# 樹種共通パラメータ
log_H_max_bar_ = idata_extract.logHmaxBar.data
log_a_bar_ = idata_extract.logabar.data
log_h_bar_ = idata_extract.loghbar.data
sigmas_ = idata_extract.sigmas.data
# 各樹種のパラメータ
a_ = idata_extract.a.data
h_ = idata_extract.h.data
H_max_ = idata_extract.Hmax.data

【実行結果】なし

① 樹種共通パラメータの事後統計量
$${\log \overline{H}_{\text{max}},\ \log \overline{a},\ \log \overline{h},\ \sigma}$$の事後統計量を確認します。

### 推論データの事後統計量の表示　樹種共通パラメータ
# 事後統計量の表示
var_names = ['logHmaxBar', 'logabar', 'loghbar', 'sigmas']
pm.summary(idata, hdi_prob=0.95, var_names=var_names)

【実行結果】

Excelデータのシート「Table S2」にテキストの樹種共通パラメータの事後統計量が掲載されています。
テキストのモデルで生成したMCMCサンプルを用いていると思われます。
引用いたします。

MCMCサンプルをこの表S2と同じ形式で表示して比べてみましょう。

# 樹種共通パラメータの事後分布の平均・分位点 ※Excel 表S2相当
columns = ['平均値', '2.5%', '25%', '50%', '75%', '97.5%']
summary2_df = pd.DataFrame(
    {
        'log H Max': calc_post_stats(log_H_max_bar_),
        'log a': calc_post_stats(log_a_bar_),
        'log h': calc_post_stats(log_h_bar_),
        'σ H': calc_post_stats(sigmas_[1]),
        'σ D': calc_post_stats(sigmas_[0]),
        'σ': calc_post_stats(sigmas_[2]),
        'σ a': calc_post_stats(sigmas_[3]),
        'σ h': calc_post_stats(sigmas_[4]),
        'σ HMax': calc_post_stats(sigmas_[5])
    }, index=columns
).T
display(summary2_df.round(2))

【実行結果】
概ねテキストと同じ結果が得られています。よかったです。

テキスト掲載のパラメータ推定値と照らし合わせてみます。

• 樹高の測定誤差のばらつき$${\sigma_H}$$の事後平均

  • テキストは$${0.03}$$

  • PyMCモデルは$${0.033}$$

• 直径の測定誤差の標準偏差$${\sigma_D}$$の事後平均

  • テキストは$${0.52cm}$$

  • PyMCモデルは$${0.507cm}$$

② 各樹種共通パラメータの事後統計量
$${a,\ h,\ H_{\text{max}}}$$の事後統計量を確認します。
パラメータ$${a}$$を見ましょう。

### 推論データの事後統計量の表示　◯a、h, Hmax
var_names = ['a']
pm.summary(idata, hdi_prob=0.95, var_names=var_names)

【実行結果】

パラメータ$${h}$$を見ましょう。

### 推論データの事後統計量の表示　a、◯h, Hmax
var_names = ['h']
pm.summary(idata, hdi_prob=0.95, var_names=var_names)

【実行結果】

パラメータ$${H_{\text{max}}}$$を見ましょう。

### 推論データの事後統計量の表示　a、h, ◯Hmax
var_names = ['Hmax']
pm.summary(idata, hdi_prob=0.95, var_names=var_names)

【実行結果】

Excelデータのシート「Table S3」にテキストの各樹種のパラメータの事後統計量が掲載されています。
テキストのモデルで生成したMCMCサンプルを用いていると思われます。
引用いたします。

MCMCサンプルをこの表S3と同じ形式で表示して比べてみましょう。

# 各樹種のパラメータa, h, HMax の事後分布の平均・分位点 ※Excel 表S3相当

# a, h, HMax ｎ事後統計量を算出してそれぞれデータフレーム化
tmp1_df = pd.DataFrame([calc_post_stats(s) for s in a_])
tmp2_df = pd.DataFrame([calc_post_stats(s) for s in h_])
tmp3_df = pd.DataFrame([calc_post_stats(s) for s in H_max_])
# ３つのデータフレームを結合
summary3_df = pd.concat([tmp1_df, tmp2_df, tmp3_df], axis=1)
# カラムとインデックスの設定
multi_cols = pd.MultiIndex.from_arrays([['a']*6 + ['h']*6 + ['HMax']*6, columns*3])
summary3_df.columns = multi_cols
summary3_df.index = species_labels
# データフレームの表示
display(summary3_df.round(2))

【実行結果】
概ねテキストと同じ結果が得られています。よかったです。

🌳🌳🌳

２．直径と樹高の相対成長関係

テキストの図5.3で示されるような樹種ごとの直径と高さの関係の可視化を行ってみます。

まず、直径Dが0.1～65cmのときの樹高の推定値（中央値、95%CI）を計算します。
複雑に見えますが、ベイズモデリングの数式に則っているので、「テキストのモデリング」項の数式を参照しながらコードを確認するのが良いと思います。

### 図5.3描画のための各種値算出

## 樹高計算関数の定義
def calc_height(D, a, h, H_max, H_D=1.3):
    return a * D**h * H_max / (a * D**h + H_max) + H_D

# 樹高の事後分布サンプルの算出
D = np.linspace(0.1, 65, 101)    # x軸の直径値
n_species = len(species_labels)  # 樹種の数
heights = np.array([
    [calc_height(d, a_[i], h_[i], H_max_[i]) for d in D] for i in range(n_species)
])

# 樹高の事後中央値と事後95%CI区間の算出
med = np.median(heights, axis=2)
ci_lower, ci_upper = np.quantile(heights, q=(0.025, 0.975), axis=2)

# 樹高の観測値の樹単位の平均値の算出
data_H_mean = data_H.pivot_table(
    index=['和名', 'No', '胸高直径'], values='樹高', aggfunc='mean', sort=False
).reset_index()

# データフレームの表示
display(data_H_mean)

【実行結果】
この表の値が散布図のデータ点になります。

テキスト掲載の「カナメモチ」と「コナラ」に注目します。
図 5.3 に相当する図を描画します。
観測値の散布図、樹高推論値の中央値の折れ線グラフ、95%CIの塗りつぶしを行っています。

### 推論値を用いた直径と高さの関係の描画：カナメモチとコナラ ※図5.3に相当

# 描画領域の指定
fig, axes = plt.subplots(1, 2, figsize=(8, 3.5), sharey=True, tight_layout=True)
# 樹種ごとに描画処理を繰り返す
for i, (tree, ax) in enumerate(zip(['カナメモチ', 'コナラ'], axes.flat)):
    # 樹ごとの直径・高さの観測値の描画（高さは平均値）
    ax.scatter(data_D[data_D['和名']==tree]['胸高直径'],
               data_H_mean[data_H_mean['和名']==tree]['樹高'],
               alpha=0.5, zorder=10)
    # 事後中央値の描画
    idx = int(np.where(species_labels == tree)[0])
    ax.plot(D, med[idx].squeeze(), color='blue')
    # 事後95%CI区間の描画
    ax.fill_between(D, ci_lower[idx], ci_upper[idx], color='tomato', alpha=0.3)
    # 修飾
    ax.set(title=f'{tree}')
    ax.grid(lw=0.5)
# 全体修飾
fig.supxlabel('直径(cm)', fontsize=12)
fig.supylabel('高さ(m)', fontsize=12)
plt.show()

【実行結果】

最大樹高 $${H_{\text{max}}}$$ の事後統計量も表示しましょう。

# 最大樹高 Hmax の事後統計量の表示
summary3_df.loc[['カナメモチ', 'コナラ'], 'HMax'].round(2)

【実行結果】

上の図と事後中央値（50%）に基づくテキストの分析をお借りします。

• カナメモチの最大樹高の事後中央値：

  • テキスト： 12.4 m（ExcelのS3表は 12.9 m）、このモデル： 13.0 m

  • この高さは現実的なものだそうです。

• コナラの最大樹高の事後中央値：

  • テキスト： 22.8 ｍ、このモデル： 22.9 ｍ

  • 30 m を超える観測値は外れ値的になっていて、事後中央値は 23 m 弱という推定値は、データを取得した場所が樹木の成長がよくない地質であったことを考慮すると「悪くない推定値」だそうです。

せっかくなので全樹種を確かめてみましょう。

### 推論値を用いた直径と高さの関係の描画：全樹種版

# 描画領域の指定
fig, axes = plt.subplots(5, 3, figsize=(12, 15), sharey=True, tight_layout=True)
# 樹種ごとに描画処理を繰り返す
for i, (tree, ax) in enumerate(zip(species_labels, axes.flat)):
    # 樹ごとの直径・高さの観測値の描画（高さは平均値）
    ax.scatter(data_D[data_D['和名']==tree]['胸高直径'],
               data_H_mean[data_H_mean['和名']==tree]['樹高'],
               alpha=0.5, zorder=10)
    # 事後中央値の描画
    ax.plot(D, med[i], color='blue')
    # 事後95%CI区間の描画
    ax.fill_between(D, ci_lower[i], ci_upper[i], color='tomato', alpha=0.3)
    # 修飾
    ax.set(title=f'{species_labels[i]}')
    ax.grid(lw=0.5)
# 全体修飾
fig.supxlabel('直径(cm)', fontsize=12)
fig.supylabel('高さ(m)', fontsize=12)
plt.show()

【実行結果】
青点が観測値、青線が推論値の中央値、赤塗りが推論値の95%CIです。
直径 60 cm まで表現していますが、樹種によっては観測値が 60 cm より小さい場合があり、この場合には大きな推論値は外挿となります。

最大樹高 $${H_{\text{max}}}$$ の事後統計量を表示しましょう。

# 最大樹高 Hmax の事後統計量の表示
summary3_df['HMax'].round(2)

【実行結果】

サカキ、ヒサカキ、ネズミモチといった「観測値が直径：小・樹高：低」の樹種について、最大樹高の推定値は高く推論されているように感じます。
気になって、Google の AI 検索で調べたところ…

との検索結果を得ました（2026年3月16日現在）。

【まとめ】
パラメータのMCMCサンプルデータを用いて、モデルの数式に当てはめて、目的変数である「樹高」を推論値を計算できました。
ベイズのMCMCサンプリングデータの活用方法をまた１つ学ぶことができました。

idata のファイル保存は省略します。

以上で第５章は終了です。

おわりに

予測

今回の「木の高さを推論するモデル」には、機械学習に通じるものを感じました。
観測データの「説明」・「メカニズムの把握」にとどまることなく、未知データの「予測」に用いる目的を持っているからです。

予測に寄与するベイズモデルを自在に構築できる「ベイズ使い」になりたいです。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。