# 必要なライブラリのインストール
pip install openai python-dotenv

# プロジェクトディレクトリ作成
mkdir llm_causal_inference
cd llm_causal_inference

# .envファイルを作成（APIキーを安全に管理）
# Windowsの場合: type nul > .env
# Mac/Linuxの場合: touch .env

OPENAI_API_KEY=sk-proj-xxxxxxxxxxxxxxxxxxxxxxxxx

import os
from openai import OpenAI
from dotenv import load_dotenv

# .envファイルから環境変数を読み込み
load_dotenv()

# OpenAIクライアントの初期化
client = OpenAI(
    api_key=os.environ.get("OPENAI_API_KEY")
)

def test_api_connection():
    """API接続テスト"""
    try:
        response = client.chat.completions.create(
            model="gpt-4o-mini",
            messages=[
                {"role": "system", "content": "あなたは因果推論の専門家です。"},
                {"role": "user", "content": "因果推論とは何か、簡潔に説明してください。"}
            ],
            temperature=0.7,
            max_tokens=500
        )
        
        print("✅ API接続成功！")
        print("\n【GPTの回答】")
        print(response.choices[0].message.content)
        
        # トークン使用量を確認
        print(f"\n【トークン使用量】")
        print(f"入力: {response.usage.prompt_tokens} tokens")
        print(f"出力: {response.usage.completion_tokens} tokens")
        print(f"合計: {response.usage.total_tokens} tokens")
        
        # コスト試算（gpt-4o-mini）
        input_cost = (response.usage.prompt_tokens / 1_000_000) * 0.15
        output_cost = (response.usage.completion_tokens / 1_000_000) * 0.60
        total_cost = input_cost + output_cost
        print(f"推定コスト: ${total_cost:.6f} (約{total_cost * 150:.4f}円)")
        
    except Exception as e:
        print(f"❌ エラー発生: {e}")
        print("\n【トラブルシューティング】")
        print("1. APIキーが正しく設定されているか確認")
        print("2. .envファイルがプロジェクトルートにあるか確認")
        print("3. 課金設定が完了しているか確認")

if __name__ == "__main__":
    test_api_connection()

import numpy as np
import pandas as pd
from openai import OpenAI
from dotenv import load_dotenv
import os

load_dotenv()
client = OpenAI(api_key=os.environ.get("OPENAI_API_KEY"))

# Phase 4で学んだT-Learnerの結果を想定（実際の分析結果を使う）
cate_results = {
    '高所得者層': {
        'cate': 0.12,
        'std_error': 0.03,
        'p_value': 0.0001,
        'sample_size': 850,
        'control_mean': 2.5,  # 問い合わせ数（広告なし）
        'treatment_mean': 2.8  # 問い合わせ数（広告あり）
    },
    '中所得者層': {
        'cate': 0.05,
        'std_error': 0.02,
        'p_value': 0.012,
        'sample_size': 1200,
        'control_mean': 2.2,
        'treatment_mean': 2.315
    },
    '低所得者層': {
        'cate': -0.02,
        'std_error': 0.025,
        'p_value': 0.423,
        'sample_size': 650,
        'control_mean': 1.8,
        'treatment_mean': 1.764
    }
}

def create_cate_interpretation_prompt(segment, results):
    """CATE解釈用のプロンプトを作成"""
    prompt = f"""
あなたは不動産マーケティングの因果推論専門家です。
以下のCATE（条件付き平均処理効果）分析結果を、経営層にわかりやすく説明してください。

【分析対象】
セグメント: {segment}
サンプルサイズ: {results['sample_size']}人

【分析結果】
- CATE値: {results['cate']:.3f}
- 標準誤差: {results['std_error']:.3f}
- p値: {results['p_value']:.4f}
- コントロール群平均（広告なし）: {results['control_mean']:.2f}件/月
- 処理群平均（広告あり）: {results['treatment_mean']:.2f}件/月

【説明に含めるべきポイント】
1. この結果が実務的に意味すること（効果の大きさを平易に説明）
2. 統計的有意性の評価（p値の解釈）
3. 経営判断への具体的な示唆（アクションアイテム）
4. 注意すべきリスクや制約条件

専門用語を使わず、マーケティング担当者が理解できる言葉で説明してください。
"""
    return prompt

def interpret_cate_with_llm(segment, results, model="gpt-4o-mini"):
    """LLMを使ってCATE結果を解釈"""
    prompt = create_cate_interpretation_prompt(segment, results)
    
    try:
        response = client.chat.completions.create(
            model=model,
            messages=[
                {"role": "system", "content": "あなたは因果推論とマーケティングの専門家で、複雑な統計結果を経営層にわかりやすく説明することに長けています。"},
                {"role": "user", "content": prompt}
            ],
            temperature=0.3,  # 創造性を抑え、事実ベースの説明にする
            max_tokens=800
        )
        
        interpretation = response.choices[0].message.content
        tokens_used = response.usage.total_tokens
        cost = (response.usage.prompt_tokens / 1_000_000) * 0.15 + \
               (response.usage.completion_tokens / 1_000_000) * 0.60
        
        return {
            'interpretation': interpretation,
            'tokens_used': tokens_used,
            'cost_usd': cost,
            'cost_jpy': cost * 150
        }
        
    except Exception as e:
        return {'error': str(e)}

def generate_full_cate_report():
    """全セグメントのCATE分析レポートを生成"""
    print("=" * 80)
    print("📊 CATE分析結果 自動解釈レポート")
    print("=" * 80)
    
    total_cost = 0
    
    for segment, results in cate_results.items():
        print(f"\n{'='*80}")
        print(f"🎯 セグメント: {segment}")
        print(f"{'='*80}\n")
        
        interpretation_result = interpret_cate_with_llm(segment, results)
        
        if 'error' in interpretation_result:
            print(f"❌ エラー: {interpretation_result['error']}")
            continue
        
        print(interpretation_result['interpretation'])
        print(f"\n【コスト情報】")
        print(f"トークン使用: {interpretation_result['tokens_used']}")
        print(f"コスト: ${interpretation_result['cost_usd']:.6f} (約{interpretation_result['cost_jpy']:.4f}円)")
        
        total_cost += interpretation_result['cost_usd']
    
    print(f"\n{'='*80}")
    print(f"📈 総コスト: ${total_cost:.6f} (約{total_cost * 150:.4f}円)")
    print(f"{'='*80}")

if __name__ == "__main__":
    generate_full_cate_report()

import shap
import matplotlib.pyplot as plt
from sklearn.ensemble import GradientBoostingRegressor
import io
import base64

# Phase 4で学んだSHAP分析の結果を想定
shap_values_summary = {
    'feature_importance': {
        '物件価格': 0.35,
        '駅徒歩分数': 0.22,
        '築年数': 0.18,
        '間取り': 0.12,
        '周辺環境スコア': 0.08,
        '広告出稿額': 0.05
    },
    'top_positive_effects': [
        {'feature': '物件価格', 'value': 5000, 'shap_value': 0.45, 'unit': '万円'},
        {'feature': '駅徒歩分数', 'value': 3, 'shap_value': 0.38, 'unit': '分'},
        {'feature': '広告出稿額', 'value': 50, 'shap_value': 0.15, 'unit': '万円'}
    ],
    'top_negative_effects': [
        {'feature': '築年数', 'value': 25, 'shap_value': -0.32, 'unit': '年'},
        {'feature': '駅徒歩分数', 'value': 15, 'shap_value': -0.28, 'unit': '分'}
    ]
}

def create_shap_interpretation_prompt(shap_summary):
    """SHAP値解釈用のプロンプトを作成"""
    # 特徴量重要度の整形
    importance_text = "\n".join([
        f"- {feature}: {importance:.2%}"
        for feature, importance in shap_summary['feature_importance'].items()
    ])
    
    # 正の効果の整形
    positive_text = "\n".join([
        f"- {effect['feature']}={effect['value']}{effect['unit']}: SHAP値 +{effect['shap_value']:.2f}"
        for effect in shap_summary['top_positive_effects']
    ])
    
    # 負の効果の整形
    negative_text = "\n".join([
        f"- {effect['feature']}={effect['value']}{effect['unit']}: SHAP値 {effect['shap_value']:.2f}"
        for effect in shap_summary['top_negative_effects']
    ])
    
    prompt = f"""
あなたは不動産マーケティングのデータサイエンティストです。
以下のSHAP（SHapley Additive exPlanations）分析結果を、マーケティングチームにわかりやすく説明してください。

【特徴量重要度（問い合わせ数への影響度）】
{importance_text}

【問い合わせ数を増やす要因（正の効果トップ3）】
{positive_text}

【問い合わせ数を減らす要因（負の効果トップ2）】
{negative_text}

【説明に含めるべきポイント】
1. 最も重要な変数は何か、その実務的な意味
2. マーケティング施策で制御可能な変数はどれか
3. 制御不可能な変数（物件特性など）にどう対処すべきか
4. 具体的なアクションプラン（3つ以上）

SHAP値の技術的な説明は不要です。「何をすべきか」に焦点を当ててください。
"""
    return prompt

def interpret_shap_with_llm(shap_summary, model="gpt-4o-mini"):
    """LLMを使ってSHAP値を解釈"""
    prompt = create_shap_interpretation_prompt(shap_summary)
    
    response = client.chat.completions.create(
        model=model,
        messages=[
            {"role": "system", "content": "あなたは機械学習の解釈可能性とマーケティング戦略の専門家です。"},
            {"role": "user", "content": prompt}
        ],
        temperature=0.4,
        max_tokens=1000
    )
    
    return response.choices[0].message.content

# 実行例
print("=" * 80)
print("🔍 SHAP分析結果の自動解釈")
print("=" * 80)
print("\n" + interpret_shap_with_llm(shap_values_summary))

def generate_causal_graph_hypothesis(variables, domain="不動産マーケティング"):
    """LLMを使って因果グラフの仮説を生成"""
    
    variables_text = "、".join(variables)
    
    prompt = f"""
あなたは{domain}の因果推論専門家です。
以下の変数間の因果関係について、DAG（有向非巡回グラフ）の形で仮説を提案してください。

【変数リスト】
{variables_text}

【指示】
1. 各変数ペアについて、因果関係がある場合は「変数A → 変数B」の形で示す
2. 因果関係がないと考えられる場合は明示的に「関係なし」と記載
3. 混乱変数（交絡因子）の可能性がある変数を指摘
4. 各因果関係について、その根拠を簡潔に説明

【注意事項】
- 「相関関係」と「因果関係」を混同しないこと
- 逆因果（Reverse Causality）の可能性を検討すること
- 時間的な前後関係を考慮すること

回答は以下の形式でお願いします：

## 因果関係の仮説
1. 変数A → 変数B（根拠: ...）
2. 変数C → 変数D（根拠: ...）
...

## 混乱変数の可能性
- 変数X: 変数AとBの両方に影響する可能性

## 注意すべき点
- ...
"""
    
    response = client.chat.completions.create(
        model="gpt-4o",  # 複雑な推論にはgpt-4oを使用
        messages=[
            {"role": "system", "content": f"あなたは{domain}と因果推論の専門家で、変数間の因果関係を慎重に分析します。"},
            {"role": "user", "content": prompt}
        ],
        temperature=0.3,
        max_tokens=1500
    )
    
    return response.choices[0].message.content

# 実行例
variables = [
    "物件価格",
    "駅徒歩分数",
    "築年数",
    "広告出稿額",
    "問い合わせ数",
    "成約数",
    "顧客所得層",
    "競合物件数"
]

print("=" * 80)
print("🔗 LLMによる因果グラフ仮説の生成")
print("=" * 80)
hypothesis = generate_causal_graph_hypothesis(variables)
print(hypothesis)

def generate_weekly_report(analysis_results):
    """週次マーケティング分析レポートを自動生成"""
    
    prompt = f"""
あなたは不動産会社のマーケティング分析チームのリーダーです。
今週実施した因果推論分析の結果をまとめ、経営層向けの週次レポートを作成してください。

【今週の分析結果サマリー】

## 1. A/Bテスト結果（Phase 1: RCT分析）
- 処理群（新LP）のコンバージョン率: {analysis_results['rct']['treatment_cvr']:.2%}
- 対照群（旧LP）のコンバージョン率: {analysis_results['rct']['control_cvr']:.2%}
- 平均処理効果（ATE）: {analysis_results['rct']['ate']:.2%}
- p値: {analysis_results['rct']['p_value']:.4f}

## 2. セグメント別効果分析（Phase 4: CATE分析）
高所得者層: +{analysis_results['cate']['high_income']:.2%}（有意）
中所得者層: +{analysis_results['cate']['mid_income']:.2%}（有意）
低所得者層: {analysis_results['cate']['low_income']:.2%}（非有意）

## 3. 時系列効果（Phase 3: DID分析）
- 介入前後の差分: {analysis_results['did']['effect']:.2%}
- 持続期間: {analysis_results['did']['duration']}週間

## 4. 重要変数（Phase 4: SHAP分析）
1位: 物件価格（寄与度35%）
2位: 駅徒歩分数（寄与度22%）
3位: 広告出稿額（寄与度5%）

【レポートの構成】

# エグゼクティブサマリー
（3-4行で今週の最重要発見を記載）

# 主要な発見
## 発見1: [タイトル]
（背景、データ、解釈、示唆を含む）

## 発見2: [タイトル]
（背景、データ、解釈、示唆を含む）

# 推奨アクション
1. 短期施策（今月実施）
2. 中期施策（次四半期）
3. 追加分析の提案

# リスクと制約条件
（分析の限界、注意すべき点）

【スタイル】
- 専門用語は最小限に
- 各発見に「So What?（だから何？）」を明確に
- アクションは具体的かつ測定可能に
- 全体で1000-1200単語
"""
    
    response = client.chat.completions.create(
        model="gpt-4o",
        messages=[
            {"role": "system", "content": "あなたは経営層への報告に長けたマーケティングアナリストです。"},
            {"role": "user", "content": prompt}
        ],
        temperature=0.4,
        max_tokens=2000
    )
    
    return response.choices[0].message.content

# サンプルデータ
sample_analysis = {
    'rct': {
        'treatment_cvr': 0.045,
        'control_cvr': 0.038,
        'ate': 0.007,
        'p_value': 0.032
    },
    'cate': {
        'high_income': 0.012,
        'mid_income': 0.005,
        'low_income': -0.002
    },
    'did': {
        'effect': 0.015,
        'duration': 4
    }
}

print("=" * 80)
print("📊 週次マーケティング分析レポート（自動生成）")
print("=" * 80)
report = generate_weekly_report(sample_analysis)
print(report)

# レポートをMarkdownファイルとして保存
with open("weekly_report.md", "w", encoding="utf-8") as f:
    f.write(f"# 週次マーケティング分析レポート\n")
    f.write(f"生成日時: {pd.Timestamp.now().strftime('%Y-%m-%d %H:%M:%S')}\n\n")
    f.write(report)

print("\n✅ レポートを 'weekly_report.md' に保存しました")

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from lingam import DirectLiNGAM
from scipy import stats

# サンプルデータ生成（不動産マーケティングの例）
np.random.seed(42)
n_samples = 500

# 真の因果構造を定義
# 駅徒歩分数 → 物件価格 → 問い合わせ数 → 成約数
# 築年数 → 物件価格
# 広告出稿額 → 問い合わせ数

# 外生変数（非ガウス分布）
e_walk = np.random.laplace(0, 1, n_samples)  # ラプラス分布
e_age = np.random.exponential(1, n_samples)  # 指数分布
e_price = np.random.laplace(0, 1, n_samples)
e_ad = np.random.exponential(1, n_samples)
e_inquiry = np.random.laplace(0, 1, n_samples)
e_contract = np.random.exponential(1, n_samples)

# 変数生成（真の因果関係に従う）
walk_minutes = 10 + e_walk  # 駅徒歩分数（分）
building_age = 15 + e_age    # 築年数（年）

# 物件価格 = 5000 - 100×駅徒歩 - 50×築年数 + 誤差
property_price = 5000 - 100 * walk_minutes - 50 * building_age + e_price

ad_spending = 30 + e_ad  # 広告出稿額（万円）

# 問い合わせ数 = 10 + 0.01×価格 + 0.5×広告 - 0.3×駅徒歩 + 誤差
inquiry_count = 10 + 0.01 * property_price + 0.5 * ad_spending - 0.3 * walk_minutes + e_inquiry

# 成約数 = 0.3×問い合わせ数 + 誤差
contract_count = 0.3 * inquiry_count + e_contract

# データフレーム作成
df = pd.DataFrame({
    '駅徒歩分数': walk_minutes,
    '築年数': building_age,
    '物件価格': property_price,
    '広告出稿額': ad_spending,
    '問い合わせ数': inquiry_count,
    '成約数': contract_count
})

# 標準化（LiNGAMの前処理として推奨）
df_standardized = (df - df.mean()) / df.std()

print("=" * 80)
print("📊 データの基本統計量")
print("=" * 80)
print(df.describe())

# 非ガウス性の検定（Jarque-Bera検定）
print("\n" + "=" * 80)
print("🔍 非ガウス性の検定（Jarque-Bera検定）")
print("=" * 80)
print("帰無仮説: データは正規分布に従う")
print("p < 0.05 の場合、非ガウス分布と判断（LiNGAMの前提を満たす）\n")

for col in df.columns:
    statistic, p_value = stats.jarque_bera(df[col])
    is_non_gaussian = "✅ 非ガウス" if p_value < 0.05 else "⚠️ ガウス的"
    print(f"{col}: p値 = {p_value:.6f} {is_non_gaussian}")

# DirectLiNGAMの実行
print("\n" + "=" * 80)
print("🔗 DirectLiNGAMによる因果推定")
print("=" * 80)

model = DirectLiNGAM()
model.fit(df_standardized.values)

# 因果順序の推定
causal_order = model.causal_order_
print("\n【推定された因果順序】")
print("（根本原因 → 中間変数 → 最終結果の順）")
for i, idx in enumerate(causal_order):
    print(f"{i+1}. {df.columns[idx]}")

# 隣接行列（因果グラフ）
adjacency_matrix = model.adjacency_matrix_
print("\n【隣接行列（因果グラフ）】")
print("行 → 列への因果効果を表す")
adj_df = pd.DataFrame(
    adjacency_matrix,
    columns=df.columns,
    index=df.columns
)
print(adj_df.round(3))

# 因果関係の可視化
def plot_causal_graph(adjacency_matrix, labels, threshold=0.1):
    """因果グラフを可視化"""
    import networkx as nx
    
    # 有向グラフの作成
    G = nx.DiGraph()
    
    # ノード追加
    for label in labels:
        G.add_node(label)
    
    # エッジ追加（閾値以上の因果効果のみ）
    for i, from_node in enumerate(labels):
        for j, to_node in enumerate(labels):
            weight = adjacency_matrix[i, j]
            if abs(weight) > threshold:
                G.add_edge(from_node, to_node, weight=weight)
    
    # 描画
    plt.figure(figsize=(12, 8))
    pos = nx.spring_layout(G, k=2, iterations=50)
    
    # ノード描画
    nx.draw_networkx_nodes(G, pos, node_color='lightblue', 
                          node_size=3000, alpha=0.9)
    
    # エッジ描画
    edges = G.edges()
    weights = [G[u][v]['weight'] for u, v in edges]
    nx.draw_networkx_edges(G, pos, edgelist=edges, 
                          width=[abs(w)*3 for w in weights],
                          edge_color=weights, edge_cmap=plt.cm.RdYlGn,
                          arrows=True, arrowsize=20, 
                          connectionstyle="arc3,rad=0.1")
    
    # ラベル描画
    nx.draw_networkx_labels(G, pos, font_size=10, font_weight='bold')
    
    # エッジラベル（因果効果の値）
    edge_labels = {(u, v): f"{G[u][v]['weight']:.2f}" 
                   for u, v in edges}
    nx.draw_networkx_edge_labels(G, pos, edge_labels, font_size=8)
    
    plt.title("LiNGAMによる推定因果グラフ", fontsize=16, fontweight='bold')
    plt.axis('off')
    plt.tight_layout()
    plt.savefig('causal_graph_lingam.png', dpi=300, bbox_inches='tight')
    plt.show()
    
    print("\n✅ 因果グラフを 'causal_graph_lingam.png' に保存しました")

plot_causal_graph(adjacency_matrix, df.columns, threshold=0.1)

# 真の因果構造との比較
print("\n" + "=" * 80)
print("🎯 真の因果構造との比較")
print("=" * 80)
print("【真の因果関係】")
print("駅徒歩分数 → 物件価格")
print("築年数 → 物件価格")
print("物件価格 → 問い合わせ数")
print("駅徒歩分数 → 問い合わせ数")
print("広告出稿額 → 問い合わせ数")
print("問い合わせ数 → 成約数")

print("\n【LiNGAMの推定精度】")
true_edges = [
    ('駅徒歩分数', '物件価格'),
    ('築年数', '物件価格'),
    ('物件価格', '問い合わせ数'),
    ('駅徒歩分数', '問い合わせ数'),
    ('広告出稿額', '問い合わせ数'),
    ('問い合わせ数', '成約数')
]

detected_edges = []
for i, from_var in enumerate(df.columns):
    for j, to_var in enumerate(df.columns):
        if abs(adjacency_matrix[i, j]) > 0.1:
            detected_edges.append((from_var, to_var))

true_positives = len(set(true_edges) & set(detected_edges))
false_positives = len(set(detected_edges) - set(true_edges))
false_negatives = len(set(true_edges) - set(detected_edges))

precision = true_positives / (true_positives + false_positives) if (true_positives + false_positives) > 0 else 0
recall = true_positives / (true_positives + false_negatives) if (true_positives + false_negatives) > 0 else 0
f1 = 2 * precision * recall / (precision + recall) if (precision + recall) > 0 else 0

print(f"適合率（Precision）: {precision:.2%}")
print(f"再現率（Recall）: {recall:.2%}")
print(f"F1スコア: {f1:.2%}")

from lingam import ICALiNGAM

# ICA-LiNGAMの実行
print("=" * 80)
print("🔗 ICA-LiNGAMによる因果推定")
print("=" * 80)

ica_model = ICALiNGAM()
ica_model.fit(df_standardized.values)

ica_causal_order = ica_model.causal_order_
ica_adjacency_matrix = ica_model.adjacency_matrix_

print("\n【ICA-LiNGAMの因果順序】")
for i, idx in enumerate(ica_causal_order):
    print(f"{i+1}. {df.columns[idx]}")

print("\n【ICA-LiNGAMの隣接行列】")
ica_adj_df = pd.DataFrame(
    ica_adjacency_matrix,
    columns=df.columns,
    index=df.columns
)
print(ica_adj_df.round(3))

# DirectLiNGAMとICA-LiNGAMの比較
print("\n" + "=" * 80)
print("⚖️ DirectLiNGAM vs ICA-LiNGAM の比較")
print("=" * 80)

comparison_table = []
for i, from_var in enumerate(df.columns):
    for j, to_var in enumerate(df.columns):
        if abs(adjacency_matrix[i, j]) > 0.1 or abs(ica_adjacency_matrix[i, j]) > 0.1:
            comparison_table.append({
                '因果関係': f'{from_var} → {to_var}',
                'DirectLiNGAM': f'{adjacency_matrix[i, j]:.3f}',
                'ICA-LiNGAM': f'{ica_adjacency_matrix[i, j]:.3f}',
                '差分': f'{abs(adjacency_matrix[i, j] - ica_adjacency_matrix[i, j]):.3f}'
            })

comparison_df = pd.DataFrame(comparison_table)
print(comparison_df.to_string(index=False))

print("\n【推奨】")
print("- サンプルサイズが大きい（n > 500）: DirectLiNGAM")
print("- サンプルサイズが小さい（n < 500）: ICA-LiNGAM")
print("- 両方を実行して結果の一致度を確認するのがベストプラクティス")

from lingam import VARLiNGAM
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# 時系列データの生成（不動産マーケティングの週次データ）
np.random.seed(42)
n_periods = 100  # 100週間
lags = 2  # ラグ次数

# 外生変数（非ガウス分布）
e_ad = np.random.laplace(0, 1, n_periods)
e_inquiry = np.random.laplace(0, 1, n_periods)
e_contract = np.random.exponential(1, n_periods)

# 時系列データ生成
ad_spending = np.zeros(n_periods)
inquiry_count = np.zeros(n_periods)
contract_count = np.zeros(n_periods)

# 初期値
ad_spending[0] = 30 + e_ad[0]
inquiry_count[0] = 10 + e_inquiry[0]
contract_count[0] = 3 + e_contract[0]

ad_spending[1] = 30 + 0.3 * ad_spending[0] + e_ad[1]
inquiry_count[1] = 10 + 0.5 * ad_spending[1] + 0.2 * inquiry_count[0] + e_inquiry[1]
contract_count[1] = 3 + 0.3 * inquiry_count[1] + 0.1 * contract_count[0] + e_contract[1]

# t >= 2の時系列データ生成
for t in range(2, n_periods):
    # 真の因果構造
    # 広告 t = 0.3×広告 t-1 + 誤差
    ad_spending[t] = 0.3 * ad_spending[t-1] + e_ad[t]
    
    # 問い合わせ t = 0.5×広告 t + 0.2×問い合わせ t-1 + 0.1×広告 t-2 + 誤差
    # （同時点因果: 広告 t → 問い合わせ t）
    # （ラグ因果: 問い合わせ t-1 → 問い合わせ t、広告 t-2 → 問い合わせ t）
    inquiry_count[t] = (0.5 * ad_spending[t] + 
                       0.2 * inquiry_count[t-1] + 
                       0.1 * ad_spending[t-2] + 
                       e_inquiry[t])
    
    # 成約 t = 0.3×問い合わせ t + 0.1×成約 t-1 + 誤差
    # （同時点因果: 問い合わせ t → 成約 t）
    contract_count[t] = (0.3 * inquiry_count[t] + 
                        0.1 * contract_count[t-1] + 
                        e_contract[t])

# データフレーム作成
ts_df = pd.DataFrame({
    '広告出稿額': ad_spending,
    '問い合わせ数': inquiry_count,
    '成約数': contract_count
})

print("=" * 80)
print("📊 時系列データの基本統計量（週次データ 100週間）")
print("=" * 80)
print(ts_df.describe())

# VAR-LiNGAMの実行
print("\n" + "=" * 80)
print("🔗 VAR-LiNGAMによる時系列因果推定")
print("=" * 80)

var_model = VARLiNGAM(lags=2)
var_model.fit(ts_df.values)

# 同時点の因果関係（instantaneous effects）
print("\n【同時点の因果関係（B行列）】")
print("同じ時点t内での変数間の因果効果")
instantaneous_effects = var_model.adjacency_matrices_[0]
inst_df = pd.DataFrame(
    instantaneous_effects,
    columns=ts_df.columns,
    index=ts_df.columns
)
print(inst_df.round(3))

# ラグ1期の因果関係
print("\n【ラグ1期の因果関係（A1行列）】")
print("t-1時点からt時点への因果効果")
lag1_effects = var_model.adjacency_matrices_[1]
lag1_df = pd.DataFrame(
    lag1_effects,
    columns=ts_df.columns,
    index=ts_df.columns
)
print(lag1_df.round(3))

# ラグ2期の因果関係
print("\n【ラグ2期の因果関係（A2行列）】")
print("t-2時点からt時点への因果効果")
lag2_effects = var_model.adjacency_matrices_[2]
lag2_df = pd.DataFrame(
    lag2_effects,
    columns=ts_df.columns,
    index=ts
    _df.columns,
    index=ts_df.columns
)
print(lag2_df.round(3))

# Phase 3のグレンジャー因果性との比較
print("\n" + "=" * 80)
print("📊 Phase 3との比較: グレンジャー因果性 vs VAR-LiNGAM")
print("=" * 80)

print("\n【グレンジャー因果性（Phase 3）】")
print("- 過去の値が将来の予測に役立つかを検定")
print("- 因果の「方向」は明確だが、同時点の因果は不明")
print("- 例: 「広告 t-1 が問い合わせ t を予測する」")

print("\n【VAR-LiNGAM（Phase 5）】")
print("- ラグ因果 + 同時点因果の両方を推定")
print("- 複雑な時系列因果ネットワークを発見")
print("- 例: 「広告 t → 問い合わせ t」（同時点）+ 「問い合わせ t-1 → 問い合わせ t」（ラグ）")

print("\n【実務での使い分け】")
print("✅ 予測が目的 → グレンジャー因果性（Phase 3）で十分")
print("✅ 因果メカニズムの理解が目的 → VAR-LiNGAM（Phase 5）が必須")

# 因果効果のヒートマップ可視化
def plot_var_lingam_heatmaps(matrices, labels, titles):
    """VAR-LiNGAMの結果をヒートマップで可視化"""
    fig, axes = plt.subplots(1, len(matrices), figsize=(16, 5))
    
    for i, (matrix, title) in enumerate(zip(matrices, titles)):
        sns.heatmap(matrix, annot=True, fmt='.3f', cmap='RdBu_r',
                   center=0, xticklabels=labels, yticklabels=labels,
                   ax=axes[i], cbar_kws={'label': '因果効果'})
        axes[i].set_title(title, fontsize=14, fontweight='bold')
        axes[i].set_xlabel('原因変数', fontsize=12)
        axes[i].set_ylabel('結果変数', fontsize=12)
    
    plt.tight_layout()
    plt.savefig('var_lingam_heatmaps.png', dpi=300, bbox_inches='tight')
    plt.show()
    
    print("\n✅ ヒートマップを 'var_lingam_heatmaps.png' に保存しました")

plot_var_lingam_heatmaps(
    [instantaneous_effects, lag1_effects, lag2_effects],
    ts_df.columns,
    ['同時点因果 (t → t)', 'ラグ1期因果 (t-1 → t)', 'ラグ2期因果 (t-2 → t)']
)

# LLMで結果を解釈
def interpret_var_lingam_with_llm(inst_df, lag1_df, lag2_df):
    """VAR-LiNGAMの結果をLLMで解釈"""
    prompt = f"""
あなたは時系列因果推論の専門家です。
以下のVAR-LiNGAM分析結果を、マーケティングチームに説明してください。

【同時点の因果関係（今週内の効果）】
{inst_df.to_string()}

【ラグ1期の因果関係（先週→今週の効果）】
{lag1_df.to_string()}

【ラグ2期の因果関係（2週前→今週の効果）】
{lag2_df.to_string()}

【説明に含めるべきポイント】
1. 最も強い因果関係は何か（即時効果 vs 遅延効果）
2. 広告の効果が現れるタイミング（即時 or 1週後 or 2週後）
3. 問い合わせから成約までのタイムラグ
4. 実務的な示唆（広告のタイミング、効果測定期間など）

【注意】
- 数値が0.05以下は「効果なし」と判断してください
- 実務で活用できる具体的なアクションを提案してください
"""
    
    response = client.chat.completions.create(
        model="gpt-4o-mini",
        messages=[
            {"role": "system", "content": "あなたは時系列因果分析とマーケティング戦略の専門家です。"},
            {"role": "user", "content": prompt}
        ],
        temperature=0.3,
        max_tokens=1000
    )
    
    return response.choices[0].message.content

print("\n" + "=" * 80)
print("🤖 VAR-LiNGAM結果のLLM解釈")
print("=" * 80)
print(interpret_var_lingam_with_llm(inst_df, lag1_df, lag2_df))

from causallearn.search.ConstraintBased.PC import pc
from causallearn.search.ScoreBased.GES import ges
from causallearn.utils.GraphUtils import GraphUtils
import numpy as np

# Phase 5で使用したデータを再利用
data = df_standardized.values

print("=" * 80)
print("🔍 因果探索アルゴリズムの比較")
print("=" * 80)

# 1. PC algorithm（制約ベース）
print("\n【1. PC algorithm（制約ベース）】")
print("条件付き独立性検定を用いて因果グラフを推定...")

cg_pc = pc(data, alpha=0.05, indep_test='fisherz')
print("\n推定されたグラフ（PC algorithm）:")
print(cg_pc.G)

# 2. GES algorithm（スコアベース）
print("\n【2. GES algorithm（スコアベース）】")
print("BICスコアを最適化して因果グラフを推定...")

record_ges = ges(data, score_func='local_score_BIC')
print("\n推定されたグラフ（GES algorithm）:")
print(record_ges['G'])

# 3. DirectLiNGAM（関数的因果モデル）
print("\n【3. DirectLiNGAM（関数的因果モデル）】")
print("非ガウス性を利用して因果方向を特定...")
print("\n推定された因果順序:")
for i, idx in enumerate(causal_order):
    print(f"{i+1}. {df.columns[idx]}")

# 比較表の作成
comparison_results = {
    'アルゴリズム': ['PC', 'GES', 'DirectLiNGAM'],
    '計算時間（相対）': ['遅い', '中程度', '高速'],
    '因果方向の特定': ['×（CPDAGまで）', '×（CPDAGまで）', '○（DAGまで）'],
    'ガウス性の要求': ['不要', '不要', '必要（非ガウス）'],
    '線形性の要求': ['不要', '不要', '必要'],
    '実務での推奨度': ['低', '低', '高']
}

comparison_df = pd.DataFrame(comparison_results)
print("\n" + "=" * 80)
print("📊 因果探索アルゴリズム比較表")
print("=" * 80)
print(comparison_df.to_string(index=False))

print("\n【用語解説】")
print("CPDAG (Completed Partially Directed Acyclic Graph): 部分的に方向が定まったグラフ")
print("  → 一部のエッジの方向が不定（A-BなのかA→BなのかB→Aなのか不明）")
print("DAG (Directed Acyclic Graph): 完全に方向が定まったグラフ")
print("  → すべてのエッジの方向が確定")

def validate_causal_graph_with_llm(adjacency_matrix, variables, domain="不動産マーケティング"):
    """LLMで因果グラフの妥当性を検証"""
    
    # 推定された因果関係を抽出
    causal_edges = []
    for i, from_var in enumerate(variables):
        for j, to_var in enumerate(variables):
            if abs(adjacency_matrix[i, j]) > 0.1:
                causal_edges.append({
                    'from': from_var,
                    'to': to_var,
                    'strength': adjacency_matrix[i, j]
                })
    
    edges_text = "\n".join([
        f"- {edge['from']} → {edge['to']} (係数: {edge['strength']:.3f})"
        for edge in causal_edges
    ])
    
    prompt = f"""
あなたは{domain}と因果推論の専門家です。
以下は、データから推定された因果関係のリストです。
これらの因果関係について、ドメイン知識の観点から妥当性を評価してください。

【推定された因果関係】
{edges_text}

【評価項目】
1. 各因果関係について
   - ✅ 妥当（ドメイン知識と一致）
   - ⚠️ 要検討（一部疑問あり）
   - ❌ 疑わしい（逆因果や交絡の可能性）

2. 欠けている可能性のある重要な因果関係

3. 逆因果（Reverse Causality）の可能性がある関係

4. 共通原因（交絡因子）の存在が疑われる関係

5. 総合的な評価と改善提案

客観的かつ批判的に評価してください。
"""
    
    response = client.chat.completions.create(
        model="gpt-4o",  # 複雑な推論にはgpt-4oを使用
        messages=[
            {"role": "system", "content": f"あなたは{domain}の専門家で、因果推論の結果を批判的に評価できます。"},
            {"role": "user", "content": prompt}
        ],
        temperature=0.3,
        max_tokens=1500
    )
    
    return response.choices[0].message.content

print("\n" + "=" * 80)
print("🔍 LLMによる因果グラフ妥当性検証")
print("=" * 80)
validation_result = validate_causal_graph_with_llm(adjacency_matrix, df.columns)
print(validation_result)

# 検証結果をファイルに保存
with open("causal_graph_validation.md", "w", encoding="utf-8") as f:
    f.write("# 因果グラフ妥当性検証レポート\n\n")
    f.write(f"生成日時: {pd.Timestamp.now().strftime('%Y-%m-%d %H:%M:%S')}\n\n")
    f.write("## 推定された因果関係\n\n")
    f.write(edges_text + "\n\n")
    f.write("## LLMによる妥当性評価\n\n")
    f.write(validation_result)

print("\n✅ 検証結果を 'causal_graph_validation.md' に保存しました")

import numpy as np
import pandas as pd
import pymc as pm
from statsmodels.tsa.api import VAR
from econml.metalearners import TLearner
from sklearn.ensemble import RandomForestRegressor
from lingam import DirectLiNGAM, VARLiNGAM
from openai import OpenAI
import os
from dotenv import load_dotenv

load_dotenv()
client = OpenAI(api_key=os.environ.get("OPENAI_API_KEY"))

print("=" * 80)
print("🎯 Phase 1-5 統合プロジェクト: 不動産マーケティング戦略の包括的分析")
print("=" * 80)

# サンプルデータの読み込み（実務では実データを使用）
# ここではPhase 1-4で使ったデータを発展させた形で準備
np.random.seed(42)
n = 1000

# 基本変数
df_integrated = pd.DataFrame({
    '顧客ID': range(n),
    '所得層': np.random.choice(['高所得', '中所得', '低所得'], n, p=[0.3, 0.5, 0.2]),
    '年齢': np.random.normal(45, 15, n).clip(20, 80),
    '物件価格': np.random.normal(5000, 1500, n).clip(2000, 10000),
    '駅徒歩分数': np.random.exponential(8, n).clip(1, 30),
    '築年数': np.random.exponential(15, n).clip(0, 50),
    '広告露出': np.random.binomial(1, 0.5, n),  # 処置変数
    '広告出稿額': np.random.normal(30, 10, n).clip(0, 100),
    '問い合わせ': np.random.poisson(2.5, n),
    '成約': np.random.binomial(1, 0.15, n)
})

# 真の因果効果を持つデータに調整
df_integrated.loc[df_integrated['広告露出'] == 1, '問い合わせ'] += np.random.poisson(0.5, 
    (df_integrated['広告露出'] == 1).sum())

print("\n【データ概要】")
print(f"サンプルサイズ: {len(df_integrated)}")
print(f"処置群: {df_integrated['広告露出'].sum()}人")
print(f"対照群: {len(df_integrated) - df_integrated['広告露出'].sum()}人")
print(f"\n基本統計量:\n{df_integrated.describe()}")

# ============================================================
# Phase 1: 因果推論基礎（RCT分析）
# ============================================================
print("\n" + "=" * 80)
print("📊 Phase 1: RCT分析 - 平均処理効果（ATE）の推定")
print("=" * 80)

treated = df_integrated[df_integrated['広告露出'] == 1]['問い合わせ']
control = df_integrated[df_integrated['広告露出'] == 0]['問い合わせ']

ate = treated.mean() - control.mean()
ate_se = np.sqrt(treated.var()/len(treated) + control.var()/len(control))
ate_ci_lower = ate - 1.96 * ate_se
ate_ci_upper = ate + 1.96 * ate_se

print(f"\n処置群平均: {treated.mean():.3f}")
print(f"対照群平均: {control.mean():.3f}")
print(f"ATE: {ate:.3f} (95% CI: [{ate_ci_lower:.3f}, {ate_ci_upper:.3f}])")

from scipy import stats
t_stat, p_value = stats.ttest_ind(treated, control)
print(f"t統計量: {t_stat:.3f}, p値: {p_value:.4f}")

phase1_result = {
    'ate': ate,
    'ate_ci': (ate_ci_lower, ate_ci_upper),
    'p_value': p_value
}

# ============================================================
# Phase 2: ベイズ統計（階層モデルで所得層別の効果を推定）
# ============================================================
print("\n" + "=" * 80)
print("📊 Phase 2: ベイズ階層モデル - 所得層別の効果推定")
print("=" * 80)

# 所得層をカテゴリカル変数に変換
income_map = {'低所得': 0, '中所得': 1, '高所得': 2}
df_integrated['所得層_idx'] = df_integrated['所得層'].map(income_map)

with pm.Model() as hierarchical_model:
    # 階層的事前分布
    mu_alpha = pm.Normal('mu_alpha', mu=2, sigma=1)
    sigma_alpha = pm.HalfNormal('sigma_alpha', sigma=1)
    
    mu_beta = pm.Normal('mu_beta', mu=0.5, sigma=0.5)
    sigma_beta = pm.HalfNormal('sigma_beta', sigma=0.5)
    
    # 所得層別パラメータ
    alpha = pm.Normal('alpha', mu=mu_alpha, sigma=sigma_alpha, shape=3)
    beta = pm.Normal('beta', mu=mu_beta, sigma=sigma_beta, shape=3)
    
    # 線形予測子
    mu = alpha[df_integrated['所得層_idx']] + beta[df_integrated['所得層_idx']] * df_integrated['広告露出']
    
    # 尤度
    sigma = pm.HalfNormal('sigma', sigma=1)
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=df_integrated['問い合わせ'])
    
    # サンプリング（実務では draws=2000, tune=1000 を推奨）
    trace = pm.sample(1000, tune=500, return_inferencedata=True, random_seed=42)

# 事後分布のサマリー
print("\n【所得層別の広告効果（beta）】")
beta_summary = pm.summary(trace, var_names=['beta'])
print(beta_summary)

phase2_result = {
    'low_income_effect': beta_summary.loc['beta[0]', 'mean'],
    'mid_income_effect': beta_summary.loc['beta[1]', 'mean'],
    'high_income_effect': beta_summary.loc['beta[2]', 'mean']
}

# ============================================================
# Phase 3: 時系列分析（VARモデルで相互依存関係を分析）
# ============================================================
print("\n" + "=" * 80)
print("📊 Phase 3: VAR分析 - 広告・問い合わせ・成約の相互依存関係")
print("=" * 80)

# 時系列データの準備（週次集計を想定）
# 実務ではタイムスタンプ付きデータを使用
ts_data = pd.DataFrame({
    '広告出稿額': df_integrated.groupby(df_integrated.index // 10)['広告出稿額'].mean(),
    '問い合わせ数': df_integrated.groupby(df_integrated.index // 10)['問い合わせ'].sum(),
    '成約数': df_integrated.groupby(df_integrated.index // 10)['成約'].sum()
})

# VAR推定
var_model = VAR(ts_data)
var_result = var_model.fit(maxlags=2, ic='aic')

print("\n【VAR(2)モデルの推定結果】")
print(var_result.summary())

# グレンジャー因果性検定
print("\n【グレンジャー因果性検定】")
from statsmodels.tsa.stattools import grangercausalitytests

for target in ['問い合わせ数', '成約数']:
    print(f"\n広告出稿額 → {target} のグレンジャー因果性:")
    test_data = ts_data[['広告出稿額', target]]
    result = grangercausalitytests(test_data, maxlag=2, verbose=False)
    p_value_lag1 = result[1][0]['ssr_ftest'][1]
    print(f"  ラグ1: p値 = {p_value_lag1:.4f}")

phase3_result = {
    'var_summary': var_result.summary(),
    'granger_ad_to_inquiry': p_value_lag1
}

# ============================================================
# Phase 4: 機械学習×因果推論（T-LearnerでCATE推定）
# ============================================================
print("\n" + "=" * 80)
print("📊 Phase 4: T-Learner CATE分析 - セグメント別効果推定")
print("=" * 80)

# 特徴量の準備
X = df_integrated[['年齢', '物件価格', '駅徒歩分数', '築年数', '広告出稿額']]
T = df_integrated['広告露出']
y = df_integrated['問い合わせ']

# T-Learnerの訓練
t_learner = TLearner(
    models=RandomForestRegressor(n_estimators=100, random_seed=42)
)
t_learner.fit(y, T, X=X)

# CATEの予測
cate_pred = t_learner.effect(X)
df_integrated['CATE'] = cate_pred

# 所得層別のCATE
print("\n【所得層別のCATE（平均）】")
cate_by_income = df_integrated.groupby('所得層')['CATE'].agg(['mean', 'std', 'count'])
print(cate_by_income)

phase4_result = {
    'cate_high': cate_by_income.loc['高所得', 'mean'],
    'cate_mid': cate_by_income.loc['中所得', 'mean'],
    'cate_low': cate_by_income.loc['低所得', 'mean']
}

# ============================================================
# Phase 5: 因果探索（DirectLiNGAMで因果ネットワーク発見）
# ============================================================
print("\n" + "=" * 80)
print("📊 Phase 5: DirectLiNGAM - 因果ネットワークの発見")
print("=" * 80)

# 因果探索用のデータ準備
causal_vars = ['物件価格', '駅徒歩分数', '築年数', '広告出稿額', '問い合わせ', '成約']
causal_data = df_integrated[causal_vars].copy()

# 標準化
causal_data_std = (causal_data - causal_data.mean()) / causal_data.std()

# DirectLiNGAM実行
lingam_model = DirectLiNGAM()
lingam_model.fit(causal_data_std.values)

causal_order = lingam_model.causal_order_
adjacency_matrix = lingam_model.adjacency_matrix_

print("\n【推定された因果順序】")
for i, idx in enumerate(causal_order):
    print(f"{i+1}. {causal_vars[idx]}")

print("\n【主要な因果関係（係数 > 0.15）】")
for i, from_var in enumerate(causal_vars):
    for j, to_var in enumerate(causal_vars):
        coef = adjacency_matrix[i, j]
        if abs(coef) > 0.15:
            print(f"{from_var} → {to_var}: {coef:.3f}")

phase5_result = {
    'causal_order': [causal_vars[idx] for idx in causal_order],
    'adjacency_matrix': adjacency_matrix
}

# ============================================================
# LLMで統合レポート生成
# ============================================================
print("\n" + "=" * 80)
print("🤖 LLMによる統合分析レポートの自動生成")
print("=" * 80)

def generate_integrated_report(phase1, phase2, phase3, phase4, phase5):
    """Phase 1-5の統合レポートをLLMで生成"""
    
    prompt = f"""
あなたは不動産会社のマーケティング分析チームのリーダーです。
過去18週間にわたり、5つのPhaseで因果推論分析を実施しました。
その結果を統合し、経営層向けの最終レポートを作成してください。

【Phase 1: RCT分析結果】
- 平均処理効果（ATE）: {phase1['ate']:.3f}
- 95%信頼区間: [{phase1['ate_ci'][0]:.3f}, {phase1['ate_ci'][1]:.3f}]
- p値: {phase1['p_value']:.4f}
→ 広告露出により、平均して問い合わせが{phase1['ate']:.2f}件増加

【Phase 2: ベイズ階層モデル結果】
- 低所得層の効果: {phase2['low_income_effect']:.3f}
- 中所得層の効果: {phase2['mid_income_effect']:.3f}
- 高所得層の効果: {phase2['high_income_effect']:.3f}
→ 所得層によって広告効果に差があることが判明

【Phase 3: VAR時系列分析結果】
- 広告 → 問い合わせのグレンジャー因果性: p = {phase3['granger_ad_to_inquiry']:.4f}
→ 広告の効果は1週間のラグで現れる

【Phase 4: CATE分析結果】
- 高所得層のCATE: {phase4['cate_high']:.3f}
- 中所得層のCATE: {phase4['cate_mid']:.3f}
- 低所得層のCATE: {phase4['cate_low']:.3f}
→ 高所得層への広告が最も効果的

【Phase 5: 因果探索結果】
- 推定された因果順序: {' → '.join(phase5['causal_order'][:4])} ...
- 発見された主要な因果関係:
  * 駅徒歩分数 → 物件価格（負の効果）
  * 広告出稿額 → 問い合わせ数（正の効果）
  * 問い合わせ数 → 成約数（正の効果）

【レポートの構成】

# エグゼクティブサマリー
（3-4行で、18週間の分析から得られた最重要発見を記載）

# 主要な発見
## 発見1: 広告効果の全体像
（Phase 1-3の結果を統合して説明）

## 発見2: セグメント別の戦略
（Phase 2, 4の結果から、どの層に注力すべきかを提言）

## 発見3: 因果メカニズムの解明
（Phase 5の因果探索結果を、実務的に解釈）

# 推奨アクション
1. 短期施策（今月実施）
   - 高所得層への広告予算を○○%増額
   - ○○な物件の訴求を強化
   
2. 中期施策（次四半期）
   - ○○のA/Bテストを実施
   - ○○指標のモニタリング体制構築
   
3. 長期施策（来年度）
   - ○○データ収集の体制整備
   - ○○の予測モデル構築

# 投資対効果（ROI）の試算
（現状の広告予算をどう再配分すれば、問い合わせ数が何%増えるか）

# リスクと制約条件
（分析の限界、注意すべき点）

【スタイル】
- 専門用語は最小限に（Phase番号は使わない）
- 具体的な数値とアクションを明示
- 全体で1500-2000単語
"""
    
    response = client.chat.completions.create(
        model="gpt-4o",
        messages=[
            {"role": "system", "content": "あなたは18週間の学習成果を統合し、実務的な提言ができるマーケティング分析の専門家です。"},
            {"role": "user", "content": prompt}
        ],
        temperature=0.4,
        max_tokens=3000
    )
    
    return response.choices[0].message.content

# 統合レポート生成
integrated_report = generate_integrated_report(
    phase1_result, phase2_result, phase3_result, phase4_result, phase5_result
)

print("\n" + integrated_report)

# レポートをMarkdownファイルとして保存
with open("integrated_final_report.md", "w", encoding="utf-8") as f:
    f.write("# 不動産マーケティング戦略 包括的分析レポート\n")
    f.write(f"## Phase 1-5 統合分析結果\n\n")
    f.write(f"生成日時: {pd.Timestamp.now().strftime('%Y-%m-%d %H:%M:%S')}\n\n")
    f.write(integrated_report)

print("\n" + "=" * 80)
print("✅ 統合レポートを 'integrated_final_report.md' に保存しました")
print("=" * 80)

def generate_final_learning_report():
    """18週間の学習成果を最終レポートとして生成"""
    
    prompt = """
あなたは教育カリキュラムの評価専門家です。
以下の18週間の学習プログラムを完了した学習者の最終レポートを作成してください。

【学習プログラム概要】
期間: 18週間（Phase 0-5）
目標: 因果推論×ベイズ統計×時系列×機械学習×LLMを統合したマーケティングサイエンティストになる

【Phase別学習内容】
- Phase 0 (2週): Python環境構築、統計基礎復習
- Phase 1 (4週): RCT、DID、傾向スコア、回帰不連続
- Phase 2 (4週): ベイズ推論、PyMC、階層モデル
- Phase 3 (3週): ARIMA、VAR、グレンジャー因果性
- Phase 4 (4週): Meta-Learners、DML、SHAP、CATE
- Phase 5 (3週): OpenAI API、LiNGAM、自動レポート生成

【完了した主要プロジェクト】
1. RCTでのA/Bテスト効果測定
2. ベイズ階層モデルでの所得層別分析
3. VARモデルでの時系列因果分析
4. T-LearnerでのCATE推定
5. DirectLiNGAMでの因果ネットワーク発見
6. Phase 1-5統合プロジェクト

【レポートの構成】

# 学習成果サマリー

## 1. 獲得した主要スキル
（Phase別に5-7個のスキルを列挙）

## 2. 実装したプロジェクト
（各Phaseの代表的なプロジェクトを簡潔に説明）

## 3. 克服した困難
（つまずいた点とそれをどう乗り越えたか）

## 4. 最も価値があった学び
（18週間で最も重要だった3つの学び）

# 実務適用計画

## 短期（1-3ヶ月）
- 現在の業務で即座に適用できる手法
- 小規模なA/Bテストやデータ分析の実施

## 中期（3-6ヶ月）
- 部門横断的なデータ分析プロジェクト
- 自動レポートシステムの構築

## 長期（6-12ヶ月）
- 企業全体の因果推論基盤の構築
- 社内勉強会の開催

# 継続学習の方向性

## 1. さらに深めるべき分野
（自分の弱点や興味のある分野）

## 2. 次に学ぶべき技術
（因果推論の最新手法、新しいMLモデル等）

## 3. 推奨リソース
（書籍、オンラインコース、コミュニティ）

# 振り返りと感謝

（18週間の学習を通じて感じたこと、成長の実感）

【スタイル】
- 誠実で前向きなトーン
- 具体的なスキルと成果物を明示
- 全体で1500-2000単語
"""
    
    response = client.chat.completions.create(
        model="gpt-4o",
        messages=[
            {"role": "system", "content": "あなたは学習者の成長を適切に評価し、励ますことができる教育専門家です。"},
            {"role": "user", "content": prompt}
        ],
        temperature=0.5,
        max_tokens=3000
    )
    
    return response.choices[0].message.content

print("=" * 80)
print("📄 18週間学習成果 最終レポート")
print("=" * 80)

final_report = generate_final_learning_report()
print(final_report)

# レポートを保存
with open("final_learning_report.md", "w", encoding="utf-8") as f:
    f.write("# 18週間 因果推論マスタープログラム 最終レポート\n\n")
    f.write(f"完了日: {pd.Timestamp.now().strftime('%Y-%m-%d')}\n\n")
    f.write(final_report)

print("\n✅ 最終レポートを 'final_learning_report.md' に保存しました")