文献：LSM uncertainty analysis

PyLandslide: A Python tool for landslide susceptibility mapping and uncertainty analysis - ScienceDirect

AI要約

背景
地すべりリスク評価において、専門家の主観的判断による重み付けが一般的だが、バイアスや不確実性が問題。
機械学習を用いた統計的アプローチが発展しているが、重みの不確実性の定量化が不十分。
既存ツール(LSM Tool Pack、LSAT PM等)は存在するが、不確実性分析機能が限定的。

手法
1.ジニ不純度に基づく特徴量重要度(重み)の計算
1)ランダムフォレスト内の各決定木において、ある特徴量で分岐する際のジニ不純度の減少量を計算
2)全ての木について、各特徴量による不純度減少量を合計
3)各特徴量の重要度 = その特徴量による不純度減少の合計 / 全特徴量による不純度減少の総和

2.不確実性分析の手法
20,000回の反復計算を実施
各反復で地すべり・非地すべり地点の80%をランダムサンプリング
accuracy 0.75以上のモデルのみの特徴量重要度（feature importance）を記録
重要度の範囲を集計し、限られたインベントリデータに起因する不確実性を定量化

3.LSI（Landslide Susceptibility Index）の算出
ある地点が地すべりに対してどれだけ脆弱であるかを数値化した指標
PyLandslideでは、複数の要因を重み付けして統合することで、空間的な地すべり発生リスクを評価

LSI = Σ(Wi × Si)

LSI: 地すべり感受性指数
Wi: 要因iの重み（weight）
Si: 要因iのスコア層（score layer）
n: 考慮する要因の総数

スコア層（Si）の算出方法
連続データ型要因（降水量、傾斜、道路距離など）：
分位数による分類: データを11個の分位数（quantile）に分類
地すべり件数の計算: 各クラス内で発生した過去の地すべり件数を集計

降水量・傾斜: そのクラス以下で発生した地すべりの累積パーセンテージ
道路距離: そのクラス以上で発生した地すべりの累積パーセンテージ
降水量が17~80mmのクラス: 80mm以下で発生した地すべりの累積割合がスコア
道路距離が143~286mのクラス: 143m以上の距離で発生した地すべりの割合がスコア

カテゴリカルデータ型要因（土地被覆、岩相など）：
各カテゴリ内で発生した地すべり件数を計算
各カテゴリのスコア = そのカテゴリ内で発生した地すべりのパーセンテージ

標準化
全ての要因について、スコアを0~100の範囲に線形正規化:
0: その要因クラスの最低寄与度
100: その要因クラスの最高寄与度

LSI分類基準
Very low: 0 ≤ LSI ≤ 20
Low: 20 < LSI ≤ 40
Moderate: 40 < LSI ≤ 60
High: 60 < LSI ≤ 80
Extremely high: 80 < LSI

感度分析の実施回数
歴史的降水条件について200回のランダムな重みの抽出を実施
9つの気候予測について各200回ずつ実施
内訳：
歴史的降水条件（1981-2023年）：200回
将来予測（2041-2050年）：9つの気候シナリオ × 各200回 = 1,800回
合計：200 + 1,800 = 2,000回

9つの気候シナリオ：
3つの気候モデル（ACCESS-ESM1-5、MRI-ESM2-0、UKESM1-0-LL）
3つの共有社会経済経路（SSP126、SSP245、SSP585）
3モデル × 3シナリオ = 9つの組み合わせ
この結果、LSIも範囲として表現される

結果
イタリアにおける6要因の重要度(重み)
道路からの距離: 0.43〜0.52(中央値0.47)
・最も影響が強い
・過去の地すべりの85%が道路から143m以内で発生

傾斜: 0.16〜0.23(中央値0.20)
・2番目に重要
・過去の地すべりの70%が傾斜9%以上の地域で発生

土地被覆: 0.10〜0.13(中央値0.12)
地質: 0.07〜0.09(中央値0.08)
降水量: 0.05〜0.07(中央値0.06)
TWI(地形湿潤指数): 0.05〜0.08(中央値0.06)

モデル精度
訓練データでの全体精度: 0.80〜0.82
テストデータでの全体精度: 0.75〜0.80
過学習は認められない

ｌSI評価
歴史的期間(1981-2023):
「極めて高い」: 7.8〜9.5%
「高い」: 23.8〜26.8%

将来予測(2041-2050):
「極めて高い」: 5.3〜7.6%(減少傾向)
「高い」: 21.5〜28.3%

イタリア北西部と南部で大幅に低下
北東部で増加

考察
人為的要因(道路)が最大の影響
道路建設に伴う斜面切土や人間活動が斜面安定性を低下

不確実性の重要性
重みの不確実性範囲を考慮することで、より堅牢な意思決定が可能
インフラ投資計画において、重みの不確実性を考慮した堅牢な地すべりリスク評価が可能

研究の限界
気候変動の影響評価
降水量変化のみを考慮、温度・山火事・植生変化は未考慮

入力データの精度
ITALICAデータベースの位置精度は最大5.6kmの誤差
より高精度なデータが望ましい

閾値の設定
accuracy0.75の閾値は本研究で設定したが、絶対的な基準ではない
ユーザーのニーズに応じて調整可能

国外では、LSM作成において機械学習を用いるアプローチは既に一般化していると言えますが、そこに不確実性評価が入っています。といっても、重要度の示し得る範囲を利用しているだけですから、新たなツールは不要です。
この文献のLSIの定義が正解かどうかはわかりませんが、少なくとも確率を導入する方針については同意です。

文献：Probabilistic rainfall thresholds for landslide occurrence using a Bayesian approach

 Probabilistic rainfall thresholds for landslide occurrence using a Bayesian approach - Berti - 2012 - Journal of Geophysical Research: Earth Surface - Wiley Online Library

AI要約＋α

背景
従来の地すべり予測手法は、地すべりを引き起こした降雨イベントのみを考慮し、決定論的閾値を提供してきた。しかし、斜面安定性は降雨だけでなく、多数の要因の組み合わせで決まるため、同じ降雨条件でも異なる結果(地すべり発生/非発生)が生じる。経験的モデル:降雨強度-継続時間(ID)閾値が一般的だが、不確実性の定量化が不十分である。

研究の対象地域
イタリア・エミリア=ロマーニャ州の山岳地域(12,000 km²)
1939年以降の4,000件以上の地すべり履歴データ
176箇所の雨量計による日降水量データ

手法
1. ベイズの基本手法(1次元解析)
P(A|B) = [P(B|A) × P(A)] / P(B)
P(A|B): 事後確率 - 降雨Bが発生したときに地すべりAが起こる確率(求めたい値)
P(B|A): 尤度 - 地すべりが発生したときに降雨Bが観測される確率
P(A): 事前確率 - 降雨に関係なく地すべりが発生する確率
P(B): 周辺確率 - 地すべりに関係なく降雨Bが観測される確率

実際の計算では、以下のように相対度数で近似:
P(A) ≈ NA / NR  (地すべり発生数 / 全降雨イベント数)
P(B) ≈ NB / NR  (降雨Bの発生数 / 全降雨イベント数)
P(B|A) ≈ N(B|A) / NA  (降雨Bで発生した地すべり数 / 全地すべり数)P
(A|B) ≈ N(B|A) / NB

一変量例（強度 I>40 mm/日）
・事前 landslide 確率：P(A)=5/20=0.25
・強度大の割合：P(B)=P(I>40)=9/20=0.45
・発生時の強度大割合（尤度）：P(B|A)=4/5=0.8

→ P(A|I>40)=P(B|A)·P(A)/P(B)=0.8·0.25/0.45≃0.44

降雨B発生時の地すべり確率は44%である。これは尤度P(B|A)=80%とは異なることに注意が必要。事前確率と周辺確率を考慮する必要がある。

2. 2次元ベイズ解析
2つの変数B(降雨強度I)とC(降雨継続時間D)を考慮:
P(A|B,C) = [P(B,C|A) × P(A)] / P(B,C)
ここで、P(B,C)は2変数の同時確率(ある強度と継続時間の組み合わせが観測される確率)を表す。

各セルでの確率計算
二変量例（I>50 mm/day かつ D<=1 day）
・事前 landslide 確率：P(A) = 5/20 = 0.25
・降雨条件を満たす割合：P(B,C) = P(I,D) = 4/20 = 0.20
・発生時の当該領域割合（尤度）：全地すべり5件中、この条件で発生2件
   P(B,C|A) = P(I,D|A) = 2/5 = 0.40

→ P(A|I>50, D≤1)= P(A|B,C) = [P(B,C|A) × P(A)] / P(B,C) =0.40·0.25/0.20=0.50

I-D平面全体の確率分布ヒストグラムを作成。
強度と継続時間の相互作用効果を把握。
確率の空間分布パターンを視覚化。
高リスク領域の明確な識別。

結果
1次元ベイズ解析の結果
有意な変数:総降雨量E、継続時間D、平均強度I
P(B|A)とP(B)の分布が明確に異なる
P(A|B)が事前確率P(A)=0.005を大きく上回る
特に降雨強度Iが最も有意:I>100mm/dayでP(A|I)=0.28

非有意な変数:先行降雨AE14, AE30
P(B|A) ≈ P(B)
P(A|B) ≈ P(A)

研究地域では先行降雨は地すべり発生と相関が低い。
降雨の激しさ(総量、継続時間、強度)とともに地すべり確率が増加。
ただし、極端な値では確率が減少する傾向(サンプルサイズの問題と定義バイアスによる)。

2次元ベイズ解析の結果
特定のI-D値で確率が急激に増加。これはシステムの状態の根本的変化、物理的閾値の存在を示唆。不確実性を含めた警報システムの構築に有用。

考察

1. ベイズ手法の利点
統計的厳密性：事前確率と周辺確率を考慮:尤度P(B|A)だけでなく、降雨の頻度P(B)と地すべり発生率P(A)を組み込む。
認知バイアスの回避:人間の直感的判断は事前確率を無視しがち(「80%の地すべりがI>50で発生」≠「I>50で地すべり確率80%」)。
不確実性の定量化:0〜1の連続値で確率を表現、決定論的手法の曖昧性(「閾値超過時に何が起こるか?」)を解消、信頼区間により推定の信頼性を評価

2. 先行降雨の非有意性
一般に細粒土壌では先行降雨が斜面安定性に重要とされるが、本研究では相関が低い。研究地域の地すべりは主にイベント降雨に対する急速な水文応答が支配的。60%の地すべりが降雨終了時またはその直後に発生。
深層地すべりでも、長期的な要因よりも短期的な降雨イベントが引き金となる。

3. 物理的閾値の示唆
I-D平面の特定領域で地すべり確率が急激に上昇。これはシステムの状態変化(安定→不安定)を示唆。
確率的手法であっても、背後に物理的メカニズム(臨界間隙水圧など)が存在。確率分布のパターンから、真の物理的閾値の存在を推測可能。

4. 方法論的課題と解決策
多発地すべりの扱い
方針:同一降雨による複数地すべりは1イベントとしてカウント。
理由:P(A) = NA/NR > 1を避け、統計的整合性を保つ。

バイアスへの対応
問題:誘因降雨は地すべり発生日で打ち切られるが、非誘因降雨は降雨終了まで継続。
影響:極端な長継続時間での確率推定に影響。
緩和:信頼区間の提示、データの60%が降雨終了時に発生することで影響を軽減。

Landslide発生降雨をベイズ手法で評価した少し古い文献です。基本的には、国内で研究されている降雨数ベースの評価と変わりません。
国内では短期雨量と土壌雨量指数の２軸で評価されます。スネーク曲線を書いてピークで評価したりするのですが、原点から離れるにつれてレアな雨（P(B,C) ≒０）となり、確率が極端に高くなります。これ、補正したいですよね。

Synthetic Control

Synthetic Control は、ベイズ統計を用いた因果推論です。機械学習ではなく、MCMCです。

1. モデル設定（合成コントロールの基本形）

ある「介入クラス」A について、時点 t の平均テスト点数を y_A,t とする。また、介入を受けていない対照クラスを j = 1,…,J とし、それぞれの平均点を y_j,t と書く。
合成コントロールの考え方は「A クラスは、対照クラスの“重み付き平均”で近似できる」とみなすことである。そこで、時点 t における A クラスの期待される点数（＝介入がなかったときの基準）を

μ_t = w_1 y_1,t + … + w_J y_J,t = Σ_j w_j y_j,t

と表す。w_j は 0 以上で和が 1 になるような重みとする。
実際の点数 y_A,t は、体調・偶然・測定誤差などにより μ_t からずれる。このズレを「平均 0、標準偏差 σ の正規分布に従う誤差」と仮定して

y_A,t = μ_t + 誤差_t
誤差_t ~ Normal(0, σ)

とおく。確率変数として書けば

y_A,t|w, σ ~ Normal(μ_t, σ)

である。

2. ベイズの定理と記号

推定したいパラメータを

θ = (w_1,…,w_J, σ)

とまとめて書く。ベイズの定理は以下となる。

p(θ|データ) = p(データ|θ) p(θ) / p(データ)

p(θ|データ)：事後分布（データを見た後に、θ がどの程度になりそうか）
p(データ|θ)：尤度（θ が真なら、このデータはどれくらい出やすいか）
p(θ)：事前分布（データを見る前に、θ がどのような値を取りやすいと考えるか）
p(データ)：正規化定数（全部の θ について積分した値）

実際の計算では p(データ) を明示的に計算する必要はなく、

p(θ|データ) ∝ p(データ|θ) p(θ)

という比例関係として扱うことが多い。

3. 事前分布 p(θ)

重みベクトル w = (w_1,…,w_J) について

w ~ Dirichlet(α_1,…,α_J)

とおく。誤差のばらつき σ については

σ ~ HalfNormal(τ)

とする。これらをまとめると、事前分布は

p(θ) = p(w_1,…,w_J, σ)
= Dirichlet(w | α_1,…,α_J) × HalfNormal(σ | τ)

となる。ここでは w と σ を互いに独立と仮定している。

4. 尤度 p(データ|θ)

合成コントロールの重み w は、「介入前」で学習する。介入前の時点を t = 1,…,T0 とし、A クラスの平均点：y_A,t、各対照クラスの平均点：y_j,t が観測されているとする。モデルでは、各 t について

y_A,t|θ ~ Normal(μ_t, σ)
μ_t = Σ_j w_j y_j,t

であるから、ある θ が与えられたとき、「その θ のもとで、実際に観測された系列 {y_A,1,…,y_A,T0} がどれくらい出やすいか」を表すのが尤度である。
時点 t ごとに独立と仮定すると、介入前データ全体に対する尤度は

p(データ|θ) = Π_{t=1}^{T0} p(y_A,t|θ)
= Π_{t=1}^{T0} Normal(y_A,t ; μ_t, σ)

となる。ここで μ_t は w と対照クラスの観測値 {y_j,t} から一意に決まる。
・μ_t に近い y_A,t が多ければ、この θ の尤度は大きくなる
・μ_t から遠く離れた y_A,t が多ければ、尤度は小さくなる
という直感と対応している。

5. 事後分布 p(θ|データ)

ベイズの定理より、介入前データを見た後の θ の分布（事後分布）は

p(θ | データ) ∝ p(データ | θ) p(θ)
= {Π_{t=1}^{T0} Normal(y_A,t ; μ_t, σ)}
　× Dirichlet(w | α_1,…,α_J) × HalfNormal(σ | τ)

で与えられる。これは介入前の A クラスの点数系列をどれだけよく再現するか（尤度）、もともと重み w や σ をどう想定していたか（事前）の両方をバランスさせた結果である。
この事後分布を式のまま解析的に求めることは難しいため、MCMC（たとえば NUTS というアルゴリズム）を用いて

θ^(1), θ^(2), …, θ^(S) = (w_1^(s),…,w_J^(s), σ^(s)), s=1,…,S

という多数のサンプルを生成し、p(θ|データ) を近似する。それぞれの θ^(s) が「あり得る合成コントロールの重みと誤差の大きさ」を表しており、全体として「どのような w がどの程度ありそうか」を数値的に把握できる。

6. Synthetic Control による介入効果推定

6.1 介入「前」のデータから事後分布 p(θ|データ) を推定
上で述べたように、介入前 t=1,…,T0 のデータを用いて事後分布 p(θ | データ) を求める。

6.2 介入「後」のシンセティック点数 μ_t^(s) を計算
介入後 t = T0+1,…,T の期間に注目する。この期間についても、各対照クラス j の平均点 y_j,t が観測されているとする。それぞれの事後サンプル θ^(s) に対して、時点 t ごとに

μ_post,t^(s) = Σ_j w_j^(s) y_j,t

を計算する。これは「サンプル s に対応する合成コントロールの重みを用いて、介入後の A クラスの“カウンターファクチュアル（介入がなかった場合の想定点数）”を計算したもの」である。したがって、{μ_post,t^(s)}_{s=1,…,S} は「介入がなかったとしたとき、時点 t における A クラスの点数はどのくらいになり得るか」を表す事後予測分布となる。

6.3 介入効果の事後分布を得る
実際に観測された介入後の A クラスの平均点を y_A,t（t > T0）とすると、サンプル s における時点 t の介入効果は

effect_t^(s) = y_A,t − μ_post,t^(s)

と定義できる。これは
・正なら「実際の点数の方が高い（介入がプラス効果）」
・負なら「実際の点数の方が低い（介入がマイナス効果）」
ことを意味する。

時点 t ごとに {effect_t^(s)} の平均
 → その時点の平均介入効果 2.5%点と 97.5%点
 → 95% 信用区間（事後分布の幅）を
求めれば、「いつ」「どれくらい」介入効果が表れたかをベイズ的に評価できる。
さらに、介入後の全期間について effect_t^(s) を平均すれば
・介入後期間全体の平均効果
・その信頼区間（信用区間）
も同じ要領で推定できる。

このようにして、従来の Synthetic Control（複数の対照ユニットの重み付き平均）を「事前分布」「尤度」「事後分布」の枠組みの中にきちんと位置づけ、ベイズ的な不確実性評価と一体化した形で介入効果を推定する。

7. Pythonコード（AI提案）

import pandas as pd
import numpy as np
import pymc as pm
import arviz as az

# --- データ準備 ---
# yA:  shape=(T,)       介入クラスA の時系列
# Y_donors: shape=(T,J) 対照クラス群の時系列
# T0: 介入前の最後の時点（0始まりのインデックス）

T, J = Y_donors.shape
yA_pre       = yA[: T0+1]            # 介入前データ
Y_donor_pre  = Y_donors[: T0+1, :]
yA_post      = yA[T0+1 :]            # 介入後データ
Y_donor_post = Y_donors[T0+1 :, :]

# --- 6.1 介入前データでモデルフィッティング ---
with pm.Model() as sc_model:
    # 事前分布
    w     = pm.Dirichlet("w",     a=np.ones(J))   # w_j ≥ 0, Σw_j=1
    sigma = pm.HalfNormal("sigma", sigma=1.0)     # σ>0

    # 合成コントロールの期待値 μ_t
    mu_lin = pm.math.dot(Y_donor_pre, w)
    mu_pre = pm.Deterministic("mu_pre", mu_lin)

    # 尤度
    pm.Normal("y_obs", mu=mu_pre, sigma=sigma,
              observed=yA_pre)

    # サンプリング
    trace = pm.sample(
        draws=1000,          # 事後サンプル 1000
        tune=1000,           # チューニング 1000
        chains=4,            # チェーン数
        cores=4,             # 並列数
        target_accept=0.9,   # NUTS の受諾率
        max_treedepth=12,
        random_seed=0,
        progressbar=True
    )

# --- 6.2 介入後のカウンターファクチュアルを計算 ---
# ポスターサンプルから w を (S,J) にフラット化
w_samples = trace.posterior["w"] \
    .stack(sample=("chain", "draw")) \
    .values  # shape = (S, J)

# μ_post[s,t] = Σ_j w_samples[s,j] * Y_donor_post[t,j]
# 結果 shape = (S, T_post)
mu_post = (Y_donor_post @ w_samples.T).T

# --- 6.3 介入効果を計算 ---
# effect[s,t] = 実測 yA_post[t] − μ_post[s,t]
effect = yA_post[None, :] - mu_post  # shape = (S, T_post)

# 時点別平均効果と95%信用区間
effect_mean_t = effect.mean(axis=0)                # (T_post,)
effect_hdi_t  = az.hdi(effect, hdi_prob=0.95)      # (T_post, 2)

# 介入後期間全体の平均効果とその信用区間
avg_effect     = effect.mean()                     # scalar
avg_effect_hdi = az.hdi(effect.flatten(), hdi_prob=0.95)  # (2,)

# --- 出力例 ---
print("時点別平均効果    :", effect_mean_t)
print("時点別95%信用区間:", effect_hdi_t)
print("全期間平均効果    :", avg_effect)
print("全期間95%信用区間:", avg_effect_hdi)

Landslide Susceptibility Map using 3D LEM

RegionGrow3D: A Deterministic Analysis for Characterizing Discrete Three‐Dimensional Landslide Source Areas on a Regional Scale - Mathews - 2024 - Journal of Geophysical Research: Earth Surface - Wiley Online Library

In this study, the RegionGrow3D (RG3D) model is developed to broadly simulate the area, volume, and location of landslides on a regional scale (≥1,000 km2) using 3D, limit-equilibrium (LE)-based slope stability modeling. Furthermore, RG3D is incorporated into a susceptibility framework that quantifies landsliding uncertainty using a distribution of soil shear strengths and their associated probabilities, back-calculated from inventoried landslides using 3D LE-based landslide forensics.

動画もあります。
Landslide susceptibility map generated using RegionGrow3D (youtube.com) 

インベントリがあれば、そこから逆算してΦを出す、その確率密度の分布から累積分布を書いて susceptibility とみなす、その確率でどの大きさの地すべりがどこに分布するかを推定する。このような土質力学ベースの計算が可能な Matlab ツールをUSGSが提供しています。

文献では、すべり層厚を決めるのに隆起量などをパラメータとした別のシミュレーションを利用していますが、この辺りはインベントリの平均層厚でも構わないでしょう。

すべる・すべらないのバイナリを拡張し確率を持たせた点、thresholdに応じて地すべりの大きさを扱える点が魅力です。日本は完全に置いてけぼりです。

コルモゴロフ・スミルノフ検定

2群のデータを比較するためにヒストグラムを作成しました。
それを眺めながら、分布形状が似ている、似ていないを判断したいのですが、主観だと説明性に欠けます。ということで統計的検定の出番です。

以前、2標本（対応ナシ）の条件で Mann-Whitney's U test を使用しました。
https://phreeqc.blogspot.com/2020/12/mann-whitneys-u-test.html

今回は2群のサンプル数が異なるとともに、サンプルが非常に多く何らかの確率分布を仮定できそうでした。
調べてみると、ありますね。
・2標本のコルモゴロフ・スミルノフ検定
two-sample Kolmogorov-Smirnov(K-S) test

手順は、①2群を同じ階級幅で区分して累積度数を求める。②累積度数を累積確率に変換し、その差の絶対値を計算。③絶対値の最大値を使って検定統計量を求め、④棄却限界値と比較し判定する。棄却限界値は自由度ｆ=2のχ2分布から求めます。

これは、好都合。ヒストグラムを全面積で正規化し、縦軸を確率密度で記載していたので処理は半分終わっています。あとは累積確率に直して差の最大値を求めるだけのようなもの。χ2分布の適用も、おかしくはない分布ですから、これで決定です。というか、2群のヒストグラムを書く際に、ルーチンワークとして求めておけば良い内容です。Python コードを残しておきましょう。

def hist_plot(col,x1,x2,label,bins):
    fig = plt.figure(figsize=(15,5))
    ax = fig.add_subplot(1,1,1)

    freq,bin_list,patch=ax.hist([x1,x2], density=True, alpha=0.5 ,edgecolor='black', bins=bins, label=label)
    ax.set_xlabel(col,fontsize=14)
    ax.set_ylabel('確率密度',fontsize=14)
    ax.grid(which='major',color='gray',linestyle='-')
    ax.grid(which='minor',color='gray',linestyle='--')
    ax.legend(borderpad = 0.8,framealpha=1.0,shadow=True)
    
    #2標本コルモゴロフ・スミルノフ検定（χ２分布、各標本数40以上、対応の無い2標本）
    df_p=pd.DataFrame({col:bin_list[1:],label[0]+'_確率密度':freq[0],label[1]+'_確率密度':freq[1]})
    df_p_sum=df_p.cumsum()*bin_list[1]
    df_p_sum=df_p_sum.drop(col,axis=1)
    df_p_sum.columns  = [label[0]+'_累積確率', label[1]+'_累積確率']
    df_p_sum['D']=abs(df_p_sum[label[0]+'_累積確率']-df_p_sum[label[1]+'_累積確率'])
    df_p=pd.concat([df_p, df_p_sum], axis=1)

    n1=len(x1)
    n2=len(x2)
    Dmax=df_p['D'].max()
    Dmax2=Dmax*Dmax
    T=4*Dmax2*(n1*n2)/(n1+n2)
    if T<5.99146:
        comment='有意水準5％で、差があるとは言えない。'
    else:
        comment='有意水準5％で、差があると言える。'        
    display(df_freq,'Dmax',Dmax,'T',T,comment)

Dが小さいと、差があるとは言えない（H0を棄却しない）

H0：帰無仮説（差がない）
H1：対立仮説（差がある）

σレベル

σレベルが出てきました。

(x-平均)/標準偏差

あれ？標準化ですよね。
調べてみると、標準得点の z score と呼ばれるそうで、偏差値の仲間でした。統計と機械学習で呼び方というか、使われ方が違うのかな？

確率年よりも平易で処理が楽そう、というのが第一印象。
Today's Earth のリスク値としても利用されているようです。

有名なのでしょうね。避けずに学んでおけばよかった。

確率密度関数のフィッテイング

確率密度関数のデータへのフィッテイングは SciPy の各分布に実装されています。

最尤推定法のようですが、一度に多数の関数を比較するようにはなっていません。
何かないか？と探してみると、すぐにヒット。皆さん、考えることは同じようです。
https://stackoverflow.com/questions/6620471/fitting-empirical-distribution-to-theoretical-ones-with-scipy-python
distfit は動きませんでした。が、tmthydvnprt & Melquíades Ochoa のコードが使いやすいですね。最小二乗法を用いた総当たり検索です。

手元のデータに当てはめると、数時間かかっても計算し続けていましたので、関数を減らして再実施。無事に最良の関数を選んでくれました。ありがたい。

以前、確率降水量の計算でヒストグラムを作った際に、残差の小さい関数を選択すれば良いのではないかと考えていました。
https://phreeqc.blogspot.com/2019/06/blog-post_20.html
ひとまず、どのようなデータでも計算できるようになりました。

あとは最適化手法でしょうか。

地下水 と MCMC その2

Using Markov Chain Monte Carlo to quantify parameter uncertainty and its effect on predictions of a groundwater flow model - ScienceDirect 

数千回もの浸透流計算が必要となるパラメータ推定手法は、実務向きではありません。https://phreeqc.blogspot.com/2020/12/mcmc.html
この文献では、既存の知見により浸透流計算を省略しています。具体的には、尤度関数内の二乗誤差、誤差分散計算の代替えとして、涵養量と透水係数の比 R/K との相関性を利用します。

overcome this obstacle, 300 realizations of Hassan et al.’s(2002) model were studied and it was found that the terms needed in the likelihood function, namely the sum of squared errors and the error variance are highly correlated to the recharge–conductivity ratio, R/K.

この文献の前に、300 realizations の計算があるのでしょう。それに気づいて、途中で細かく読むのをやめてしまいました。R/K などの相関関係は他の現場で使えたとしても、やはり事前の計算による確認が必要なのでしょうね。

おっ！と思ったり、うーん、と思ったりした文献でした。

 

Mann-Whitney's U test

Mann-Whitney's U test を実施。

scipy.stats.mannwhitneyu を利用していたのですが、ランキングのつけ方が書かれていません。Σ の計算に n(n+1)/2 が含まれていますので、１位、２位、２位、４位のような同順位のつけ方はまずい。そこは正しく計算してくれているのだろうと思いながら、初めて使うので確認。

ソースでは rankdata が使用されています。

ranked = rankdata(np.concatenate((x, y)))

https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rankdata.html?highlight=rank#scipy.stats.rankdata

>>> from scipy.stats import rankdata
>>> rankdata([0, 2, 3, 2])
array([ 1. ,  2.5,  4. ,  2.5])
>>> rankdata([0, 2, 3, 2], method='min')
array([ 1,  2,  4,  2])
>>> rankdata([0, 2, 3, 2], method='max')
array([ 1,  3,  4,  3])
>>> rankdata([0, 2, 3, 2], method='dense')
array([ 1,  2,  3,  2])
>>> rankdata([0, 2, 3, 2], method='ordinal')
array([ 1,  2,  4,  3]) 

はい、一番上の平均 (default is ‘average’) が利用されていました。OKです。

比較する母集団の平均、標準偏差、サンプル数の違いはどのように影響するのでしょうか？簡単な正規分布のプラス側を取り出してチェックしてみました（プロットは正規化済み）。シード値は設けず、何度かテストを繰り返して両側・有意水準5％にかかわる傾向を見てみました。

x1=np.random.normal(1, 10, 1000)
x2=np.random.normal(1, 10, 1000)
x1=x1[x1 > 0]
x2=x2[x2 > 0]   

当然ですが、棄却されない分布です（稀に棄却されるのが気になりますが)。

2群が同じ平均・標準偏差で、一方のサンプルサイズを1000から100000に増やしたケース：検定結果に大きく影響せず（棄却され難い）。

2群が同じ平均・サンプルサイズで、一方の標準偏差を10から20に増やしたケース：検定結果に影響（棄却され易い）。

2群が同じ標準偏差・サンプルサイズで、一方の平均を1から２に増やしたケース：検定結果に大きく影響せず（棄却され難い）。

 

2群が同じ標準偏差・サンプルサイズで、一方の平均を1から3に増やしたケース：検定結果に影響（棄却され易い）。

2群が同じ標準偏差で、一方の平均を1から3、サンプルサイズを1000から100000に増やしたケース：検定結果に影響（棄却され易い）。

ランダムサンプリングをしているため、棄却・棄却されないという結果は固定になりません。それを繰り返し実行しただけですが、比較的容易に傾向が見えてきました。
平均2、３は見た目が微妙ですね。やや大きいかな？というような違いでも「有意差がないとは言えない」という傾向を出せるようです。人によってはやや大きい、あるいは同じように見える、などと意見が分かれるような分布ですが、ある考え方に従い再現性のある答えを出してくれるのはありがたい。

最近では、土研交流研究員の方が景観に配慮した護岸に関する類似性の評価に、この検定を利用されていたように思います。アイデア次第でいろいろな分野に適用できそうです。

好きではなかった統計的手法ですが、その合理性をもっと利用すべきでしょうとも考えるようになってきました。

地下水 と MCMC

ようやくMCMC を概ね理解し動かせるようになったので、地下水モデル10章の不確実性に関する部分へ少しだけ立ち戻ることに。
https://phreeqc.blogspot.com/2019/07/4.html

MCMC関連の文献を探してみると、簡単な計算をされている事例がありました。第2著者以降は日本の方ですね。
Julien Boulange, Hirozumi Watanabe, Shinpei Akai (2017) A Markov Chain Monte Carlo technique for parameter estimation and　inference in pesticide fate and transport modeling

• 気候、水収支、土壌、農薬特性などの40以上のパラメータから、農薬濃度の予測精度に大きな影響を与えることが報告されている以下の4パラメータを選択
  
• 農薬溶解速度（kdiss）、水田土壌中の農薬の一次分解速度（kbio）、農薬の脱離速度（kdes）、農薬分配係数（kd）
• これらの事前確率分布を作成（最大・最小間で一様分布）
• 尤度は正規分布(E(θ)は観測値と予測値の平均2乗誤差)
• MCMC で事後確率分布を作成

尤度を計算する前に採択されたパラメータでシミュを回す必要があります。大きなモデルを扱う実務向きではないですね。

実務で扱う地下水モデルで不確実性を検討するには、PEST一択なのでしょうか？確かに、市販のパッケージソフトには組み込まれており、使い勝手も良かった印象はあります。加えて、逆解析を古くから扱っているせいか、Bayes も良く出てきます。海外では当たり前なのかな？

もう少し、文献等を読んでみましょう。

スパースモデリング

EXCELをよく使っていた頃、データに多項式をフィッティングする場面で何次式を使うかは感覚的でした。

低次だと合いませんし、高次だとノイズまで合わせてしまいます（いわゆる過学習）。ちょうどよいところを探す方法は、、、古くからありましたね。

LASSO
https://scikit-learn.org/stable/modules/linear_model.html#lasso
\begin{align*}
\min_{w}{ \frac{1}{2n_{\text{samples}}} ||y - X w||_2 ^ 2 + \alpha ||w||_1}\end{align*}

この損失関数はモデルの適合性を表す誤差項と、複雑さを表すペナルティー項で構成されています。αは求めたいパラメータwを制御する上位のパラメータ（ハイパーパラメータ）。αを大きくするほど大きなL1ノルムに対してペナルティーがかかる仕組みです。するとL1ノルムは小さくならざるを得なくなり、スパースな解が導かれる（が適合性は低くなりやすい）という仕組み。逆にαを小さくするほどL1ノルムを大きく取れるようになり、適合性が向上しやすい（が複雑で、過学習に陥りやすい）。αをどのように調整・採用するかがキモです。

その判断基準は以下の通り。

・情報量基準（AIC、BICなど）
・cross‐validation

情報量基準の利用では客観性に優れるものの感覚的にずれていると感じることがあるでしょう。一方、cv は理屈でなく現実的。ベースを統計にしているか、機械学習にしているかなどで好き嫌いが分かれそうです。

スパースを仮定した場合、非ゼロのwを増やすと組み込まれるXも高次まで多く含まれて適合性は良くなる（が複雑になる）といった表現が理想です。が、L0ノルムを扱うと組み合わせ（計算コスト）が膨大になるとのこと。そこでα＊L1ノルムで代用しているそうです。
ちなみに、L1、L2 は機械学習でも多用されていましたね。回帰だと以下の通りです。

・L0正則化・・・非ゼロノルムの個数
・L1正則化・・・L1ノルム利用、LASSO回帰
・L2正則化・・・L2ノルム利用、Ridge回帰
・L1＋L2正則化・・・Elastic Net

scikit-learn を利用すれば、これらのモデリングは容易。ありがたいですね。

 ==========================================
ラグランジュの未定乗数法を適用した場合の図や説明も web 上には多くあります。そちらの方が直感的で理解し易いでしょう。

関連内容はコチラ↓
https://phreeqc.blogspot.com/2019/06/blog-post_19.html
https://phreeqc.blogspot.com/2020/02/blog-post.html

MCMC

MCMCを実装するライブラリはいくつかあるようです。

PyMC3:Theano
PyMC4:TensorFlow
Pyro:PyTorch
NumPyro:JAX

手元の環境で動いたのはPyroのみ。環境構築に関してはシビアなようです。

web上のサンプルを見てみますと、書き方はほぼPyMC3と同じでした。動かすのは簡単（traceplot がないのは残念。seaborn で書きました）。

が、写経し終わると理解できていない箇所が露わに。

 y_model = pyro.sample('y_model', dist.Normal(a*x + b, sd), obs=y)

ここの obs の計算内容がわかりません。ここがキモなのはわかっているのですが、マニュアルに詳細が書かれていませんし、ソースも追えません。尤度×事前確率にどのようにかかわっているのでしょう？
おそらく、理解している方々が必要として使われているので、学び方が逆なのでしょうね。

ま、他の部分は理解できていることが分かりました。メモしておきましょう。

 **************************************
20201123 MCMC メモ
 
データ駆動科学

• データを出発点と考える=因果律を遡る
• 従来：変形係数を確定(モデルを確定)→応力、ひずみがノイズを含み確率的分布（→安全率を含めた設計）
• データ駆動：応力、ひずみが確定→変形係数が確率的分布（逆解析的発想）
• データｘ,y（ひずみ、応力）が得られたとき、パラメータa（変形係数）が取り得る確率分布（事後確率）を求めたい。
• 必要なのは、①パラメータaを与えたときにyが取り得る確率分布(尤度)、②aの確立分布（事前確立）

取得データD[x,y]

• 確定x：ひずみx
• 確定y：応力y

予測分布モデル　（例）物理モデル+ノイズ
　物理モデル y=ax+b

• 未確定a:平均μa,標準偏差σaとして設定(事前確率P(a)=N(μa,σa2)を仮定)
• 未確定b:平均μb,標準偏差σbとして設定(事前確率P(b)=N(μb,σb2)を仮定)
• 未確定σn:一様分布(U([n])を仮定)
  ←いずれもデータDを見て取りうる範囲を仮定
  

　ノイズ

• 重畳ノイズを正規分布として設定（尤度P(y_model)=N(ax+b,σn2)を仮定
  

予測分布モデルから θ[a:変形係数、b：切片、sd:標準偏差]の事後確率 p(θ|D)を求める。

• p(θ|D) = p(D|θ)p(θ)/p(D) ∝ p(D|θ)p(θ)
  
• 事後確率 ∝ 尤度×事前確率
• 事前確率、尤度に正規分布を仮定した場合、事後分布を解析的に求めることが可能（共役事前分布）。
• しかし、一般的には計算が困難。→MCMC法などで帰納的に事後確率（最大）を求める。

マルコフ連鎖モンテカルロ法（MCMC法）

• モンテカルロ法：乱数を使ったシミュレーション手法
• マルコフ連鎖：現在の状態が、前時刻の状態のみに依存するモデル→メトロポリス法などの利用

1. θ0を仮定
2. 乱数εを加え、θ1を作成。
3. 確率密度の比ｒを計算r=p(D|θ1)p(θ1)/ p(D|θ0)p(θ0)
   ノイズを正規分布と仮定した場合、p(D|θ) ∝ exp(-NE(θ)/σ2)
   あるθ0、θ1における exp(-NE(θ)/σ2)p(θ) を求める
   ←obs=D を利用 ※点での推定
   
4. 確率密度の比ｒを評価
   r>1でθ1を採用
   r<=1かつ乱数U（1,0）＜ｒでθ1を採用、U>=ｒでx0を採用
5. １～４の繰り返し←exp(-NE(θ)/σ2)p(θ)の分布形状を推定

 

他手法との違い

• 最小二乗法:残差の平方和 E(θ) を最小化
  
• 最尤推定:a,b,σについて導関数=0の位置（パラメータの値）を求める(global minimum とは限らない)
  
• ベイズ推定:パラメータを確率的に扱う→分布(形)を求める 

******************************************

20201129

MAP推定：パラメータ値の同定のみの場合（パラメータの分布を求めない）

• 作成した予測分布モデル＋観測データ（obs）に対し、MAP（maximum a posterior :最大事後確率）推定でパラメータを同定
• pyro.condition で同定パラメータをモデルに反映

 

20201202

変分推論：分布も推定

• svi = SVI(model, guide, optimizer, loss=Trace_ELBO())でモデル構造を作成
  
• step method でデータDをモデルに渡して調整←obsの役目
• predictive で予測分布を作成（モンテカルロ近似）
  

概ね、obs の役割について理解しました。
数式は理解できていません。残念ながら後回しです。

 

20201204

数式もOK（おそらく）。

これで MCMC（メトロポリス法利用）は理解できたと思います。

PyMC

PyMCを使いたくて苦闘。

MCMCを理解する過程で必要となったのですが、これが動きません。

PyMC3 が動いたと教えていただいたライブラリVer.の組み合わせを真似し、仮想環境を構築したのですが、ダメでした。PIP を混ぜてPyMC4にもしてみましたが、ダメ。Anaconda を最新版に入れなおしてもダメ、gcc 関連のパスを通してもダメ。お手上げでした。

よく引っかかったのは Theano。これ、Deep Learning にも流用できる数値計算ライブラリとして有名だったのですが、既に開発が終わって時が経っています。
ちなみに、PyMC4 は TF を利用。TFも２年ほど前はよく利用されていましたが、今は PyTorch に首位の座を奪われているように感じます。

バックエンドが変更になるたびにコードが動かなくなるのは避けたいですね。Docker しかないかな？

統計モデリング

先日、以下の図書をベースにした統計モデリングにかかわる話を聞きました。

樋口「予測に生かす統計モデリングの基本」

以前読んだ「地下水モデル」にデータ同化が記載されており、どのようなものかを知りたいと購入していた図書でした。が、意味が分からず寝かせていました。

その後、逆解析に粒子フィルタを利用された方の話を聞いたこともあり、この機会に引っ張り出して読み直しました。

ま、結局、最後まで読んでも具体的にどのように役立つのか、どのようにシミュレーションに組み込めば良いかは、相変わらず理解できませんでした。

基本コードが手元にあるので、いずれ迫られたら読んでみましょう。

以下、備忘録です。6章まで。

*********************************

加法定理の一般化、同時確率から確率変数を1つ取り除く作業が周辺化

→周辺確率を求める作業

P（A1）が周辺確率。

P（A|B)が「Bが所与のもとでのAの」条件付確率

周辺化

P（A)=∑B P(A,B)

書き下すと（加法定理）

P(A1)=P（A1,B1)+P（A1,B0)

P(A0)=P(A0,B1)+P（A0,B0)

全体：100、A1：30、B1：60、A1∩B1：10の時、

P(A1)=10/100+20/100=30/100

乗法定理

P(A,B)=P(A|B)p(B)=P(B|A)p(A)

書き下すと

P（A1,B1)=P(A1|B1)p(B1)=P（B1,A1)P(A1)

10/100=10/60*60/100=10/30*30/100

ベイズの定理（乗法定理を周辺確率で割る）

P(A|B)=p(B|A)p(A)/P(B)

=p(B|A)p(A)/∑A P(B,A)∵周辺化

=p(B|A)p(A)/∑A(P(B|A)p(A))∵乗法定理

*********************************

状態空間モデル：システムモデルと観測モデルの連立

線形：ガウス状態空間モデル

システムモデル（計算更新時のシステム誤差 vt を含む）

xt=Ftxt-1+Gtvt　 vt～N(0,Qt)

観測モデル（観測誤差wtを含む）

yt=Htxt+wt　 wt～N(0,Rt)

非線形：非線形・非ガウス状態空間モデル（データ同化で利用）

xt=ft(xt-1,vt)　 vt～p(v|θsys)

yt=ht(xt,wt) 　wt～p(w|θobs)

θは未知→最尤法で推定

マルコフ性1：xt～p(xt|xt-1)

マルコフ性2：yt～p(yt|xt)

*********************************

予測分布:p(xt|y1:t-1)

フィルタ分布:p(xt|y1:t) 

平滑化分布:p(xt|y1:T)

情報処理での「フィルタ」は信号処理のフィルタとは別

情報処理での「平滑化」はデータ解析の平滑化とは別。回顧による知識発見が目的。工学系の応用問題は、リアルタイムに適切な処理を施すことが重要で、平滑化分布から得る意味があまりない。

固定ラグ平滑化

固定点平滑化：データ取得毎に初期分布の改善（データ同化で利用）

※固定点平滑化は地下水シミュよりも、力学シミュで役立ちそうです。が、具体的な処理方法が思い浮かぶまで理解できていません。うーん。

仮説検定 用語と手順

上田拓治「44の例題で学ぶ 統計的 検定と推定の解き方」より（一部、変更しています）。

両側検定
・ＡとＢに違いがあるか？など。
片側検定
・ＡはＢより大きいか？多いか？など。


Ｈ0：帰無仮説
・差はない
・検定統計量の実現値pが棄却域内（有意水準（例えば5%）未満）でＨ0棄却、棄却域外でＨ0採択。
※H0棄却時の誤り（第1種の誤り）を犯す確率は5%未満（有意水準）。
・Ｈ0を棄却するとき（Ｈ1を採択するとき）は「差がある」。
・Ｈ0を採択するときは「差があるとは言えない」という表現にとどめる。どちらが正しいかは不明。
Ｈ1：対立仮設
・差がある


パラメトリック検定
・母集団の分布型を決める母数（平均、分散、標準偏差など）について仮定を設ける検定法
・平均値の差、分散・分散比の検定など。
ノンパラメトリック検定
・仮定を置かない検定法
・比率（binary）の差、適合度、独立性、順位の検定など。


仮説検定の手順
1.帰無仮説Ｈ0の設定（差はない、棄却が前提）
2.対立仮設Ｈ1の設定（差がある、採択が前提）
3.条件（基本統計量、サンプルサイズ、検定目的）の確認
4.条件にあった検定手法の選定
5.検定統計量T、検定統計量の実現値pと有意水準α、棄却域Rの比較
6.結論

母集団の分布

上田拓治「44の例題で学ぶ 統計的 検定と推定の解き方」より

まずは分布の関係から。
それぞれ何かの頭文字を取っていると思うのですが、まったくわかりません。Fだけ大文字の理由も。

z分布（標準正規分布）
・平均0
・分散1
・標準化に使用。
・σが既知、もしくは未知でも大標本（100以上）の時の母平均の検定や推定に用いられる。
※ラプラスの中心極限定理：正規分布にかかわらず、nが大きな場合はほぼ正規分布とみなしてよいという定理。
※大数の法則：多くの事象が集まるとそこに法則性が見えてくることを意味する。
→標本平均と母平均の差がεを超える確率はnを十分に大きくとると0に近づく。

ｔ分布
・平均0
・分散は自由度ｆが小さい場合にｚ分布に比べ広がる。サンプル数１００以上でｚ分布とほぼ一致（無限大で一致）
・σが未知かつ小標本（100未満）の検定や推定に用いられる。

χ2分布
・z分布に従う確率変数の2乗
・母分散の検定や推定、適合性や独立性の検定に用いられる。

Ｆ分布
・2種類の分散比から成り立つ分布
・母分散の比の検定や推定に応用される。
・Ｆ値は1以上なので、分散の大きい方を分子として検定統計量を計算する。
・2群の差異を検定したい場合、結果が等分散の場合はｔ検定、非等分散の場合はウェルチｔ検定に続く。

判別分析

他支店の方から、SPSS を買ったと伺いました。

判別分析などを実施する目的で購入されたとのこと。
描画は遅いものの計算部分は高速で、数百万行のデータでも容易に扱えるそうです。さすがに商用ソフトは最近のハードに最適化されているのでしょう。逆に、数百万行を高速に扱えないソフトは売れないでしょうね。

聞いたところ、数量化2類が入っていないとのことでした。
調べてみると、ずいぶん前にオプション販売は終了しています。ま、判別分析が入っているのなら前処理としてダミーデータを作成（ integer encoding ）しておけば良いだけなので、オプションとしての需要がなくなったのでしょう。あるいは、GBM など性能の良い判別機が出ていますので、アルゴリズム自体をあえて利用する機会が減ったのかも。

最初、判別分析は SVM など Classification Algorithm の総称を指しているのかと思って聞いていたのですが、異なっていたようです。scikit-learn で言うところの LDA, QDA アルゴリズムを指していたようです。見る限り、分布のピークを分離しているようですね（利用を迫られた段階で調べます）。
https://scikit-learn.org/stable/modules/lda_qda.html

GUI で高速というのはありがたいでしょう。備えていなくても迫られた際にすぐに答えが出せる点で実務向き。
個人的には Python で十分なのですが、どのような機能があるのか知りたいところです。

RFEM & FORM

昨日と同じ文献です。
Man-Yu Wanga et al.(2020) Probabilistic stability analyses of multi-stage soil slopes by bivariate random fields and finite element methods

RFEM はガウス分布に従ったモンテカルロシミュレーションによって不均質性を表現し、FEM や FDM モデルとした後に SSR 法により安全率を求めるという手法（空間にかかわるパラメータの扱い方がまだ理解できていません）。
知らなかったのですが、比較的多くの事例がありますね。

RFEM に至る流れは、おそらく以下の通り。
3つのレベルで確率を考慮する手法が紹介されています。FORM は通過点でしょう。
https://ieaghg.org/docs/General_Docs/5th_Risk_Assess/VaughanGriffithsieaghg_secured.pdf

FORM において標準正規分布を盲目的に当てはめるのが乱暴で、全体の分布を把握しないといけないというのが直感的にわかる図、言葉も含まれています。
The "safer" slope has a higher "probability of failure"!

FORM は国総研資料第379号でもを扱われています。正規分布でない2次元の確率分布を正規分布で近似し、性能関数に直行する1軸で評価する手順を示されています。取り扱いを簡易にしたいということなのでしょう。ちなみに、扱われた例では second-order でも精度は変わらないという結果でした。
http://www.nilim.go.jp/lab/bcg/siryou/tnn/tnn0379.htm


10年以上前の資料まで追いかけてみましたが、知らないことが多くありました。確率論を避けていたからかもしれません。
ありがたいことに10年以上前では取り扱いが難しい不規則な2次元確率密度分布でも、今ではPython のおかげで誰でも容易に取り扱うことができるようになっています。2変量を標準化しなくても、そのままカーネル密度推定で表現できますし、ある領域化の確率密度を積分することが可能です（手を動かしたことはないですが）。
試してみましょうか。

安定計算と崩壊確率

災害や崩壊にかかわる危険度として相対値なり確率なりの導入を試みる事例を時折目にします。

斜面安定計算でも確率を扱えます。

安定計算に強くかかわるキーパラメータに ｃ・φ があります。この2軸上で確率密度分布を表現するには多くの試験が必要になります。
円弧すべりを考える場合に1軸＋簡易CUを多数実施し、ｃの1次元分布を把握する。それをカーネル密度分布などで連続関数にしてFs=1.0未満の累積確率を求める。このような手順を踏めば、ある条件下での崩壊確率は得られるでしょう。
港湾分野では他分野よりも対象とする土質が比較的均質で室内試験を多く実施する傾向にあるため、上記のような確率の考え方を導入することは容易です（平成19年版設計事例集でも FORM として触れられていました。最新版では消えましたが）。

これが斜面だと難しい。基本は水も土質も均質ではない条件下で、円弧にならないすべり面形状が前提のため、多数の観測や室内試験が必要となります。この点で難色を示される方が多くいらっしゃるでしょう。仮に、試験結果が多数得られたとしても、基準にない確率の考え方を受け入れてもらうことへの抵抗が予想されます。基準に載っており、簡単で、実績のある手法以外を導入し難い気持ちはある程度理解できます。

多量のN値があれば、その確率密度推定結果をφの分布に変換することは、場合によってはありかもしれません。が、数がないので標準正規分布といった仮定を持ち込むのは、乱暴でしょう（数値実験なら良いのですが）。ちなみに、数値実験では標準正規分布ではなく、対数正規分布が多く使われているようです（∵non negative）。例えば以下の文献。
Man-Yu Wanga et al.(2020) Probabilistic stability analyses of multi-stage soil slopes by bivariate random fields and finite element methods

仮に2次元の確率密度分布が得られたところで、それを扱うにも精度の問題が生じます。安全率では計画安全率の0.1%の不足は認められません。確率密度の積分で、その精度が必要とは思えないのですが、基準などで「確率〇％以上を要求」などと決まってしまえば、安全率同様に0.1％の誤差は許されなくなります。計算時間がかかっても、精度の良い積分方法を選択せざるを得ないでしょう。


確率を論じるにはデータが必要になり、それには費用と時間を要します（リスク発現時の対応費用を考えると、微々たる金額ですが）。
こうなると純粋な技術論ではなくマネージメントの領域になります。その判断を経験に依存する分野では、確率論との併用化を進めることは難しいでしょう。

確率を導入する流れに反対ではありません。が、実務上、クリアすべき問題点は多くあるように感じます。

統計学の今後

統計学ではデータセットを training（＋ varidation）, test に区分することはなかったでしょう（況や、CV をや）。

このあたり、機械学習では標準。結果を出すために必要な手順です。
重回帰でもできるのですが、どの統計の図書を見ても理論ばかりで、この流れには触れられていませんでした。予測に対して理論よりも結果（性能）重視の立場が機械学習、理論武装（仕様）重視が統計というような感じを受けます。
単にプログラミングしやすくなったというのも後発に影響しているでしょう。ハードルが低くなり、統計でも CV 等が当たり前になってくるのでしょうか。今、統計のプロはどうされているのでしょう。

機械学習における特徴量選択手法の呼び方 wrapper method は明らかにプログラミングの立場からですね。もちろん、手順に則した呼び方もあります。 Sequential Feature Selection (Sequential Forward Selection , Sequential Backward Elimination, Sequential Forward Floating Selection (SFFS), Sequential Backward Floating Selection (SBFS)), Exhaustive Feature Selection 等。
手元の統計学の図書では SFFS を stepwise と記載されています。手順に関わる呼び方はある程度共通なのかもしれません。が、統計で wrapper method と呼ぶ発想はないでしょう。ある程度の住み分けができているのでしょうね。

どちらの立場からでも良いと思います。両方知っているのが理想でしょう。
PCA、ラフ集合、回帰など、統計の世界にも入り始めましたので、今後も機会があれば学んでいきましょう。