T CrB Scope

Probability transition of the next outburst of T CrB during 2026

概要

本実装では AAVSO (アメリカ変光星観測者協会) の長期観測アーカイブ data/aavso_1866_2026.txt を用いて再帰新星 T CrB (かんむり座 T 星) の次回新生爆発の発生確率を 2026 年内の日次分布として近似している。ここで出している値は物理シミュレーションではなく、1945-1946 年の前駆期と現代の B, V, TG, B-V 系列を比較した統計的・解析学的モデルである。

但し現行版ではこれに加えて ASAS-SN Sky Patrol の最新 light curve が取得できた場合に限り極く弱い nowcast 補正を掛ける。補正は主モデルを書き換えるものではなく、archive-only の日次確率に対して数 % 程度の bounded multiplier を掛ける補助層である。

成果物

• 確率プロット: assets/prob_graph.png
• 生成スクリプト: scripts/plot_2026.py
• 基礎解析: src/tcrb/analysis/eruption_study.py
• live 補助取得: src/tcrb/service/monitoring.py

実行するには:

uv sync
source .venv/bin/activate
python scripts/plot_2026.py

データ前処理

入力は AAVSO のテキストアーカイブであり、ローダーは以下を行う。

• 対象名の正規化
• JD -> UTC 変換
• Magnitude の数値化
• T CrB, 000-BBW-825, 1555+26 等のエイリアス統合

各バンド時系列は日単位 median に纏めた上で全日付格子に再標本化し、Savitzky-Golay 平滑化を施す。平滑系列を $x_t$ とする。

平滑系列は各バンドごとに基準期 median を引いて、正の異常量のみを使う。

$$a_t = \max(x_t - b, 0)$$

ここで $b$ は基準期 median である。現代系列では基準期を主に 2021-01-01 から 2022-12-31 に置いている。

アナログ推定

最終確率は、1945-1946 visual 前駆期を参照系列として、3 種類のアナログ推定を作り合成する。

1. Shape Match

現代系列の直近 180 日を正規化し、1945-1946 の参照系列の全ての 180 日窓と比較する。現代系列を $c_i$、歴史窓を $h_i$ とすると、評価量は

$$\mathrm{RMSE} = \sqrt{\frac{1}{n} \sum_i (c_i - h_i)^2}$$

である。$\mathrm{RMSE}$ 最小の歴史窓終端を $t_h$、1946 爆発日を $t_*$ とすると、推定日 $\hat{t}$ は

$$\hat{t} = t_{\mathrm{now}} + (t_* - t_h)$$

で与える。

2. Phase Ratio

現代異常量の現在値とピーク値の比

$$r = \frac{a_{\mathrm{now}}}{\max_t a_t}$$

を求め、1945-1946 の減衰枝で同じ振幅比に最も近い日を探す。そこから爆発日までの残り日数を現代へ写像する。

3. State Match

位相状態を 3 変数ベクトルで表す。

$$s = (f, g, c)$$

各成分は

• $f$: 異常量のピーク比 $f = \frac{a_{\mathrm{now}}}{\max_t a_t}$
• $g$: 直近傾きのピーク規格化 $g = \frac{\mathrm{mean},\Delta a_t}{\max_t a_t}$
• $c$: 累積異常量の進捗 $c = \frac{\sum_{u \le t} a_u}{\max_t \sum_{u \le t} a_u}$

歴史系列の各時点にも同じベクトル $s_h(t)$ を作り、

$$d(t) = \sqrt{\sum_k w_k (s_k - s_{h,k}(t))^2}$$

を最小化する。ここで重みは実装上 $w = (1.2, 1.0, 0.8)$ である。最小距離の歴史時点を現在の位相対応点と見做し、そこから爆発日までの残り日数を写像する。

手法較正

各手法は 1945-07-01 から 1946-01-15 までの rolling cutoff で自己バックテストし、各カットオフ時点から 1946-02-09 をどれだけ誤差なく予測できるかを測る。

各手法の平均絶対誤差を $\mathrm{MAE}_m$ とすると、素の手法スコアは

$$q_m = \frac{1}{\mathrm{MAE}_m + 30}$$

で定義する。$\mathrm{phase}$ は同じ前駆系列に対する自己一致が強く、過学習しやすいので penalty を掛ける。最終重みは

$$\alpha_m = \frac{q_m}{\sum_j q_j}$$

で正規化する。

Latent Phase Offset

$B/TG/V$ はより achromatic な前駆相、$B-V$ は chromatic な遅れとして分離する。

まず $\mathrm{state}$ 推定のうち $B$, $TG$, $V$ を achromatic 群、$B-V$ を chromatic 群とし、achromatic 側の重み付き中心日を

$$t_A = \frac{\sum_i w_i t_i}{\sum_i w_i}$$

chromatic 側の中心日を $t_C$ とする。両者の差

$$\Delta = t_C - t_A$$

を latent phase offset と見做す。

更に現在の $B-V$ 赤化速度と $B/TG$ 回復速度を比較し、offset をどの程度保つかを動的に決める。30 日・90 日の平均傾きを混合して、

$$r_{\mathrm{red}} = 0.65,r_{30}^{B-V} + 0.35,r_{90}^{B-V}$$

$$r_{\mathrm{rec}} = \frac{\max(-r_{30}^{B}, 0) + \max(-r_{30}^{TG}, 0)}{2}$$

に 90 日成分も同様に混合したものを作る。これにより、色の遅れがまだ優勢か、あるいは achromatic 側が先に回復しているかを判断する。

Achromatic Trigger と Chromatic Drag

最終モデルは単一 latent date を直接使わず、2 変数のハザードで合成する。

Achromatic Trigger

$B$, $TG$, $V$ の $\mathrm{shape}$, $\mathrm{phase}$, $\mathrm{state}$ シナリオを achromatic 群とし、重み付き中心日 $\mu_T$、重み付き spread $\sigma_T$、総重み $W_T$ を作る。trigger hazard は

$$h_T(t) = W_T,\phi(t; \mu_T, \sigma_T)$$

で表す。ここで $\phi$ は正規密度である。

Chromatic Drag

$B-V$ の $\mathrm{shape}$, $\mathrm{phase}$, $\mathrm{state}$ シナリオから drag 中心日 $\mu_D$、spread $\sigma_D$、総重み $W_D$ を作る。

但し drag 自体は単一の CDF ではなく、$B-V$ 異常量の将来投影から時変に作る。現在の色異常量を $A_0$、drag center までの増加率を $r_{\mathrm{red}}$ とすると、drag center までは

$$A(t) = A_0 + r_{\mathrm{red}} t$$

とする。drag center 後は 2 段の piecewise decay に切り替える。遅い減衰率を $\lambda_s$、速い減衰率を $\lambda_f$、肩の長さを $\tau_s$ とすると、

$$A(t) = A_{\mathrm{peak}} - \lambda_s (t - \mu_D) \quad \text{for } 0 \le t - \mu_D \le \tau_s$$

$$A(t) = A_{\mathrm{shoulder}} - \lambda_f (t - \mu_D - \tau_s) \quad \text{for } t - \mu_D > \tau_s$$

である。

この projected anomaly から drag gate を

$$g_D(t) = \frac{1}{1 + A(t) / s_D}$$

とし、更に完全にはゼロにならない floor を設けて

$$\tilde{g}_D(t) = g_0 + (1 - g_0) g_D(t)$$

とする。$s_D$ は drag 強度に応じた anomaly scale、$g_0$ は drag floor である。

Live Nowcast 補正

scripts/plot_2026.py は、src/tcrb/service/monitoring.py の build_monitoring_snapshot() を通じて ASAS-SN の最新 light curve を取得できた場合のみ、弱い live 補正を作る。現在の実装は Web UI と同じ get_lightcurve/<object_id> 経路を用いる。取得に失敗した場合は自動的に無効化され、完全に archive-only の分布へフォールバックする。

live 系列からは、最近 45 日の平均等級、そこから 180 日以上前の基準 median を引いた anomaly、同じ 45 日窓での線形傾き $s_{\mathrm{live}}$ を作る。異常量を $a_{\mathrm{live}}$、信頼度を $c_{\mathrm{live}}$ とすると、補正振幅 $\epsilon$ は

$$\epsilon = \mathrm{clip}\left(0.03 + 0.04 c_{\mathrm{live}} + 0.06,\min\left(\frac{|a_{\mathrm{live}}|}{0.6}, 1\right),, 0.03,, 0.07\right)$$

で与える。従って補正倍率は常に

$$1 - \epsilon \le M_{\mathrm{live}}(t) \le 1 + \epsilon$$

を満たし、主モデルを大きくは動かさない。

補正の向きは anomaly と傾きから

$$d_{\mathrm{live}} = 0.65,\mathrm{clip}(a_{\mathrm{live}} / 0.5, -1, 1) + 0.35,\mathrm{clip}(s_{\mathrm{live}} / 0.01, -1, 1)$$

で定める。これを未来 45 日付近を中心にしたガウス型プロファイル $k(t)$ に乗せて

$$M_{\mathrm{live}}(t) = \mathrm{clip}(1 + \epsilon, d_{\mathrm{live}}, k(t), 1 - \epsilon, 1 + \epsilon)$$

を作る。

ここで重要なのは、ASAS-SN は prior でも主成分でもなく、あくまで nowcast である点である。よって light curve が得られても、補正は最大で数 % に制限される。

最終日次確率

シナリオ混合密度を

$$m(t) = \sum_j w_j,\phi(t; \mu_j, \sigma_j)$$

trigger hazard を $h_T(t)$、drag gate を $\tilde{g}_D(t)$ とすると、最終密度は

$$p_{\mathrm{raw}}(t) = 0.35,m(t) + 0.65,h_T(t),\tilde{g}_D(t)$$

で与える。ASAS-SN live 補正が有効な場合は更に

$$p_{\mathrm{live}}(t) = p_{\mathrm{raw}}(t),M_{\mathrm{live}}(t)$$

を使い、無効な場合は $p_{\mathrm{live}}(t) = p_{\mathrm{raw}}(t)$ とする。最後に

$$p(t) = \frac{p_{\mathrm{live}}(t)}{\sum_u p_{\mathrm{live}}(u)}$$

として 2026 年内の離散日次確率へ正規化する。

現在の assets/prob_graph.png では、破線が archive-only、実線が live-adjusted の日次確率を表す。

現在の解釈

現時点の実装では、概ね次の構図になっている。

• achromatic trigger 中心: 2026-05-20
• chromatic drag 中心: 2026-07-07
• latent offset: 約 75 日
• 最頻域: 2026-05-17 から 2026-05-20

従って、モデル上は achromatic にはかなり成熟した前駆状態にあり、しかし $B-V$ の色異常がまだ drag として残っているため、確率は 5 月後半へ押し広げられていると解釈している。

2026-03-30 時点の実行では、ASAS-SN の get_lightcurve/137439769436 から 1280 点の light curve 自体は取得できるが、系列の最新点は 2024-05-25 であり、2026 年の予測窓に対して古すぎる。そのため実装では stale な live signal と判定し、live 補正を自動で無効化する。現在の判定閾値は予測窓開始の 120 日より前であり、従ってこの場合のプロットは archive-only と live-adjusted が一致する。取得系列の利用可能フィルタは V と g である。

ライセンス

このリポジトリにおけるソースコードは MIT License に従います。

謝辞 / Acknowledgments

This project uses observational data and related services from the AAVSO International Database (AID), the AAVSO Variable Star Index (VSX), the AAVSO Variable Star Plotter (VSP), ASAS-SN public data services, and the SIMBAD database.

I gratefully acknowledge the contributions of the AAVSO observer community, whose photometric data and metadata resources were used in this project and made available through the AAVSO’s scientific archives.

This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.

This repository is an independent, non-peer-reviewed project created for personal research and software development purposes.