Parameters:
-----------
x : array-like
広告費（入力値）
half_saturation : float
半飽和点（効果が50%に達する広告費）
slope : float
曲率（>1で急速な飽和、<1で緩やかな飽和） """
return x**slope / (half_saturation**slope + x**slope)

# 実データでの使用例
import numpy as np
import matplotlib.pyplot as plt

# シミュレーションデータ
ad_spend = np.linspace(0, 10000, 100)

# 異なるパラメータでの飽和曲線
saturated_1 = hill_saturation(ad_spend, half_saturation=3000, slope=2)
saturated_2 = hill_saturation(ad_spend, half_saturation=5000, slope=1)
saturated_3 = hill_saturation(ad_spend, half_saturation=3000, slope=0.5)

# 可視化
plt.figure(figsize=(10, 6))
plt.plot(ad_spend, saturated_1, label='k=3000, s=2 (急速な飽和)')
plt.plot(ad_spend, saturated_2, label='k=5000, s=1 (標準)')
plt.plot(ad_spend, saturated_3, label='k=3000, s=0.5 (緩やかな飽和)')
plt.xlabel('広告費')
plt.ylabel('飽和後の効果')
plt.title('Hill飽和関数の比較')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

Parameters:
-----------
x : array-like
広告費の時系列データ
decay_rate : float
減衰率（0 < decay_rate < 1） max_lag : int 最大遅延期間 """
adstocked = np.zeros_like(x, dtype=float)
adstocked[0] = x[0]

for t in range(1, len(x)):
adstocked[t] = x[t] + decay_rate * adstocked[t-1]

return adstocked

def weibull_adstock(x, shape, scale, max_lag=8):
"""
Weibull Adstock変換

Parameters:
-----------
x : array-like
広告費の時系列データ
shape : float
形状パラメータ（k）
scale : float
スケールパラメータ（λ）
max_lag : int
最大遅延期間
"""
# Weibull重みの計算
lags = np.arange(0, max_lag + 1)
weights = weibull_min.pdf(lags, c=shape, scale=scale)
weights = weights / weights.sum() # 正規化

# Convolution（畳み込み）でAdstock効果を適用
adstocked = np.convolve(x, weights, mode='same')

return adstocked

# 実例：週次広告費データ
weeks = 20
ad_spend = np.zeros(weeks)
ad_spend[5] = 100 # 第5週に100万円の広告を出稿

# 異なるAdstockモデルの比較
geometric_05 = geometric_adstock(ad_spend, decay_rate=0.5)
geometric_08 = geometric_adstock(ad_spend, decay_rate=0.8)
weibull_effect = weibull_adstock(ad_spend, shape=2, scale=3)

# 可視化
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))

# Geometricの比較
ax1.bar(range(weeks), ad_spend, alpha=0.3, label='元の広告費')
ax1.plot(geometric_05, marker='o', label='Geometric (α=0.5)')
ax1.plot(geometric_08, marker='s', label='Geometric (α=0.8)')
ax1.set_xlabel('週')
ax1.set_ylabel('Adstock効果')
ax1.set_title('Geometric Adstockの減衰率比較')
ax1.legend()
ax1.grid(True, alpha=0.3)

# WeibullとGeometricの比較
ax2.bar(range(weeks), ad_spend, alpha=0.3, label='元の広告費')
ax2.plot(geometric_08, marker='o', label='Geometric (α=0.8)')
ax2.plot(weibull_effect, marker='^', label='Weibull (k=2, λ=3)')
ax2.set_xlabel('週')
ax2.set_ylabel('Adstock効果')
ax2.set_title('GeometricとWeibull Adstockの比較')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

data = pd.DataFrame({
'date': pd.date_range('2023-01-01', periods=n_weeks, freq='W'),
'sales': np.random.poisson(50, n_weeks) + np.random.randn(n_weeks) * 5, # 週次問い合わせ数
'google_ads': np.random.uniform(100, 500, n_weeks), # Google広告費（千円）
'sns_ads': np.random.uniform(50, 300, n_weeks), # SNS広告費（千円）
'content_marketing': np.random.uniform(0, 200, n_weeks), # コンテンツ制作費（千円）
})

# 季節性の追加（オフィス移転は3-4月、9-10月が繁忙期）
data['month'] = data['date'].dt.month
data['is_busy_season'] = ((data['month'].isin([3, 4, 9, 10]))).astype(int)

print("データ準備完了：")
print(data.head())

# Step 2: MMMモデルの構築
mmm = DelayedSaturatedMMM(
date_column="date",
channel_columns=["google_ads", "sns_ads", "content_marketing"],
control_columns=["is_busy_season"], # 季節性をコントロール変数に
adstock="geometric", # AdstockタイプはGeometric
saturation="logistic", # 飽和関数はLogistic
yearly_seasonality=4, # 年次季節性（四半期）
)

# Step 3: 事前分布の設定（Phase 2の知識を活用）
# ドメイン知識に基づいた事前分布の設定
with pm.Model() as model:
# Adstockパラメータの事前分布
# Google広告は即時効果が強い → 低い減衰率
alpha_google = pm.Beta("alpha_google", alpha=2, beta=5) # 平均0.3付近

# SNS広告も即時効果
alpha_sns = pm.Beta("alpha_sns", alpha=2, beta=5)

# コンテンツマーケティングは長期効果 → 高い減衰率
alpha_content = pm.Beta("alpha_content", alpha=5, beta=2) # 平均0.7付近

# 飽和パラメータの事前分布
# 半飽和点（効果が50%に達する広告費）
lam_google = pm.Gamma("lam_google", alpha=2, beta=0.01) # 平均200千円
lam_sns = pm.Gamma("lam_sns", alpha=2, beta=0.015) # 平均133千円
lam_content = pm.Gamma("lam_content", alpha=2, beta=0.02) # 平均100千円

print("事前分布の設定完了")

# Step 4: モデルのフィッティング（MCMCサンプリング）
mmm.fit(
X=data,
y=data["sales"],
chains=4, # MCMCチェーン数
draws=2000, # サンプル数
tune=1000, # バーンイン期間
target_accept=0.9, # 受容率（高いほど精度向上）
random_seed=42,
)

print("MCMCサンプリング完了")

# Step 5: 収束診断
# R-hat統計量（1.01未満が望ましい）
print("\n収束診断（R-hat）：")
print(mmm.fit_result.summary())

# トレースプロット（目視確認）
az.plot_trace(mmm.fit_result)
plt.tight_layout()
plt.show()

# Step 6: モデルの評価
# 事後予測チェック（Posterior Predictive Check）
mmm.plot_posterior_predictive()
plt.show()

# R²スコア
print(f"\nモデルのR²スコア: {mmm.score(X=data, y=data['sales']):.3f}")

print("チャネル別貢献度：")
print(contributions.head())

# Step 2: ROIの計算
def calculate_roi_metrics(mmm_model, data, value_per_conversion=1000000):
"""
チャネル別のROI指標を計算

Parameters:
-----------
mmm_model : DelayedSaturatedMMM
フィッティング済みのMMMモデル
data : pd.DataFrame
元データ
value_per_conversion : float
1コンバージョンあたりの価値（円）
例：不動産仲介なら成約1件=100万円

Returns:
--------
pd.DataFrame
チャネル別のROI指標
"""
results = []

for channel in mmm_model.channel_columns:
# 総広告費
total_spend = data[channel].sum() * 1000 # 千円→円

# 総貢献数（コンバージョン数）
total_contribution = contributions[channel].sum()

# 総貢献額
total_value = total_contribution * value_per_conversion

# 平均ROI
avg_roi = (total_value / total_spend) - 1 if total_spend > 0 else 0

# CPAおよびCPC（参考値）
cpa = total_spend / total_contribution if total_contribution > 0 else np.inf

results.append({
'チャネル': channel,
'総広告費（円）': total_spend,
'総貢献数': total_contribution,
'総貢献額（円）': total_value,
'平均ROI（%）': avg_roi * 100,
'CPA（円）': cpa,
})

return pd.DataFrame(results)

# ROI計算の実行
roi_metrics = calculate_roi_metrics(
mmm,
data,
value_per_conversion=1000000 # オフィス仲介1件=100万円と仮定
)

print("\n=== チャネル別ROI指標 ===")
print(roi_metrics.to_string(index=False))

# Step 3: 可視化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 1. チャネル別貢献度の時系列
contributions.plot(ax=axes[0, 0], title='チャネル別貢献度（時系列）')
axes[0, 0].set_ylabel('貢献度（コンバージョン数）')
axes[0, 0].legend(title='チャネル')
axes[0, 0].grid(True, alpha=0.3)

# 2. チャネル別ROI比較（棒グラフ）
sns.barplot(
data=roi_metrics,
x='チャネル',
y='平均ROI（%）',
ax=axes[0, 1],
palette='viridis'
)
axes[0, 1].set_title('チャネル別 平均ROI')
axes[0, 1].set_ylabel('ROI（%）')
axes[0, 1].axhline(y=0, color='r', linestyle='--', alpha=0.5)
axes[0, 1].grid(True, alpha=0.3)

# 3. 広告費 vs 貢献額の散布図
for channel in mmm.channel_columns:
axes[1, 0].scatter(
data[channel],
contributions[channel],
label=channel,
alpha=0.6
)
axes[1, 0].set_xlabel('広告費（千円）')
axes[1, 0].set_ylabel('貢献度（コンバージョン数）')
axes[1, 0].set_title('広告費と貢献度の関係')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# 4. 貢献度の構成比（円グラフ）
total_contributions = roi_metrics['総貢献数']
axes[1, 1].pie(
total_contributions,
labels=roi_metrics['チャネル'],
autopct='%1.1f%%',
startangle=90
)
axes[1, 1].set_title('チャネル別貢献度の構成比')

plt.tight_layout()
plt.show()

# Step 4: 限界ROI（mROI）の計算
def calculate_marginal_roi(mmm_model, data, channel, delta=10):
"""
特定チャネルの限界ROIを計算

Parameters:
-----------
mmm_model : DelayedSaturatedMMM
フィッティング済みのMMMモデル
data : pd.DataFrame
元データ
channel : str
チャネル名
delta : float
広告費の増分（千円）

Returns:
--------
float
限界ROI（次の1単位の投資に対するリターン）
"""
# 現在の予測
current_pred = mmm_model.predict(data)

# 広告費を少し増やした場合の予測
data_modified = data.copy()
data_modified[channel] = data_modified[channel] + delta
modified_pred = mmm_model.predict(data_modified)

# 限界貢献度
marginal_contribution = (modified_pred - current_pred).mean()

# 限界ROI
marginal_roi = (marginal_contribution / delta) - 1

return marginal_roi

# 各チャネルのmROI計算
print("\n=== 限界ROI（mROI） ===")
for channel in mmm.channel_columns:
mroi = calculate_marginal_roi(mmm, data, channel)
print(f"{channel}: {mroi*100:.2f}%")

Parameters:
-----------
mmm_model : DelayedSaturatedMMM
フィッティング済みのMMMモデル
data : pd.DataFrame
元データ
channels : list
分析対象のチャネルリスト

Returns:
--------
pd.DataFrame
チャネル間の相乗効果マトリックス
"""
results = []

# すべてのチャネルペアについて分析
for ch1, ch2 in combinations(channels, 2):
# ベースライン（両方ゼロ）
data_baseline = data.copy()
data_baseline[ch1] = 0
data_baseline[ch2] = 0
pred_baseline = mmm_model.predict(data_baseline).mean()

# チャネル1のみ
data_ch1 = data.copy()
data_ch1[ch2] = 0
pred_ch1 = mmm_model.predict(data_ch1).mean()
effect_ch1 = pred_ch1 - pred_baseline

# チャネル2のみ
data_ch2 = data.copy()
data_ch2[ch1] = 0
pred_ch2 = mmm_model.predict(data_ch2).mean()
effect_ch2 = pred_ch2 - pred_baseline

# 両方あり
pred_both = mmm_model.predict(data).mean()
effect_both = pred_both - pred_baseline

# 相乗効果の計算
expected_additive = effect_ch1 + effect_ch2
synergy = effect_both - expected_additive
synergy_pct = (synergy / expected_additive) * 100 if expected_additive > 0 else 0

results.append({
'チャネル1': ch1,
'チャネル2': ch2,
'効果1': effect_ch1,
'効果2': effect_ch2,
'期待効果（加算）': expected_additive,
'実際の効果': effect_both,
'相乗効果': synergy,
'相乗効果率（%）': synergy_pct,
})

return pd.DataFrame(results)

# 相乗効果分析の実行
synergy_results = analyze_channel_synergy(
mmm,
data,
mmm.channel_columns
)

print("\n=== チャネル間の相乗効果分析 ===")
print(synergy_results.to_string(index=False))

# ヒートマップで可視化
# 相乗効果率のマトリックスを作成
channels = mmm.channel_columns
synergy_matrix = np.zeros((len(channels), len(channels)))

for i, ch1 in enumerate(channels):
for j, ch2 in enumerate(channels):
if i < j: # 該当する相乗効果を取得
row = synergy_results[
((synergy_results['チャネル1'] == ch1) & (synergy_results['チャネル2'] == ch2)) |
((synergy_results['チャネル1'] == ch2) & (synergy_results['チャネル2'] == ch1))
]
if not row.empty:
synergy_matrix[i, j] = row['相乗効果率（%）'].values[0]
synergy_matrix[j, i] = synergy_matrix[i, j] # 対称

plt.figure(figsize=(10, 8))
sns.heatmap(
synergy_matrix,
annot=True,
fmt='.1f',
cmap='RdYlGn',
center=0,
xticklabels=channels,
yticklabels=channels,
cbar_kws={'label': '相乗効果率（%）'}
)
plt.title('チャネル間の相乗効果マトリックス')
plt.tight_layout()
plt.show()

def generate_mmm_report(roi_metrics, synergy_results, business_context):
"""
LLMを使ってMMM分析結果を経営層向けレポートに変換

Parameters:
-----------
roi_metrics : pd.DataFrame
チャネル別ROI指標
synergy_results : pd.DataFrame
チャネル間の相乗効果分析結果
business_context : str
ビジネスコンテキスト（業種、目標等）

Returns:
--------
str
生成されたレポート
"""

# データをテキスト形式に変換
roi_text = roi_metrics.to_string(index=False)
synergy_text = synergy_results.to_string(index=False)

prompt = f"""
あなたは経験豊富なマーケティングサイエンティストです。
以下のMarketing Mix Modeling（MMM）分析結果を、経営層向けの戦略的レポートに変換してください。

【ビジネスコンテキスト】
{business_context}

【チャネル別ROI分析】
{roi_text}

【チャネル間の相乗効果分析】
{synergy_text}

【レポート要件】
1. エグゼクティブサマリー（3行以内）
2. 主要な発見事項（箇条書き3-5点）
3. 具体的な推奨アクション（予算配分の提案を含む）
4. 期待される効果（定量的に）
5. リスクと対策

専門用語は最小限にし、ビジネスインパクトを明確に伝えてください。
"""

response = client.chat.completions.create(
model="gpt-4",
messages=[
{"role": "system", "content": "経営層向けマーケティングレポート作成の専門家"},
{"role": "user", "content"
": prompt}
],
temperature=0.7,
)

return response.choices[0].message.content

# レポート生成の実行
business_context = """
業種：不動産仲介（オフィス物件専門）
目標：月間問い合わせ数を50件から70件に増加（+40%）
予算制約：月額総広告費800万円まで
期間：次の四半期（3ヶ月）
"""

report = generate_mmm_report(roi_metrics, synergy_results, business_context)

print("\n" + "="*80)
print(" マーケティングミックスモデリング 戦略レポート")
print("="*80 + "\n")
print(report)
print("\n" + "="*80)

# レポートをMarkdownファイルとして保存
with open('mmm_strategic_report.md', 'w', encoding='utf-8') as f:
f.write(f"# マーケティングミックスモデリング 戦略レポート\n\n")
f.write(f"**生成日時**: {pd.Timestamp.now().strftime('%Y年%m月%d日 %H:%M')}\n\n")
f.write(report)

print("\nレポートを 'mmm_strategic_report.md' として保存しました")

np.random.seed(42)
n = 1000

# 交絡因子
company_size = np.random.choice(['small', 'medium', 'large'], n, p=[0.5, 0.3, 0.2])
is_busy = np.random.binomial(1, 0.3, n) # 繁忙期フラグ

# 処置（Google広告出稿）は企業規模と繁忙期に影響される
treatment_prob = 0.3 + 0.2 * (company_size == 'large') + 0.1 * is_busy
treatment = np.random.binomial(1, treatment_prob)

# 結果（問い合わせ数）
# 真の因果効果 = +5件
inquiries = (
10 + # ベースライン
5 * treatment + # 処置効果（真の値）
8 * (company_size == 'large') + # 企業規模の影響
4 * (company_size == 'medium') +
6 * is_busy + # 繁忙期の影響
np.random.randn(n) * 2 # ノイズ
)

data = pd.DataFrame({
'treatment': treatment,
'inquiries': inquiries,
'company_size_medium': (company_size == 'medium').astype(int),
'company_size_large': (company_size == 'large').astype(int),
'is_busy': is_busy,
})

print("サンプルデータ：")
print(data.head())
print(f"\n真の因果効果: +5件")

# ==========================================
# Step 1: Model - 因果グラフ（DAG）の定義
# ==========================================

# グラフをGML形式で定義
causal_graph = """
digraph {
company_size_medium;
company_size_large;
is_busy;
treatment;
inquiries;

company_size_medium -> treatment;
company_size_large -> treatment;
is_busy -> treatment;

company_size_medium -> inquiries;
company_size_large -> inquiries;
is_busy -> inquiries;

treatment -> inquiries;
}
"""

# CausalModelの作成
model = CausalModel(
data=data,
treatment='treatment',
outcome='inquiries',
graph=causal_graph,
common_causes=['company_size_medium', 'company_size_large', 'is_busy']
)

print("\n【Step 1: Model】")
print("因果グラフを定義しました")
model.view_model() # グラフを可視化（要graphviz）

# ==========================================
# Step 2: Identify - 因果効果の識別
# ==========================================

identified_estimand = model.identify_effect(proceed_when_unidentifiable=True)

print("\n【Step 2: Identify】")
print(identified_estimand)

# 識別戦略の確認
print("\n識別された推定量（Estimand）：")
print(f"バックドア基準が適用可能: {identified_estimand.estimands['backdoor'] is not None}")

estimate_psm = model.estimate_effect(
identified_estimand,
method_name="backdoor.propensity_score_matching",
method_params={
'propensity_score_model': 'logit', # ロジスティック回帰
}
)

print(f"PSMによる推定結果: {estimate_psm.value:.2f} (真の値: 5.00)")
print(estimate_psm)

# ==========================================
# 手法2: Inverse Propensity Weighting（IPW）
# ==========================================
print("\n【Step 3: Estimate - IPW】")

estimate_ipw = model.estimate_effect(
identified_estimand,
method_name="backdoor.propensity_score_weighting"
)

print(f"IPWによる推定結果: {estimate_ipw.value:.2f} (真の値: 5.00)")

# ==========================================
# 手法3: Linear Regression（回帰分析）
# ==========================================
print("\n【Step 3: Estimate - Linear Regression】")

estimate_reg = model.estimate_effect(
identified_estimand,
method_name="backdoor.linear_regression"
)

print(f"回帰分析による推定結果: {estimate_reg.value:.2f} (真の値: 5.00)")

# ==========================================
# 手法4: Double Machine Learning（DML）
# ==========================================
# EconMLとの統合
from econml.dml import CausalForestDML

print("\n【Step 3: Estimate - DML (EconML統合)】")

estimate_dml = model.estimate_effect(
identified_estimand,
method_name="backdoor.econml.dml.CausalForestDML",
method_params={
'init_params': {
'n_estimators': 100,
'min_samples_leaf': 10,
}
}
)

print(f"DMLによる推定結果: {estimate_dml.value:.2f} (真の値: 5.00)")

# ==========================================
# すべての手法の比較
# ==========================================
results_comparison = pd.DataFrame({
'手法': ['PSM', 'IPW', '回帰分析', 'DML'],
'推定値': [
estimate_psm.value,
estimate_ipw.value,
estimate_reg.value,
estimate_dml.value
],
'真の値との差': [
abs(estimate_psm.value - 5.0),
abs(estimate_ipw.value - 5.0),
abs(estimate_reg.value - 5.0),
abs(estimate_dml.value - 5.0)
]
})

print("\n=== 推定手法の比較 ===")
print(results_comparison.to_string(index=False))

# 可視化
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(10, 6))
ax.barh(results_comparison['手法'], results_comparison['推定値'], color='skyblue')
ax.axvline(x=5.0, color='red', linestyle='--', label='真の因果効果')
ax.set_xlabel('推定された因果効果')
ax.set_title('DoWhyによる複数手法の比較')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

refute_placebo = model.refute_estimate(
identified_estimand,
estimate_psm,
method_name="placebo_treatment_refuter",
placebo_type="permute" # 処置をランダムに並び替え
)

print(refute_placebo)

# ==========================================
# Refutation 2: Random Common Cause
# ==========================================
print("\n【Step 4: Refute - Random Common Cause】")
print("ランダム変数を交絡因子として追加しても、推定値は変わらないはずです")

refute_random = model.refute_estimate(
identified_estimand,
estimate_psm,
method_name="random_common_cause"
)

print(refute_random)

# ==========================================
# Refutation 3: Data Subset Validation
# ==========================================
print("\n【Step 4: Refute - Data Subset】")
print("データの一部を除外しても、推定値は安定しているはずです")

refute_subset = model.refute_estimate(
identified_estimand,
estimate_psm,
method_name="data_subset_refuter",
subset_fraction=0.8 # 80%のデータで再推定
)

print(refute_subset)

# ==========================================
# Refutation 4: Bootstrap (信頼区間)
# ==========================================
print("\n【Step 4: Refute - Bootstrap】")
print("ブートストラップで信頼区間を計算します")

refute_bootstrap = model.refute_estimate(
identified_estimand,
estimate_psm,
method_name="bootstrap_refuter",
num_simulations=100 # ブートストラップ回数
)

print(refute_bootstrap)

# ==========================================
# すべてのRefutation結果をまとめる
# ==========================================
refutation_summary = pd.DataFrame({
'テスト': [
'元の推定',
'Placebo Treatment',
'Random Common Cause',
'Data Subset',
'Bootstrap'
],
'推定値': [
estimate_psm.value,
refute_placebo.new_effect,
refute_random.new_effect,
refute_subset.new_effect,
refute_bootstrap.new_effect
],
'判定': [
'—',
'Pass' if abs(refute_placebo.new_effect) < 1.0 else 'Fail',
'Pass' if abs(refute_random.new_effect - estimate_psm.value) < 1.0 else 'Fail',
'Pass' if abs(refute_subset.new_effect - estimate_psm.value) < 1.0 else 'Fail',
'Pass'
]
})

print("\n=== Refutation（頑健性チェック）サマリー ===")
print(refutation_summary.to_string(index=False))

# 可視化
fig, ax = plt.subplots(figsize=(12, 6))

colors = ['green' if x == 'Pass' else 'gray' if x == '—' else 'red'
for x in refutation_summary['判定']]

ax.barh(refutation_summary['テスト'], refutation_summary['推定値'], color=colors, alpha=0.7)
ax.axvline(x=5.0, color='red', linestyle='--', linewidth=2, label='真の因果効果')
ax.axvline(x=0, color='black', linestyle='-', linewidth=1, alpha=0.3)
ax.set_xlabel('推定された因果効果')
ax.set_title('DoWhyの頑健性チェック結果')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# ベースモデル（すべてのMeta-LearnersでXGBoostを使用）
base_model = XGBRegressor(
n_estimators=100,
max_depth=3,
learning_rate=0.1,
random_state=42
)

# ==========================================
# S-Learner
# ==========================================
print("\n【S-Learner】")
s_learner = BaseSLearner(learner=base_model)

# フィッティング
s_learner.fit(X=X, treatment=treatment, y=y)

# CATEの推定
cate_s = s_learner.predict(X=X)

# ATEの推定
ate_s = s_learner.estimate_ate(X=X, treatment=treatment, y=y)
print(f"S-Learner ATE: {ate_s[0]:.2f} (真の値: 5.00)")

# ==========================================
# T-Learner
# ==========================================
print("\n【T-Learner】")
t_learner = BaseTLearner(learner=base_model)
t_learner.fit(X=X, treatment=treatment, y=y)
cate_t = t_learner.predict(X=X)
ate_t = t_learner.estimate_ate(X=X, treatment=treatment, y=y)
print(f"T-Learner ATE: {ate_t[0]:.2f} (真の値: 5.00)")

# ==========================================
# X-Learner
# ==========================================
print("\n【X-Learner】")
x_learner = BaseXLearner(learner=base_model)
x_learner.fit(X=X, treatment=treatment, y=y)
cate_x = x_learner.predict(X=X)
ate_x = x_learner.estimate_ate(X=X, treatment=treatment, y=y)
print(f"X-Learner ATE: {ate_x[0]:.2f} (真の値: 5.00)")

# ==========================================
# R-Learner
# ==========================================
print("\n【R-Learner】")
r_learner = BaseRLearner(learner=base_model)
r_learner.fit(X=X, treatment=treatment, y=y)
cate_r = r_learner.predict(X=X)
ate_r = r_learner.estimate_ate(X=X, treatment=treatment, y=y)
print(f"R-Learner ATE: {ate_r[0]:.2f} (真の値: 5.00)")

# ==========================================
# DR-Learner（最も頑健）
# ==========================================
print("\n【DR-Learner】")
dr_learner = BaseDRLearner(learner=base_model)
dr_learner.fit(X=X, treatment=treatment, y=y)
cate_dr = dr_learner.predict(X=X)
ate_dr = dr_learner.estimate_ate(X=X, treatment=treatment, y=y)
print(f"DR-Learner ATE: {ate_dr[0]:.2f} (真の値: 5.00)")

# ==========================================
# すべてのMeta-Learnersの比較
# ==========================================
meta_learner_results = pd.DataFrame({
'手法': ['S-Learner', 'T-Learner', 'X-Learner', 'R-Learner', 'DR-Learner'],
'ATE': [ate_s[0], ate_t[0], ate_x[0], ate_r[0], ate_dr[0]],
'誤差': [
abs(ate_s[0] - 5.0),
abs(ate_t[0] - 5.0),
abs(ate_x[0] - 5.0),
abs(ate_r[0] - 5.0),
abs(ate_dr[0] - 5.0)
]
})

print("\n=== Meta-Learnersの比較 ===")
print(meta_learner_results.to_string(index=False))

# CATEの分布を可視化
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

for i, (name, cate) in enumerate([
('S-Learner', cate_s),
('T-Learner', cate_t),
('X-Learner', cate_x),
('R-Learner', cate_r),
('DR-Learner', cate_dr)
]):
ax = axes[i // 3, i % 3]
ax.hist(cate, bins=30, alpha=0.7, edgecolor='black')
ax.axvline(x=5.0, color='red', linestyle='--', label='真の効果')
ax.axvline(x=cate.mean(), color='blue', linestyle='-', label=f'平均={cate.mean():.2f}')
ax.set_title(name)
ax.set_xlabel('CATE')
ax.set_ylabel('頻度')
ax.legend()
ax.grid(True, alpha=0.3)

# 最後のサブプロットにATE比較
ax = axes[1, 2]
ax.barh(meta_learner_results['手法'], meta_learner_results['ATE'], color='skyblue')
ax.axvline(x=5.0, color='red', linestyle='--', label='真の値')
ax.set_xlabel('ATE')
ax.set_title('Meta-LearnersのATE比較')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

X = data[['company_size_medium', 'company_size_large', 'is_busy']].values
treatment = data['treatment'].values
y = data['conversion'].values

# トレーニングとテストに分割
from sklearn.model_selection import train_test_split

X_train, X_test, t_train, t_test, y_train, y_test = train_test_split(
X, treatment, y, test_size=0.3, random_state=42
)

# ==========================================
# Uplift Tree
# ==========================================
print("\n【Uplift Tree】")

uplift_tree = UpliftTreeClassifier(
max_depth=5,
min_samples_leaf=100,
evaluationFunction='KL' # KL divergence
)

uplift_tree.fit(X_train, treatment=t_train, y=y_train)

# Uplift予測
uplift_tree_pred = uplift_tree.predict(X_test)

print(f"Uplift Treeの予測完了。予測数: {len(uplift_tree_pred)}")

# ==========================================
# Uplift Random Forest
# ==========================================
print("\n【Uplift Random Forest】")

uplift_rf = UpliftRandomForestClassifier(
n_estimators=100,
max_depth=5,
min_samples_leaf=100,
evaluationFunction='KL',
random_state=42
)

uplift_rf.fit(X_train, treatment=t_train, y=y_train)

# Uplift予測
uplift_rf_pred = uplift_rf.predict(X_test)

print(f"Uplift Random Forestの予測完了。予測数: {len(uplift_rf_pred)}")

# ==========================================
# Meta-Learner（比較用）
# ==========================================
from causalml.inference.meta import BaseXLearner
from xgboost import XGBClassifier

x_learner = BaseXLearner(learner=XGBClassifier(random_state=42))
x_learner.fit(X=X_train, treatment=t_train, y=y_train)
uplift_xlearner_pred = x_learner.predict(X=X_test).flatten()

# ==========================================
# Uplift Curve（Gain Chart）で評価
# ==========================================
fig, ax = plt.subplots(figsize=(12, 8))

# Uplift Treeの曲線
plot_gain(
y_test,
uplift_tree_pred,
t_test,
plot_type='gain',
ax=ax
)

ax.set_title('Uplift Gain Curve（累積効果曲線）', fontsize=14)
ax.set_xlabel('ターゲット人口の割合（%）', fontsize=12)
ax.set_ylabel('累積Uplift', fontsize=12)
ax.legend(['Uplift Model', 'Random Targeting'], fontsize=10)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# ==========================================
# AUUC（Area Under Uplift Curve）で定量評価
# ==========================================
auuc_tree = auuc_score(y_test, uplift_tree_pred, t_test)
auuc_rf = auuc_score(y_test, uplift_rf_pred, t_test)
auuc_xlearner = auuc_score(y_test, uplift_xlearner_pred, t_test)

print("\n=== AUUC（Area Under Uplift Curve）スコア ===")
print(f"Uplift Tree: {auuc_tree:.4f}")
print(f"Uplift Random Forest: {auuc_rf:.4f}")
print(f"X-Learner: {auuc_xlearner:.4f}")

# ==========================================
# 顧客セグメンテーション
# ==========================================
def segment_customers(uplift_scores, threshold_high=0.1, threshold_low=-0.05):
"""
Upliftスコアに基づいて顧客を4つのセグメントに分類

Parameters:
-----------
uplift_scores : array-like
各顧客のUpliftスコア
threshold_high : float
Persuadablesの閾値
threshold_low : float
Do Not Disturbsの閾値

Returns:
--------
pd.DataFrame
セグメント分類結果
"""
segments = []

for score in uplift_scores:
if score > threshold_high:
segments.append('Persuadables')
elif score < threshold_low: segments.append('Do Not Disturbs')
elif score >= 0:
segments.append('Sure Things')
else:
segments.append('Lost Causes')

return segments

# セグメンテーション実行
test_data = pd.DataFrame(X_test, columns=['company_size_medium', 'company_size_large', 'is_busy'])
test_data['uplift_score'] = uplift_rf_pred
test_data['segment'] = segment_customers(uplift_rf_pred)

# セグメント別の分布
segment_counts = test_data['segment'].value_counts()

print("\n=== 顧客セグメント分布 ===")
print(segment_counts)
print(f"\nPersuadables（ターゲット層）: {segment_counts.get('Persuadables', 0)}人 ({segment_counts.get('Persuadables', 0)/len(test_data)*100:.1f}%)")

# 可視化
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# Upliftスコアの分布
ax1.hist(test_data['uplift_score'], bins=30, edgecolor='black', alpha=0.7)
ax1.axvline(x=0.1, color='green', linestyle='--', label='Persuadables閾値')
ax1.axvline(x=-0.05, color='red', linestyle='--', label='Do Not Disturbs閾値')
ax1.set_xlabel('Upliftスコア')
ax1.set_ylabel('頻度')
ax1.set_title('Upliftスコアの分布')
ax1.legend()
ax1.grid(True, alpha=0.3)

# セグメント別の円グラフ
colors = {
'Persuadables': '#4caf50',
'Sure Things': '#2196f3',
'Lost Causes': '#9e9e9e',
'Do Not Disturbs': '#f44336'
}
segment_colors = [colors[seg] for seg in segment_counts.index]

ax2.pie(
segment_counts.values,
labels=segment_counts.index,
autopct='%1.1f%%',
colors=segment_colors,
startangle=90
)
ax2.set_title('顧客セグメントの構成比')

plt.tight_layout()
plt.show()

# データ準備（前述のデータを使用）
causal_graph = """
digraph {
company_size_medium -> treatment;
company_size_large -> treatment;
is_busy -> treatment;

company_size_medium -> conversion;
company_size_large -> conversion;
is_busy -> conversion;

treatment -> conversion;
}
"""

model = CausalModel(
data=data,
treatment='treatment',
outcome='conversion',
graph=causal_graph,
common_causes=['company_size_medium', 'company_size_large', 'is_busy']
)

# Identify & Estimate
identified_estimand = model.identify_effect()
estimate = model.estimate_effect(
identified_estimand,
method_name="backdoor.propensity_score_matching"
)

print(f"\nDoWhyによるATE推定: {estimate.value:.4f}")

# 頑健性チェック（簡易版）
refute_result = model.refute_estimate(
identified_estimand,
estimate,
method_name="random_common_cause"
)

print(f"頑健性チェック: {refute_result.new_effect:.4f}")

if abs(refute_result.new_effect - estimate.value) < 0.05:
print("✅ 頑健性チェックPass → CATEの推定に進む")
proceed_to_cate = True
else:
print("❌ 頑健性チェックFail → 因果グラフの見直しが必要")
proceed_to_cate = False

# ==========================================
# ステップ2: CausalMLでCATE推定とUplift Modeling
# ==========================================
if proceed_to_cate:
print("\n【ステップ2: CausalMLでCATE推定】")

# データ準備
X = data[['company_size_medium', 'company_size_large', 'is_busy']].values
treatment = data['treatment'].values
y = data['conversion'].values

# X-Learnerでカスタマイズされた効果を推定
x_learner = BaseXLearner(learner=XGBClassifier(random_state=42))
x_learner.fit(X=X, treatment=treatment, y=y)
cate = x_learner.predict(X=X).flatten()

# 結果をデータフレームに追加
data['cate'] = cate

print(f"CATE推定完了。平均CATE: {cate.mean():.4f}")

# ==========================================
# ステップ3: セグメンテーションと意思決定
# ==========================================
print("\n【ステップ3: ビジネス意思決定】")

# 高CATEセグメントを特定
data['segment'] = pd.cut(
data['cate'],
bins=[-np.inf, 0, 0.1, np.inf],
labels=['Low', 'Medium', 'High']
)

# セグメント別の分析
segment_analysis = data.groupby('segment').agg({
'cate': ['count', 'mean'],
'conversion': 'mean'
}).round(4)

print("\nセグメント別分析：")
print(segment_analysis)

# 推奨アクション
high_cate_count = (data['segment'</span >] == 'High').sum()
total_count = len(data)

print(f"\n【推奨アクション】")
print(f"高CATE層: {high_cate_count}社 ({high_cate_count/total_count*100:.1f}%)")
print(f"→ この層に広告予算を集中させることで、ROIを最大化できます")

# ==========================================
# ステップ4: LLMで最終レポート生成（Phase 5統合）
# ==========================================
print("\n【ステップ4: LLMでレポート生成】")

from openai import OpenAI

client = OpenAI(api_key="your-api-key-here")

report_prompt = f"""
あなたはマーケティングサイエンティストです。以下の因果推論分析結果を経営層向けレポートに変換してください。

【DoWhyによるATE分析】
- 平均処置効果（ATE）: {estimate.value:.4f}
- 頑健性チェック: Pass

【CausalMLによるCATE分析】
- 平均CATE: {cate.mean():.4f}
- 高CATE層: {high_cate_count}社 ({high_cate_count/total_count*100:.1f}%)

【セグメント別分析】
{segment_analysis.to_string()}

【要件】
1. 3行以内のエグゼクティブサマリー
2. 具体的な予算配分提案
3. 期待されるROI改善効果
"""

response = client.chat.completions.create(
model="gpt-4",
messages=[
{"role": "system", "content": "経営層向けレポート作成の専門家"},
{"role": "user", "content": report_prompt}
],
temperature=0.7
)

print("\n" + "="*80)
print(" 因果推論分析 最終レポート")
print("="*80 + "\n")
print(response.choices[0].message.content)
print("\n" + "="*80)

def channel_response(spend, channel_params):
"""
チャネル別のレスポンス関数（MMMで推定）

Parameters:
-----------
spend : float
広告費（千円）
channel_params : dict
チャネルのパラメータ（adstock、saturation）

Returns:
--------
float
予測されるコンバージョン数
"""
# Adstock（キャリーオーバー効果）
decay = channel_params['decay']
adstocked = spend # 簡易版：1期間のみ

# Saturation（飽和効果） - Hill関数
k = channel_params['half_saturation']
s = channel_params['slope']
saturated = adstocked**s / (k**s + adstocked**s)

# スケール調整
scale = channel_params['scale']

return saturated * scale

# MMMから推定されたチャネルパラメータ（例）
channel_params = {
'google_ads': {
'decay': 0.3,
'half_saturation': 300, # 300千円で50%の効果
'slope': 2.0,
'scale': 15 # 最大15件のコンバージョン
},
'sns_ads': {
'decay': 0.4,
'half_saturation': 150,
'slope': 1.5,
'scale': 10
},
'content_marketing': {
'decay': 0.8, # 長期効果
'half_saturation': 100,
'slope': 1.2,
'scale': 20 # 最も高い効果
}
}

# ==========================================
# ステップ2: 最適化問題の定義
# ==========================================

def objective(x):
"""
目的関数：総コンバージョン数を最大化
（最小化問題に変換するため負の値を返す）

Parameters:
-----------
x : array-like
各チャネルの広告費 [google, sns, content]

Returns:
--------
float
負の総コンバージョン数
"""
channels = ['google_ads', 'sns_ads', 'content_marketing']
total_conversions = sum(
channel_response(x[i], channel_params[channels[i]])
for i in range(len(channels))
)
return -total_conversions # 最大化→最小化

# ==========================================
# ステップ3: 制約条件の設定
# ==========================================

# 総予算制約
total_budget = 800 # 800千円（80万円）

# 予算制約：x1 + x2 + x3 <= 800
budget_constraint = LinearConstraint(
np.ones(3), # [1, 1, 1]
0, # 下限
total_budget # 上限
)

# 各チャネルの予算範囲
# [最小出稿額、最大出稿額]
bounds = Bounds(
[50, 50, 0], # 最小：Google 50k、SNS 50k、Content 0k
[500, 300, 400] # 最大：Google 500k、SNS 300k、Content 400k
)

# ==========================================
# ステップ4: 最適化の実行
# ==========================================

# 初期値（現在の予算配分）
x0 = np.array([300, 200, 100]) # Google 300k、SNS 200k、Content 100k

print("\n【現在の予算配分】")
print(f"Google広告: {x0[0]:.0f}千円")
print(f"SNS広告: {x0[1]:.0f}千円")
print(f"コンテンツ: {x0[2]:.0f}千円")
print(f"総予算: {x0.sum():.0f}千円")

current_conversions = -objective(x0)
print(f"予測コンバージョン数: {current_conversions:.1f}件\n")

# 最適化実行
result = minimize(
objective,
x0,
method='SLSQP', # Sequential Least Squares Programming
bounds=bounds,
constraints=[budget_constraint],
options={'ftol': 1e-9, 'maxiter': 1000}
)

print("\n【最適化後の予算配分】")
print(f"Google広告: {result.x[0]:.0f}千円 ({(result.x[0]-x0[0])/x0[0]*100:+.1f}%)")
print(f"SNS広告: {result.x[1]:.0f}千円 ({(result.x[1]-x0[1])/x0[1]*100:+.1f}%)")
print(f"コンテンツ: {result.x[2]:.0f}千円 ({(result.x[2]-x0[2])/x0[2]*100:+.1f}%)")
print(f"総予算: {result.x.sum():.0f}千円")

optimal_conversions = -result.fun
print(f"\n予測コンバージョン数: {optimal_conversions:.1f}件")
print(f"改善: +{optimal_conversions - current_conversions:.1f}件 ({(optimal_conversions/current_conversions-1)*100:.1f}%)")

# ==========================================
# ステップ5: 結果の可視化
# ==========================================
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

channels = ['Google広告', 'SNS広告', 'コンテンツ']

# 1. 予算配分の比較（棒グラフ）
x_pos = np.arange(len(channels))
width = 0.35

axes[0, 0].bar(x_pos - width/2, x0, width, label='現在', alpha=0.8)
axes[0, 0].bar(x_pos + width/2, result.x, width, label='最適化後', alpha=0.8)
axes[0, 0].set_xlabel('チャネル')
axes[0, 0].set_ylabel('予算（千円）')
axes[0, 0].set_title('予算配分の比較')
axes[0, 0].set_xticks(x_pos)
axes[0, 0].set_xticklabels(channels)
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# 2. レスポンス曲線
spend_range = np.linspace(0, 500, 100)
for i, (ch_name, ch_key) in enumerate(zip(channels, channel_params.keys())):
responses = [channel_response(s, channel_params[ch_key]) for s in spend_range]
axes[0, 1].plot(spend_range, responses, label=ch_name)
# 現在の予算位置
axes[0, 1].scatter(x0[i], channel_response(x0[i], channel_params[ch_key]), marker='o', s=100)
# 最適化後の予算位置
axes[0, 1].scatter(result.x[i], channel_response(result.x
[i], channel_params[ch_key]), marker='*', s=200)

axes[0, 1].set_xlabel('広告費（千円）')
axes[0, 1].set_ylabel('コンバージョン数')
axes[0, 1].set_title('チャネル別レスポンス曲線')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# 3. 予算配分比率（円グラフ）
axes[1, 0].pie(result.x, labels=channels, autopct='%1.1f%%', startangle=90)
axes[1, 0].set_title('最適化後の予算配分比率')

# 4. ROI比較
current_roi = [(current_conversions * 1000000 / x0.sum() / 1000) - 1] # 1CV=100万円と仮定
optimal_roi = [(optimal_conversions * 1000000 / result.x.sum() / 1000) - 1]

axes[1, 1].bar(['現在', '最適化後'], [current_roi[0]*100, optimal_roi[0]*100], color=['orange', 'green'], alpha=0.8)
axes[1, 1].set_ylabel('ROI（%）')
axes[1, 1].set_title('ROIの改善')
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

def optimize_for_budget(budget, channel_params):
"""特定の予算で最適配分を計算"""
def objective(x):
channels = list(channel_params.keys())
total = sum(channel_response(x[i], channel_params[channels[i]])
for i in range(len(channels)))
return -total

budget_constraint = LinearConstraint(np.ones(3), 0, budget)
bounds = Bounds([50, 50, 0], [500, 300, 400])
x0 = np.array([budget/3, budget/3, budget/3])

result = minimize(objective, x0, method='SLSQP',
bounds=bounds, constraints=[budget_constraint])

return result.x, -result.fun

# 複数の予算水準でシミュレーション
budget_scenarios = np.linspace(400, 1200, 9) # 40万円〜120万円
results = []

for budget in budget_scenarios:
allocation, conversions = optimize_for_budget(budget, channel_params)
roi = (conversions * 1000000 / (budget * 1000)) - 1 # 1CV=100万円
results.append({
'予算': budget,
'コンバージョン': conversions,
'ROI': roi * 100,
'Google': allocation[0],
'SNS': allocation[1],
'Content': allocation[2]
})

scenario_df = pd.DataFrame(results)

print("\n予算シナリオ分析結果：")
print(scenario_df.to_string(index=False))

# 可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# コンバージョン数とROI
ax1 = axes[0]
ax1_twin = ax1.twinx()

ax1.plot(scenario_df['予算'], scenario_df['コンバージョン'],
marker='o', color='blue', label='コンバージョン数')
ax1_twin.plot(scenario_df['予算'], scenario_df['ROI'],
marker='s', color='red', label='ROI（%）')

ax1.set_xlabel('予算（千円）')
ax1.set_ylabel('コンバージョン数', color='blue')
ax1_twin.set_ylabel('ROI（%）', color='red')
ax1.set_title('予算増減シナリオ：効果とROI')
ax1.grid(True, alpha=0.3)
ax1.legend(loc='upper left')
ax1_twin.legend(loc='upper right')

# チャネル別予算配分の変化
axes[1].stackplot(
scenario_df['予算'],
scenario_df['Google'],
scenario_df['SNS'],
scenario_df['Content'],
labels=['Google広告', 'SNS広告', 'コンテンツ'],
alpha=0.8
)
axes[1].set_xlabel('総予算（千円）')
axes[1].set_ylabel('チャネル別配分（千円）')
axes[1].set_title('予算増減に伴う最適配分の変化')
axes[1].legend(loc='upper left')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# ==========================================
# シナリオ2: チャネル削減の影響
# ==========================================
print("\n【シナリオ2: チャネル削減の影響】")

def scenario_without_channel(budget, channel_params, exclude_channel):
"""特定チャネルを除外した場合の最適配分"""
available_channels = {k: v for k, v in channel_params.items() if k != exclude_channel}
n_channels = len(available_channels)

def objective(x):
channels = list(available_channels.keys())
total = sum(channel_response(x[i], available_channels[channels[i]])
for i in range(n_channels))
return -total

budget_constraint = LinearConstraint(np.ones(n_channels), 0, budget)

# チャネルごとのbounds設定
if exclude_channel == 'google_ads':
bounds = Bounds([50, 0], [300, 400]) # SNS, Content
elif exclude_channel == 'sns_ads':
bounds = Bounds([50, 0], [500, 400]) # Google, Content
else: # content_marketing
bounds = Bounds([50, 50], [500, 300]) # Google, SNS

x0 = np.ones(n_channels) * (budget / n_channels)
result = minimize(objective, x0, method='SLSQP',
bounds=bounds, constraints=[budget_constraint])

return result.x, -result.fun

budget = 800
exclusion_scenarios = []

# 全チャネル使用（ベースライン）
baseline_alloc, baseline_conv = optimize_for_budget(budget, channel_params)
exclusion_scenarios.append({
'シナリオ': '全チャネル（ベースライン）',
'コンバージョン': baseline_conv,
'損失': 0,
'損失率': 0
})

# 各チャネルを除外
for channel in channel_params.keys():
alloc, conv = scenario_without_channel(budget, channel_params, channel)
loss = baseline_conv - conv
exclusion_scenarios.append({
'シナリオ': f'{channel} 除外',
'コンバージョン': conv,
'損失': loss,
'損失率': (loss / baseline_conv) * 100
})

exclusion_df = pd.DataFrame(exclusion_scenarios)

print("\nチャネル削減シナリオ：")
print(exclusion_df.to_string(index=False))

# ==========================================
# シナリオ3: LLMによる対話的分析
# ==========================================
print("\n【シナリオ3: LLMによる対話的What-If分析】")

from openai import OpenAI

client = OpenAI(api_key="your-api-key-here")

def llm_whatif_analysis(question, scenario_df, exclusion_df):
"""LLMを使ってWhat-Ifシナリオに回答"""

context = f"""
あなたはマーケティング予算最適化の専門家です。以下のシミュレーション結果に基づいて、ユーザーの質問に答えてください。

【予算増減シナリオ】
{scenario_df.to_string(index=False)}

【チャネル削減シナリオ】
{exclusion_df.to_string(index=False)}

現在の予算: 800千円（80万円）
現在のコンバージョン: {baseline_conv:.1f}件

ユーザーの質問: {question}

具体的な数値を使って、わかりやすく回答してください。
"""

response = client.chat.completions.create(
model="gpt-4",
messages=[
{"role": "system", "content": "マーケティング予算最適化の専門家"},
{"role": "user", "content": context}
],
temperature=0.5
)

return response.choices[0].message.content

# サンプル質問
questions = [
"予算を100万円に増やすべきですか？それとも現状維持が良いですか？",
"SNS広告をやめてその分をコンテンツに回したらどうなりますか？",
"最もコストパフォーマンスが良い予算水準はいくらですか？"
]

for i, q in enumerate(questions, 1):
print(f"\n質問{i}: {q}")
answer = llm_whatif_analysis(q, scenario_df, exclusion_df)
print(f"\n回答:\n{answer}")
print("\n" + "="*80)