ガウス過程回帰で電気抵抗の温度依存性を予測できるかやってみた

 【イントロ】

実験て楽しいですよね。

特に初めて作った物質の物性を初めて測定する瞬間、

世界で誰も見たことがない物性が測定が進むにつれて明らかになっていく時間、

最高のときめきです。

ときめきメモリアル forever with you エモーショナル©

ところがどっこい。順調に測定がすすめばよいですが、

途中で相転移なんて起こった日にはたまらない。

もっと温度変化を細かくとっておけばよかった！

と後悔するころには実験マシンタイムは終わってしまいます。

え？1回目は試し測定で荒く測定して、2回目で細かく測ればいいだけ？

はぁ・・・

鬼滅の刃©

そんな悲劇を回避するために、測定済みのデータ点からこれから測定するデータが取りうる値の幅を予測することは可能でしょうか？

そうした手法の1つがガウス過程回帰です。

本記事では、ガウス過程回帰を用いて実験データを予測することが可能か、試してみました。

【手法】

本記事では、予測する実験データとして、電気抵抗の温度依存性を例に着目しました。

電気抵抗の温度依存性の例として、以下３つについて、分析を行いました。

1. 一般的な金属
2. 半導体
3. 超伝導体

ガウス過程回帰(Gaussian Process Regression)については、持橋他「ガウス過程と機械学習 (機械学習プロフェッショナルシリーズ) 」を参考にしました。

一般的な線形回帰では、係数を測定点から決定しますが、ガウス過程回帰では、測定点自体が特定の確率分布から得られると仮定して、予測の分散付きで予測値を算出することが可能となります。このとき、測定点から近い予測点は分散が小さく、遠い点は分散が大きく予測されます。測定点と予測点の近さを表す関数はカーネルと呼ばれますが、カーネル関数を調整し、さらにカーネルのパラメータを測定データを使って決定することで、様々な関数を学習し予測することが可能になります。

今回はカーネルは標準的なRBF(Radial Basis Function)を利用しています。

コーディングについては、ChatGPTとClaudeを利用しました。

利用したノートブックのHTMLを記事の末尾に記載します。

【結果】

まず一般的な金属の例です。

一般的な金属では電子とフォノンによる散乱が電気抵抗の主な要因です。

電子フォノン散乱による電気抵抗ρの温度依存性は、高温ではT-linear、低温ではT^5となることが知られています（参考文献）

今回は、温度依存性は単純にT^5とし、残留抵抗はゼロとしています。つまり、

　ρ(T) ∝ T^5

ということです。

下図では、測定を300Kから初めて、200K、150K、100K、50Kまで測定したときの測定点と、GPRによる予測点および標準偏差をプロットしています。

200Kまで測定した時点では、100K以下で抵抗が再上昇し、ゼロ温度の予測値はかなり大きな標準偏差を持っていることがわかります。

一方の測定点が50Kに向かって増えるにつれてこの標準偏差は小さくなり、50Kまで測定した時点では、ほぼ理論値の0に収束していることがわかります。

200Kまで測定した時点では低温でアップターンする可能性もあることがわかれば、心構えができてよいですね。

図１，金属的温度依存性のGPRによる予測

次に半導体の場合です。

半導体の電気抵抗の温度依存性です。

単純には半導体の電気抵抗の温度依存性はρ∝T^(-1.5)exp(1/T)となることが知られています（参考文献）

今回は、フォノンによる散乱も加えて、

　ρ(T) ∝ T^(-1.5)exp(1/T) ＋ T^5

の形式で理論値を設定しています。

高温では温度低下に伴い電気抵抗が減少し、低温では半導体的性質にともない、抵抗が上昇する振る舞いを想定しています。

下図のように200Kまで測定した時点では、低温のアップターンを予測できていないことがわかります。さらに測定点を50Kまで下げても、予測の標準偏差は小さくなっていますが、低温のアップターンは正確に予測できていないことがわかります。

RBFカーネルでは、指数関数的に変化するふるまいは正確に予測するのが難しいようです。

図２，半導体的温度依存性のGPRによる予測

最後は超伝導体です。

今回の超伝導体は、Tc=95Kでゼロ抵抗に変化することを想定しています。

Bi2212とかYBCO的なやつですね。

また、超伝導体では超伝導転移温度以上から抵抗に変化が生じる超伝導ゆらぎという現象があるため、その振る舞いも理論値に含めています。

まとめると、高温ではT-linear、転移温度付近では超伝導ゆらぎを考慮して、転移温度以下ではゼロ抵抗としています。つまり、

　ρ(T) ∝ T (if T>Tc), 0 (T<Tc)

ということです。

下図の200Kまでの測定では、測定点ではほぼ直線にみえています。GPRの予測もほぼ直線ですが、若干下振れするような予測をしています。これが、150Kまでの測定点を増やすと、転移温度は予測できていませんが、超伝導ゆらぎによる抵抗の下振れは正しく予測できています。更に転移温度直上の100Kまでの測定点を増やすと、ほぼTcと一致する温度でゼロ抵抗をよぎる予測値を出力していることがわかります。一方で低温部分では抵抗有限値に増えると予想しており、物理的に正しい値を予測できていません。

いずれにせよ、比較的高温から超伝導ゆらぎに伴う抵抗の下振れを予測できているように見えるのは、興味深いことです。

【サマリ】

今回の記事では、ガウス過程回帰で電気抵抗の温度依存性を予測してみることを試みました。シンプルなGPRとして、RBFを使ったガウス過程回帰を実施しました。

金属的な電気抵抗の温度依存性や、超伝導ゆらぎをもつ超伝導体の抵抗の温度依存性は比較的よく予測することができ、測定中に予測を併用すれば、興味深い温度で温度測定幅を細くすることができそうです。

一方で、半導体な指数関数的に急激に変化する温度依存性はうまく予測できない結果となりました。利用したカーネルをRBFから変更することで予測が可能になる可能性はあり、将来の課題といえそうです。

簡単な例でのガウス過程回帰の実践でしたが、予測値の標準偏差を同時に可視化できるのは、単純な線形回帰では難しい手法のため、非常に有力な手法に感じられます。また、関数系が不明な物性に対しても予測値を返せる点は魅力的です。ただし、背後の物理的機構の解明にはそれだけでは不十分といえるため、利用する場面は検討が必要そうです。場面に応じた手法の選定が必要なのは世の常です。

今後さらにやってみたいことは以下のような分析です。

• 相転移を挟んだ予測はできないか
• 転移点温度以上から相転移を捉えられないか
• 転移温度以上のゆらぎの情報を捉えて、なんらかの物理パラメータを抽出できないか
• 常伝導体と思われているが、実は超極低温で超伝導体になる物質を超伝導ゆらぎを捉えてスクリーニングできないか
• 量子振動などの周期的変化をするデータにも適用できるか（カーネルを調整すればできる想定）
• ガウス過程回帰の課題は、計算量がデータ点NのN^3で増える点。計算を省力化する方法を実装してみたい

しかし、今の御時世は、アイデアがあればAIでコード実装が可能になって便利なものです。出てきた結果の正しさを検証するだけの知識が必要ですが、たたき台があるのと無いのとでは雲泥の差、学習速度を加速させていきたいものです。

【参考文献】

持橋他,「ガウス過程と機械学習 (機械学習プロフェッショナルシリーズ) 」

井野明洋, 固体物理学 I 講義ノート（PDFリンク）

格子振動による散乱（Webページ、20251014確認）

赤穂 昭太郎, ガウス過程回帰の基礎

Jinglai Li, Hongqiao Wang, Gaussian Processes Regression for Uncertainty Quantification: An Introductory Tutorial arXiv:2502.03090 [stat.CO]

Jie Wang, An Intuitive Tutorial to Gaussian Process Regression, arXiv:2009.10862 [stat.ML]

金子弘昌, Ｐｙｔｈｏｎで学ぶ実験計画法入門　ベイズ最適化によるデータ解析 (ＫＳ情報科学専門書)

[参考コード]

In [ ]:

import numpy as np
import matplotlib.pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C
import pandas as pd
from joblib import Parallel, delayed

In [ ]:

np.random.seed(0)  # 再現性のため
noise_std = 0   # 標準偏差（ノイズレベル）を設定

In [ ]:

# 最高温度
T_max = 300
# 最低温度
T_min = 1

# 電気抵抗の温度依存性
# mode = 'phonon'
# mode = 'linear'
# mode = 'sc-fluctuation'
mode = 'semiconductor'

In [ ]:

# 各条件で測定済み温度範囲を変える
measured_ranges = [
    (200, 300),  # 左上
    (150, 300),  # 右上
    (100, 300),  # 左下
    (50,  300),  # 右下
]
titles = [
    "Measured: 200–300 K",
    "Measured: 150–300 K",
    "Measured: 100–300 K",
    "Measured:  50–300 K",
]

In [ ]:

# --------------------------
# 1. 真の抵抗関数（例：金属）
# --------------------------
def true_resistivity(T, mode, noise_std=noise_std):
  if mode == 'phonon':

    R0 = 1  # 残留抵抗
    a = 1.0e-1  # フィッティング定数

    return R0 + a * T**5 + np.random.normal(0, noise_std)

  if mode == 'linear':
    a =1.0

    return a * T

  if mode == 'semiconductor':
    Eg = 0.01  # eV程度（大きいほど絶縁性が強い）
    kB = 8.617e-5  # eV/K
    a = 1
    b = 1.0e-12

    return b*T**5 + a * T**(-1.5)*np.exp(Eg / (2 * kB * T)) + np.random.normal(0, noise_std)

  if mode == 'sc-fluctuation':
    """
    完全に連続的な超伝導転移とゆらぎを含む電気抵抗関数

    Parameters:
    -----------
    T: 温度 [K] (float or ndarray)
    Tc: 超伝導転移温度 [K]
    R0, a: 高温での金属的抵抗 R0 + a*T のパラメータ
    rho_n: 転移直上での正常状態抵抗率
    fluc_strength: ゆらぎによる抵抗低下の最大強度 (0-1)
    fluc_decay_length: ゆらぎ効果の減衰長 [K]
    transition_width: 転移の幅パラメータ [K]
    noise_std: 観測ノイズ標準偏差

    Returns:
    --------
    R_eff: 有効電気抵抗（完全に連続的）
    """
    Tc=95.0
    R0=1.0
    a=1.0
    rho_n=10.0 # 転移直上での正常状態抵抗率
    fluc_strength=0.3  # ゆらぎの強度 (0-1)
    fluc_decay_length=15.0  # ゆらぎの減衰長 [K]
    transition_width=2.0  # 転移幅パラメータ
    noise_std=0.0

    T = np.asarray(T, dtype=float)

    # === 完全に滑らかな正常状態抵抗 ===
    # 高温: R0 + a*T の線形増加
    # 低温: rho_n 付近での飽和
    # 全域で連続的かつ微分可能

    # スムーズな接続のための重み関数
    transition_scale = 30.0  # 接続の滑らかさを制御
    weight_high_T = 0.5 * (1 + np.tanh((T - Tc) / transition_scale))
    weight_low_T = 1.0 - weight_high_T

    # 高温域の線形抵抗
    R_linear = R0 + a * T

    # 低温域の飽和抵抗（Tcで rho_n になるよう調整）
    R_saturation = rho_n * (1.0 + 0.01 * np.maximum(T - Tc, -50))

    # 重み付き平均で滑らかに接続
    R_normal = weight_high_T * R_linear + weight_low_T * R_saturation

    # === 超伝導転移関数（BCS型） ===
    epsilon = (T - Tc) / transition_width

    # より滑らかな転移関数（エラー関数ベース）
    from scipy import special
    try:
        # erfが利用可能な場合
        sc_suppression = 0.5 * (1 + special.erf(epsilon))
    except:
        # erfが無い場合はtanhで近似
        sc_suppression = 0.5 * (1 + np.tanh(2.5 * epsilon))

    # === 完全に滑らかな超伝導ゆらぎ効果 ===
    fluctuation_factor = 1.0

    if fluc_strength > 0:
        # 温度差
        delta_T = T - Tc

        # Ginzburg-Landau理論ベースの滑らかなゆらぎ関数
        # f(δT) = A * (δT + δ₀)^(-1/2) * exp(-δT/λ) for δT > 0
        #       = 0 for δT ≤ 0

        delta_0 = 0.5  # 発散回避パラメータ
        lambda_fluc = fluc_decay_length

        # T > Tc でのゆらぎ効果
        fluc_amplitude = np.zeros_like(T)
        mask_fluctuation = delta_T > 0

        if np.any(mask_fluctuation):
            dt_pos = delta_T[mask_fluctuation]

            # パワー則成分（GL理論）
            power_term = 1.0 / np.sqrt(dt_pos + delta_0)

            # 指数減衰成分（長距離カットオフ）
            exp_term = np.exp(-dt_pos / lambda_fluc)

            # スムーズな立ち上がり（T = Tcで0から連続的に開始）
            smooth_onset = 1.0 - np.exp(-dt_pos / 0.5)

            # 全体のゆらぎ振幅
            fluc_amplitude[mask_fluctuation] = fluc_strength * power_term * exp_term * smooth_onset

        # ゆらぎによる導電率増加 → 抵抗率減少
        fluctuation_factor = 1.0 / (1.0 + fluc_amplitude)

    # === 最終的な抵抗計算 ===
    # 超伝導転移 × ゆらぎ効果を適用
    R_eff = R_normal * sc_suppression * fluctuation_factor

    # ノイズ付加
    if noise_std > 0:
        R_eff = R_eff + np.random.normal(0, noise_std, size=R_eff.shape)
        R_eff = np.maximum(R_eff, 0.0)

    return R_eff

In [ ]:

# --------------------------
# 2. ガウス過程回帰で予測
# --------------------------
def gpr_predict(T_train, R_train, T_test):
    # カーネルの定義（適宜調整可能）
    kernel = C(1.0, (1e-5, 1e5)) * RBF(length_scale=50.0, length_scale_bounds=(1.0, 1e4))
    gpr = GaussianProcessRegressor(kernel=kernel, alpha=1e-12, normalize_y=True, n_restarts_optimizer=5)

    # 学習
    gpr.fit(T_train.reshape(-1, 1), R_train)

    # 予測
    R_pred, sigma = gpr.predict(T_test.reshape(-1, 1), return_std=True)
    return R_pred, sigma, gpr

In [ ]:

# --------------------------
# 3. 結果を可視化
# --------------------------
def plot_gpr(T_train, R_train, T_test, R_pred, sigma, ax=None, mode='phonon'):
    if ax is None:
        fig, ax = plt.subplots(figsize=(8, 5))

    ax.plot(T_train, R_train, 'ro', label='Measured')
    ax.plot(T_test, R_pred, 'b-', label='GPR Mean')
    ax.fill_between(T_test, R_pred - sigma, R_pred + sigma, alpha=0.3, color='blue', label='±1σ')
    ax.set_xlabel('Temperature (K)')
    ax.set_ylabel('Resistivity (Ohm)')

    # modeが'semiconductor'のとき、y軸をlogスケールにする
    if mode == 'semiconductor':
        ax.set_yscale('log')
        ax.set_ylabel('Resistivity (Ohm, log scale)')

    ax.set_title(f"GPR Prediction ({len(T_train)} points)")
    ax.legend()
    ax.grid(True)

In [ ]:

# ----------------------
# メイン処理
# ----------------------
def main():
   # プロット準備（2×2）
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    axes = axes.ravel()  # 2×2を1次元配列に展開

    # 4条件ループ
    for i, (T_min_measured, T_max_measured) in enumerate(measured_ranges):
        ax = axes[i]

        # 真の関数
        T_full = np.arange(T_min, T_max, 0.5)
        R_full = true_resistivity(T_full, mode)

        # 測定済み領域
        T_measured = np.arange(T_min_measured, T_max_measured, 1)
        R_measured = true_resistivity(T_measured,mode)

        # 予測領域（未測定部分）
        T_predict = np.arange(T_min, T_min_measured, 0.1)

        # ソート
        sort_idx = np.argsort(T_measured)
        T_train = T_measured[sort_idx]
        R_train = R_measured[sort_idx]

        # ---- GPRによる予測 ----
        R_pred, sigma, _ = gpr_predict(T_train, R_train, T_predict)

        # ---- プロット ----
        plot_gpr(T_train, R_train, T_predict, R_pred, sigma, ax=ax, mode=mode)

        # 真の関数も描画
        ax.plot(T_predict, true_resistivity(T_predict, mode), 'k--', label='True function')
        ax.set_title(titles[i])
        ax.legend()

    # 全体レイアウト調整
    plt.tight_layout()
    plt.show()

In [ ]:

if __name__ == '__main__':
    main()

/usr/local/lib/python3.12/dist-packages/sklearn/gaussian_process/_gpr.py:660: ConvergenceWarning: lbfgs failed to converge (status=2):
ABNORMAL: .

Increase the number of iterations (max_iter) or scale the data as shown in:
    https://scikit-learn.org/stable/modules/preprocessing.html
  _check_optimize_result("lbfgs", opt_res)

In [ ]:

## 低温拡大
T_min_measured = 30

In [ ]:

# ----------------------
# メイン処理
# ----------------------
def main_sub():
    # すべての温度と真の抵抗
    T_full = np.arange(T_min, T_max, 1)
    R_full = true_resistivity(T_full, mode)

    # 測定済み温度：150K〜300K
    T_measured = np.arange(T_min_measured, T_max, 1)
    R_measured = true_resistivity(T_measured, mode)

    # 予測したい温度：0K〜150K
    T_predict = np.arange(T_min, T_min_measured, 0.1)

    # 昇順に並び替え（GPRの学習用）
    sort_idx = np.argsort(T_measured)
    T_train = T_measured[sort_idx]
    R_train = R_measured[sort_idx]

    # GPRで予測
    R_pred, sigma, _ = gpr_predict(T_train, R_train, T_predict)

    # プロット
    fig, ax = plt.subplots(figsize=(10, 6))
    plot_gpr(T_train, R_train, T_predict, R_pred, sigma, ax)

    # 真の関数も表示（参考）
    R_true = true_resistivity(T_predict, mode)
    ax.plot(T_predict, R_true, 'k--', label='True function')
    ax.legend()

    ax.set_xlim(0, 50)                # 温度軸を0〜50 Kに制限
    ax.set_ylim(0, np.max(R_measured[T_measured <= 50]) * 1.2)  # y軸を自動で適度に調整
    ax.set_title("GPR Prediction (Low-T region)")
    ax.legend()
    ax.grid(True)

    plt.tight_layout()
    plt.show()

In [ ]:

if __name__ == '__main__':
    main_sub()