# 2022 Q5(2)  2024.8.14

import sympy as sp

# サンプル数を具体的な数値に設定
n = 3  # 例として n = 3 とする

# シンボリック変数の定義
mu = sp.symbols(f'mu_1:{n+1}')  # mu_1, mu_2, mu_3 を定義
theta, nu = sp.symbols('theta nu')
a = sp.symbols(f'a_1:{n+1}')  # a_1, a_2, a_3 を定義
b1, b2 = sp.symbols('b1 b2')

# 制約条件
constraint_a = sp.Eq(sum(a), 0)
constraint_b = sp.Eq(b1 + b2, 0)

# Z_{i1} = X_i の関係式
eq1 = sp.Eq(nu + a[0] + b1, mu[0])

# Z_{i2} = Y_i の関係式
eq2 = sp.Eq(nu + a[0] + b2, mu[0] + theta)

# b1, b2 の解を求める
solution_b = sp.solve([constraint_b, b2 - b1 - theta], (b1, b2))
b1_solution = solution_b[b1]
b2_solution = solution_b[b2]

# νとa_iを求めるための方程式
nu_a_eq = eq1.subs({b1: b1_solution})

# νを求める
nu_solution = sp.solve(nu_a_eq, nu)[0]

# a_iの解
a_solution = [mu[i] - sp.symbols('mu_bar') for i in range(n)]

# 結果を表示
display(nu_solution, a_solution, b1_solution, b2_solution)