ochilab · November 13, 2025 01:28
diff --git a/dunnett_critical_value.py b/dunnett_critical_value.py
 import numpy as np

 def dunnett_critical_value_from_groups(alpha=0.05,
                                       g_total=3,   # ★ 表に書いてある「群数 g（対照＋処理の総数）」を渡す
                                       df=20,
                                       rho=0.5,     # 等nなら 0.5
                                       n_sim=100_000,
                                       random_state=None):
    """
    Dunnett 検定の臨界値 K をモンテカルロで近似する（両側検定想定）by ChatGPT

    Parameters
    ----------
    alpha : float
        有意水準（例: 0.05）
    g_total : int
        群数 g = 対照群 + 処理群の合計数（ダネット表の「群数」）
        例: 対照 + 処理2群 → g_total = 3
    df : int
        誤差自由度（ANOVA の誤差自由度）
    rho : float
        各比較間の相関。繰り返し数がすべて等しい場合は 0.5。
    n_sim : int
        シミュレーション回数（大きいほど精度↑）
    random_state : int or None
        乱数シード

    Returns
    -------
    K : float
        Dunnett の臨界値 K（max |T_i| の (1 - alpha) 分位）
    """
    rng = np.random.default_rng(random_state)

    # ★ 実際の比較本数（処理群の数） = 総群数 - 1（対照を除いた数）
    k = g_total - 1
    if k < 1:
        raise ValueError("g_total は 2 以上（対照 + 少なくとも1処理群）にしてください。")

    # --- 相関行列 Σ の作成（対角1, 非対角rho） ---
    Sigma = np.full((k, k), rho, dtype=float)
    np.fill_diagonal(Sigma, 1.0)
    L = np.linalg.cholesky(Sigma)

    # --- 多変量 t 乱数を n_sim 回生成 ---
    N = rng.normal(size=(n_sim, k))
    Z = N @ L.T
    chi2 = rng.chisquare(df, size=n_sim)
    scale = np.sqrt(df / chi2)[:, None]
    T = Z * scale

    # --- max |T_i| の (1 - alpha) 分位を臨界値とする ---
    max_abs_T = np.max(np.abs(T), axis=1)
    K = np.quantile(max_abs_T, 1 - alpha)
    return K

 # 使い方例
 if __name__ == "__main__":
    alpha = 0.05
    g_total = 4   # 「対照 + 処理3群」の 4 群 → 表の『群数 4』
    df = 36       # 例：浜田先生の資料などで出てくる自由度 36 等

    K = dunnett_critical_value_from_groups(alpha=alpha,
                                           g_total=g_total,
                                           df=df,
                                           rho=0.5,
                                           n_sim=200_000,
                                           random_state=0)
    print(f"Dunnett 臨界値 K (α={alpha}, 群数={g_total}, df={df}) ≒ {K:.4f}")
	import numpy as np

	def dunnett_critical_value_from_groups(alpha=0.05,
	g_total=3, # ★ 表に書いてある「群数 g（対照＋処理の総数）」を渡す
	df=20,
	rho=0.5, # 等nなら 0.5
	n_sim=100_000,
	random_state=None):
	"""
	Dunnett 検定の臨界値 K をモンテカルロで近似する（両側検定想定）by ChatGPT

	Parameters
	----------
	alpha : float
	有意水準（例: 0.05）
	g_total : int
	群数 g = 対照群 + 処理群の合計数（ダネット表の「群数」）
	例: 対照 + 処理2群 → g_total = 3
	df : int
	誤差自由度（ANOVA の誤差自由度）
	rho : float
	各比較間の相関。繰り返し数がすべて等しい場合は 0.5。
	n_sim : int
	シミュレーション回数（大きいほど精度↑）
	random_state : int or None
	乱数シード

	Returns
	-------
	K : float
	Dunnett の臨界値 K（max \|T_i\| の (1 - alpha) 分位）
	"""
	rng = np.random.default_rng(random_state)

	# ★ 実際の比較本数（処理群の数） = 総群数 - 1（対照を除いた数）
	k = g_total - 1
	if k < 1:
	raise ValueError("g_total は 2 以上（対照 + 少なくとも1処理群）にしてください。")

	# --- 相関行列 Σ の作成（対角1, 非対角rho） ---
	Sigma = np.full((k, k), rho, dtype=float)
	np.fill_diagonal(Sigma, 1.0)
	L = np.linalg.cholesky(Sigma)

	# --- 多変量 t 乱数を n_sim 回生成 ---
	N = rng.normal(size=(n_sim, k))
	Z = N @ L.T
	chi2 = rng.chisquare(df, size=n_sim)
	scale = np.sqrt(df / chi2)[:, None]
	T = Z * scale

	# --- max \|T_i\| の (1 - alpha) 分位を臨界値とする ---
	max_abs_T = np.max(np.abs(T), axis=1)
	K = np.quantile(max_abs_T, 1 - alpha)
	return K

	# 使い方例
	if __name__ == "__main__":
	alpha = 0.05
	g_total = 4 # 「対照 + 処理3群」の 4 群 → 表の『群数 4』
	df = 36 # 例：浜田先生の資料などで出てくる自由度 36 等

	K = dunnett_critical_value_from_groups(alpha=alpha,
	g_total=g_total,
	df=df,
	rho=0.5,
	n_sim=200_000,
	random_state=0)
	print(f"Dunnett 臨界値 K (α={alpha}, 群数={g_total}, df={df}) ≒ {K:.4f}")
No results found