import numpy as np

# セル確率 p_ij を与える (例として 3x3 の行列を使用)
p_ij = np.array([[0.2, 0.1, 0.05],
                 [0.15, 0.25, 0.05],
                 [0.05, 0.05, 0.1]])

I = p_ij.shape[0]

# 1. 与えられた p_ij から θ_i, φ_j, λ_ij を計算
p_11 = p_ij[0, 0]  # p_11 はセル (1, 1) の確率
theta_i = np.log(p_ij[1:, 0] / p_11)  # θ_i = log(p_i1 / p_11), i = 2, ..., I
phi_j = np.log(p_ij[0, 1:] / p_11)    # φ_j = log(p_1j / p_11), j = 2, ..., I
lambda_ij = np.zeros((I - 1, I - 1))  # λ_ij を初期化

for i in range(1, I):
    for j in range(1, I):
        lambda_ij[i - 1, j - 1] = np.log(p_ij[i, j] * p_11 / (p_ij[i, 0] * p_ij[0, j]))

# 2. θ_i, φ_j, λ_ij から p_ij を再構成
p_ij_reconstructed = np.zeros_like(p_ij)
p_ij_reconstructed[0, 0] = 1 / (
    1 +
    np.sum(np.exp(theta_i)) +
    np.sum(np.exp(phi_j)) +
    np.sum(np.exp(theta_i[:, None] + phi_j[None, :] + lambda_ij))
)

# 再構成した確率の計算
for i in range(I):
    for j in range(I):
        if i == 0 and j > 0:  # p_1j
            p_ij_reconstructed[i, j] = p_ij_reconstructed[0, 0] * np.exp(phi_j[j - 1])
        elif j == 0 and i > 0:  # p_i1
            p_ij_reconstructed[i, j] = p_ij_reconstructed[0, 0] * np.exp(theta_i[i - 1])
        elif i > 0 and j > 0:  # p_ij
            p_ij_reconstructed[i, j] = p_ij_reconstructed[0, 0] * np.exp(
                theta_i[i - 1] + phi_j[j - 1] + lambda_ij[i - 1, j - 1]
            )

# 3. 結果を確認
print("元の p_ij:")
print(p_ij)
print("\n再構成した p_ij:")
print(p_ij_reconstructed)

# 一致を確認
if np.allclose(p_ij, p_ij_reconstructed, atol=1e-10):
    print("\n与えられた p_ij と再構成した p_ij は一致しました！変換は一対一です。")
else:
    print("\n再構成した p_ij が元の p_ij と一致しません。")

元の p_ij:
[[0.2  0.1  0.05]
 [0.15 0.25 0.05]
 [0.05 0.05 0.1 ]]

再構成した p_ij:
[[0.2  0.1  0.05]
 [0.15 0.25 0.05]
 [0.05 0.05 0.1 ]]

与えられた p_ij と再構成した p_ij は一致しました！変換は一対一です。