main.py

import numpy as np
from bokeh.layouts import column
from bokeh.plotting import figure, output_file, save
from loguru import logger
from scipy import stats
from scipy.stats import gamma, gaussian_kde, kstest, ttest_ind

from utilities import (
    create_vad_tensors,
    load_data_from_json,
    load_mp3_as_tensors,
    load_vad_files,
)
from vap.turntaking_events import TurnTakingEvents


def stat_test(target1, target2):
    target1_silence_ratio = sum([1 for shift in target1 if shift > 0]) / len(target1)
    target1_overlap_ratio = 1 - target1_silence_ratio
    target2_silence_ratio = sum([1 for shift in target2 if shift > 0]) / len(target2)
    target2_overlap_ratio = 1 - target2_silence_ratio

    logger.info(f"target1 turn shifts num.: {len(target1)}")
    logger.info(f"target1 turn shift max: {max(target1)}")
    logger.info(f"target1 turn shift min: {min(target1)}")
    logger.info(f"target1 turn shift 1st quartile: {np.percentile(target1, 25)}")
    logger.info(f"target1 turn shift 3rd quartile: {np.percentile(target1, 75)}")
    logger.info(f"target1 turn shift mean: {np.mean(target1)}")
    logger.info(f"target1 turn shift median: {np.median(target1)}")
    logger.info(f"target1 turn shift variance: {np.var(target1)}")
    logger.info(f"target1 silence ratio: {target1_silence_ratio}")
    logger.info(f"target1 overlap ratio: {target1_overlap_ratio}")

    logger.info(f"target2 turn shifts num.: {len(target2)}")
    logger.info(f"target2 turn shift max: {max(target2)}")
    logger.info(f"target2 turn shift min: {min(target2)}")
    logger.info(f"target2 turn shift 1st quartile: {np.percentile(target2, 25)}")
    logger.info(f"target2 turn shift 3rd quartile: {np.percentile(target2, 75)}")
    logger.info(f"target2 turn shift mean: {np.mean(target2)}")
    logger.info(f"target2 turn shift median: {np.median(target2)}")
    logger.info(f"target2 turn shift variance: {np.var(target2)}")
    logger.info(f"target2 silence ratio: {target2_silence_ratio}")
    logger.info(f"target2 overlap ratio: {target2_overlap_ratio}")

    target1_normality_p = stats.shapiro(target1).pvalue
    target2_normality_p = stats.shapiro(target2).pvalue

    logger.info(f"target1 normality test p-value: {target1_normality_p}")
    logger.info(f"target2 normality test p-value: {target2_normality_p}")

    # 分布の比較（正規分布であるかによる条件分岐）
    if target1_normality_p > 0.05 and target2_normality_p > 0.05:
        # 両グループが正規分布に従う場合：t検定
        t_stat, p_value, df = ttest_ind(target1, target2, usevar="unequal")
        logger.info(f"T-test p-value: {p_value}")
    else:
        # 正規分布でない場合：Mann-Whitney U検定
        u_stat, p_value = stats.mannwhitneyu(target1, target2)
        logger.info(f"Mann-Whitney U test p-value: {p_value}")

    # 効果量の計算 (Cohen's d)
    cohen_d = (np.mean(target1) - np.mean(target2)) / np.sqrt(
        (np.var(target1) + np.var(target2)) / 2
    )
    logger.info(f"Cohen's d: {cohen_d}")


def plot_turn_shift_distribution(shift_durations, title):
    """
    ターンシフトの分布をbokehで描画する関数
    """
    hist, edges = np.histogram(
        shift_durations, bins=np.arange(-4.0, 4.2, 0.2)
    )  # 0.2秒間隔のbin

    p = figure(
        title=title,
        x_axis_label="Turn Shift Duration (s)",
        y_axis_label="Probability Density",
        width=800,
        height=600,
        x_range=(-4, 4),
    )
    p.quad(
        top=hist,
        bottom=0,
        left=edges[:-1],
        right=edges[1:],
        line_color="white",
        alpha=0.7,
    )

    return p


def analysis_pesonal():
    data = load_data_from_json("data.json")
    personal_dict = {}

    # 統計情報をファイルごとに計算
    for d in data:
        user0_vad, user1_vad = load_vad_files(d["audio_dir"])
        vad_tensor = create_vad_tensors(
            user0_vad, user1_vad
        )  # [(発話開始時刻(秒), 発話終了時刻(秒)), ...]

        TTE = TurnTakingEvents(
            vad_tensor,
            frame_rate=100,
            pre_offset_shift=0.3,
            post_onset_shift=0.3,
            pre_offset_hold=0.3,
            post_onset_hold=0.3,
            min_silence=0.0,
            pre_silence=1.0,
            post_silence=2.0,
            backchannel_threshold=1.0,
        )
        turntaking_events = TTE.get_turntaking_events()
        expert_id = d["expert_user"]["id"]
        novice_id = d["novice_user"]["id"]

        if expert_id not in personal_dict:
            personal_dict[expert_id] = []
        if novice_id not in personal_dict:
            personal_dict[novice_id] = []

        # Expert (user1 -> user0 のターンシフト)
        for start, end in turntaking_events["turn_shift_silence_user1_to_user0"]:
            personal_dict[expert_id].append((end - start) / 100.0)  # 100 Hz -> 秒
        for start, end in turntaking_events["turn_shift_overlap_user1_to_user0"]:
            personal_dict[expert_id].append(-(end - start) / 100.0)  # 100 Hz -> 秒

        # Novice (user0 -> user1 のターンシフト)
        for start, end in turntaking_events["turn_shift_silence_user0_to_user1"]:
            personal_dict[novice_id].append((end - start) / 100.0)  # 100 Hz -> 秒
        for start, end in turntaking_events["turn_shift_overlap_user0_to_user1"]:
            personal_dict[novice_id].append(-(end - start) / 100.0)  # 100 Hz -> 秒

    plots = []
    output_file("output/turn_shift_statistics_personal.html")

    for user_id, turn_shifts in personal_dict.items():
        silence_ratio = sum([1 for shift in turn_shifts if shift > 0]) / len(
            turn_shifts
        )
        overlap_ratio = 1 - silence_ratio

        plots.append(
            plot_turn_shift_distribution(
                turn_shifts, f"Turn Shift Distribution (user_id={user_id})"
            )
        )

        # 統計情報の出力
        logger.info(f"User ID: {user_id}")
        logger.info(f"Turn shifts num.: {len(turn_shifts)}")
        logger.info(f"Turn shift mean: {np.mean(turn_shifts)}")
        logger.info(f"Turn shift median: {np.median(turn_shifts)}")
        logger.info(f"Turn shift variance: {np.var(turn_shifts)}")
        logger.info(f"Silence ratio: {silence_ratio}")
        logger.info(f"Overlap ratio: {overlap_ratio}")

    save(column(*plots))

    """
    高外向性に分類されるユーザID: 1, 2, 10, 12, 13, 16, 18, 19, 20
    中外向性に分類されるユーザID: 6, 14, 15
    低外向性に分類されるユーザID: 3, 4, 5, 7, 8, 9, 11, 17
    """

    high_extroversion = []
    low_extroversion = []
    middle_extroversion = []

    for user_id, turn_shifts in personal_dict.items():
        if user_id in [1, 2, 10, 12, 13, 16, 18, 19, 20]:
            high_extroversion.extend(turn_shifts)
        elif user_id in [3, 4, 5, 7, 8, 9, 11, 17]:
            low_extroversion.extend(turn_shifts)
        else:
            middle_extroversion.extend(turn_shifts)
    logger.info("Extroversion, target1: high, target2: low")
    stat_test(high_extroversion, low_extroversion)

    """
    高開放性に分類されるユーザID: 1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12, 14, 15, 17, 18, 19
    中開放性に分類されるユーザID: 16, 20
    低開放性に分類されるユーザID: 9, 13
    """

    high_openness = []
    low_openness = []
    middle_openness = []

    for user_id, turn_shifts in personal_dict.items():
        if user_id in [1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12, 14, 15, 17, 18, 19]:
            high_openness.extend(turn_shifts)
        elif user_id in [16, 20]:
            middle_openness.extend(turn_shifts)
        else:
            low_openness.extend(turn_shifts)
    logger.info("Openness, target1: high, target2: low")
    stat_test(high_openness, low_openness)

    """
    高協調性に分類されるユーザID: 1, 4, 5, 10, 13, 14, 15, 16, 17, 18, 20
    中協調性に分類されるユーザID: 2, 3, 6, 7, 8, 9
    低協調性に分類されるユーザID: 11, 12, 19
    """

    high_agreeableness = []
    low_agreeableness = []
    middle_agreeableness = []

    for user_id, turn_shifts in personal_dict.items():
        if user_id in [1, 4, 5, 10, 13, 14, 15, 16, 17, 18, 20]:
            high_agreeableness.extend(turn_shifts)
        elif user_id in [2, 3, 6, 7, 8, 9]:
            middle_agreeableness.extend(turn_shifts)
        else:
            low_agreeableness.extend(turn_shifts)
    logger.info("Agreeableness, target1: high, target2: low")
    stat_test(high_agreeableness, low_agreeableness)

    """
    高勤勉性に分類されるユーザID: 4, 10, 17
    中勤勉性に分類されるユーザID: 2, 12, 18, 19, 20
    低勤勉性に分類されるユーザID: 1, 3, 5, 6, 7, 8, 9, 11, 13, 14, 15, 16
    """

    high_conscientiousness = []
    low_conscientiousness = []
    middle_conscientiousness = []

    for user_id, turn_shifts in personal_dict.items():
        if user_id in [4, 10, 17]:
            high_conscientiousness.extend(turn_shifts)
        elif user_id in [2, 12, 18, 19, 20]:
            middle_conscientiousness.extend(turn_shifts)
        else:
            low_conscientiousness.extend(turn_shifts)
    logger.info("Conscientiousness, target1: high, target2: low")
    stat_test(high_conscientiousness, low_conscientiousness)

    """
    高神経症傾向に分類されるユーザID: 3, 4, 5, 6, 7, 8, 10, 11, 14, 17, 18, 20
    低神経症傾向に分類されるユーザID: 1, 2, 9, 12, 13, 15, 16, 19    
    """

    high_neuroticism = []
    low_neuroticism = []

    for user_id, turn_shifts in personal_dict.items():
        if user_id in [3, 4, 5, 6, 7, 8, 10, 11, 14, 17, 18, 20]:
            high_neuroticism.extend(turn_shifts)
        else:
            low_neuroticism.extend(turn_shifts)
    logger.info("Neuroticism, target1: high, target2: low")
    stat_test(high_neuroticism, low_neuroticism)

    plots = []
    output_file("output/turn_shift_statistics_BIG5.html")
    plots.append(
        plot_turn_shift_distribution(
            high_extroversion, "Turn Shift Distribution (high extroversion)"
        )
    )
    plots.append(
        plot_turn_shift_distribution(
            low_extroversion, "Turn Shift Distribution (low extroversion)"
        )
    )
    plots.append(
        plot_turn_shift_distribution(
            high_openness, "Turn Shift Distribution (high openness)"
        )
    )
    plots.append(
        plot_turn_shift_distribution(
            low_openness, "Turn Shift Distribution (low openness)"
        )
    )
    plots.append(
        plot_turn_shift_distribution(
            high_agreeableness, "Turn Shift Distribution (high agreeableness)"
        )
    )
    plots.append(
        plot_turn_shift_distribution(
            low_agreeableness, "Turn Shift Distribution (low agreeableness)"
        )
    )
    plots.append(
        plot_turn_shift_distribution(
            high_conscientiousness, "Turn Shift Distribution (high conscientiousness)"
        )
    )
    plots.append(
        plot_turn_shift_distribution(
            low_conscientiousness, "Turn Shift Distribution (low conscientiousness)"
        )
    )
    plots.append(
        plot_turn_shift_distribution(
            high_neuroticism, "Turn Shift Distribution (high neuroticism)"
        )
    )
    plots.append(
        plot_turn_shift_distribution(
            low_neuroticism, "Turn Shift Distribution (low neuroticism)"
        )
    )
    save(column(*plots))

    return (
        high_extroversion,
        low_extroversion,
        high_openness,
        low_openness,
        high_agreeableness,
        low_agreeableness,
        high_conscientiousness,
        low_conscientiousness,
        high_neuroticism,
        low_neuroticism,
    )


def analysis_relationship():
    data = load_data_from_json("data.json")

    friend_turn_shifts = []
    stranger_turn_shifts = []

    # 統計情報をファイルごとに計算
    for d in data:
        user0_vad, user1_vad = load_vad_files(d["audio_dir"])
        vad_tensor = create_vad_tensors(
            user0_vad, user1_vad
        )  # [(発話開始時刻(秒), 発話終了時刻(秒)), ...]
        TTE = TurnTakingEvents(
            vad_tensor,
            frame_rate=100,
            pre_offset_shift=0.3,
            post_onset_shift=0.3,
            pre_offset_hold=0.3,
            post_onset_hold=0.3,
            min_silence=0.0,
            pre_silence=1.0,
            post_silence=2.0,
            backchannel_threshold=1.0,
        )
        turntaking_events = TTE.get_turntaking_events()
        relationship = d["relationships"]  # 'friend' or 'stranger'

        if relationship == "friend":
            # Expert (user1 -> user0 のターンシフト)
            for start, end in turntaking_events["turn_shift_silence_user1_to_user0"]:
                friend_turn_shifts.append((end - start) / 100.0)  # 100 Hz -> 秒
            for start, end in turntaking_events["turn_shift_overlap_user1_to_user0"]:
                friend_turn_shifts.append(-(end - start) / 100.0)  # 100 Hz -> 秒

            # Novice (user0 -> user1 のターンシフト)
            for start, end in turntaking_events["turn_shift_silence_user0_to_user1"]:
                friend_turn_shifts.append((end - start) / 100.0)  # 100 Hz -> 秒
            for start, end in turntaking_events["turn_shift_overlap_user0_to_user1"]:
                friend_turn_shifts.append(-(end - start) / 100.0)  # 100 Hz -> 秒
        else:
            # Expert (user1 -> user0 のターンシフト)
            for start, end in turntaking_events["turn_shift_silence_user1_to_user0"]:
                stranger_turn_shifts.append((end - start) / 100.0)  # 100 Hz -> 秒
            for start, end in turntaking_events["turn_shift_overlap_user1_to_user0"]:
                stranger_turn_shifts.append(-(end - start) / 100.0)  # 100 Hz -> 秒

            # Novice (user0 -> user1 のターンシフト)
            for start, end in turntaking_events["turn_shift_silence_user0_to_user1"]:
                stranger_turn_shifts.append((end - start) / 100.0)  # 100 Hz -> 秒
            for start, end in turntaking_events["turn_shift_overlap_user0_to_user1"]:
                stranger_turn_shifts.append(-(end - start) / 100.0)  # 100 Hz -> 秒

    # グラフの作成
    output_file("output/turn_shift_statistics_relationship.html")

    plots = []
    plots.append(
        plot_turn_shift_distribution(
            friend_turn_shifts, "Turn Shift Distribution (friend)"
        )
    )
    plots.append(
        plot_turn_shift_distribution(
            stranger_turn_shifts, "Novice Turn Shift Distribution (stranger)"
        )
    )

    save(column(*plots))

    stat_test(friend_turn_shifts, stranger_turn_shifts)

    return friend_turn_shifts, stranger_turn_shifts


def analysis_speaker():
    data = load_data_from_json("data.json")

    expert_turn_shifts = []
    novice_turn_shifts = []
    expert_va_durations = []
    novice_va_durations = []

    # 統計情報をファイルごとに計算
    for d in data:
        user0_vad, user1_vad = load_vad_files(d["audio_dir"])
        vad_tensor = create_vad_tensors(
            user0_vad, user1_vad
        )  # [(発話開始時刻(秒), 発話終了時刻(秒)), ...]
        TTE = TurnTakingEvents(
            vad_tensor,
            frame_rate=100,
            pre_offset_shift=0.3,
            post_onset_shift=0.3,
            pre_offset_hold=0.3,
            post_onset_hold=0.3,
            min_silence=0.0,
            pre_silence=1.0,
            post_silence=2.0,
            backchannel_threshold=1.0,
        )
        turntaking_events = TTE.get_turntaking_events()

        for start, end in user0_vad:
            expert_va_durations.append(end - start)
        for start, end in user1_vad:
            novice_va_durations.append(end - start)

        # Expert (user1 -> user0 のターンシフト)
        for start, end in turntaking_events["turn_shift_silence_user1_to_user0"]:
            expert_turn_shifts.append((end - start) / 100.0)  # 100 Hz -> 秒
        for start, end in turntaking_events["turn_shift_overlap_user1_to_user0"]:
            expert_turn_shifts.append(-(end - start) / 100.0)  # 100 Hz -> 秒

        # Novice (user0 -> user1 のターンシフト)
        for start, end in turntaking_events["turn_shift_silence_user0_to_user1"]:
            novice_turn_shifts.append((end - start) / 100.0)  # 100 Hz -> 秒
        for start, end in turntaking_events["turn_shift_overlap_user0_to_user1"]:
            novice_turn_shifts.append(-(end - start) / 100.0)  # 100 Hz -> 秒

    # グラフの作成
    output_file("output/turn_shift_statistics_speaker.html")

    plots = []
    plots.append(
        plot_turn_shift_distribution(
            expert_turn_shifts, "Expert Turn Shift Distribution"
        )
    )
    plots.append(
        plot_turn_shift_distribution(
            novice_turn_shifts, "Novice Turn Shift Distribution"
        )
    )

    save(column(*plots))

    stat_test(expert_turn_shifts, novice_turn_shifts)

    return expert_turn_shifts, novice_turn_shifts


def analysis_all():
    data = load_data_from_json("data.json")

    audio_duration = 0
    turn_shifts = []

    # 統計情報をファイルごとに計算
    for d in data:
        user0_vad, user1_vad = load_vad_files(d["audio_dir"])
        vad_tensor = create_vad_tensors(
            user0_vad, user1_vad
        )  # [(発話開始時刻(秒), 発話終了時刻(秒)), ...]
        user0_waveform, user1_waveform = load_mp3_as_tensors(d["audio_dir"])
        TTE = TurnTakingEvents(
            vad_tensor,
            frame_rate=100,
            pre_offset_shift=0.3,
            post_onset_shift=0.3,
            pre_offset_hold=0.3,
            post_onset_hold=0.3,
            min_silence=0.0,
            pre_silence=1.0,
            post_silence=2.0,
            backchannel_threshold=1.0,
        )
        turntaking_events = TTE.get_turntaking_events()

        audio_duration += user0_waveform.size(1) / 16000  # 音声の長さ

        for start, end in turntaking_events["turn_shift_silence_user0_to_user1"]:
            turn_shifts.append((end - start) / 100.0)  # 100 Hz -> 秒

        for start, end in turntaking_events["turn_shift_silence_user1_to_user0"]:
            turn_shifts.append((end - start) / 100.0)  # 100 Hz -> 秒

        for start, end in turntaking_events["turn_shift_overlap_user0_to_user1"]:
            turn_shifts.append(-(end - start) / 100.0)  # 100 Hz -> 秒

        for start, end in turntaking_events["turn_shift_overlap_user1_to_user0"]:
            turn_shifts.append(-(end - start) / 100.0)  # 100 Hz -> 秒

    silence_ratio = sum([1 for shift in turn_shifts if shift > 0]) / len(turn_shifts)
    overlap_ratio = 1 - silence_ratio

    # グラフの作成
    output_file("output/turn_shift_statistics.html")

    plots = []
    plots.append(
        plot_turn_shift_distribution(turn_shifts, "Total Turn Shift Distribution")
    )

    save(column(*plots))

    # 統計情報の出力
    logger.info(f"Total audio duration: {audio_duration} seconds")
    logger.info(f"Total turn shifts num.: {len(turn_shifts)}")
    logger.info(f"Turn shift max: {max(turn_shifts)}")
    logger.info(f"Turn shift min: {min(turn_shifts)}")
    logger.info(f"Turn shift 1st quartile: {np.percentile(turn_shifts, 25)}")
    logger.info(f"Turn shift 3rd quartile: {np.percentile(turn_shifts, 75)}")
    logger.info(f"Turn shift mean: {np.mean(turn_shifts)}")
    logger.info(f"Turn shift median: {np.median(turn_shifts)}")
    logger.info(f"Turn shift variance: {np.var(turn_shifts)}")
    logger.info(f"Silence ratio: {silence_ratio}")
    logger.info(f"Overlap ratio: {overlap_ratio}")

    return turn_shifts


def evaluate_gamma_fit(turn_shifts, params, label=""):
    """
    KS 検定でフィット具合を評価し、logger.info する。
    """
    shape, loc, scale = params["shape"], params["loc"], params["scale"]
    D, p_value = kstest(turn_shifts, "gamma", args=(shape, loc, scale))
    logger.info(f"{label} gamma KS test D-statistic: {D:.4f}, p-value: {p_value:.4f}")


def fit_gamma_distribution(turn_shifts):
    """
    ターンシフトデータをガンマ分布に近似させる関数。

    パラメータ:
    turn_shifts (list): ターンシフトの速度データのリスト

    戻り値:
    dict: ガンマ分布のパラメータ
    """
    # データをNumPy配列に変換
    turn_shifts = np.array(turn_shifts)

    # ガンマ分布にフィットさせる
    shape, loc, scale = gamma.fit(turn_shifts)

    # パラメータを辞書形式で返す
    fitted_params = {"shape": shape, "loc": loc, "scale": scale}

    return fitted_params


def plot_gamma_fit_html(turn_shifts, gamma_params):
    """
    ターンシフトデータとガンマ分布にフィットさせた結果を同じグラフにプロットし、HTMLで出力する関数。

    パラメータ:
    turn_shifts (list): ターンシフトの速度データのリスト
    gamma_params (dict): ガンマ分布のパラメータ
    """
    # データをNumPy配列に変換
    turn_shifts = np.array(turn_shifts)

    # 元データのヒストグラムを計算
    hist, edges = np.histogram(turn_shifts, bins=50, density=True)

    # ガンマ分布の近似を計算
    x = np.linspace(min(turn_shifts), max(turn_shifts), 1000)
    shape, loc, scale = (
        gamma_params["shape"],
        gamma_params["loc"],
        gamma_params["scale"],
    )
    pdf_fitted = gamma.pdf(x, shape, loc=loc, scale=scale)

    # Bokehの出力をHTMLに設定
    output_file("output/gamma_fit.html")

    # Bokehでプロット
    p = figure(
        # title="Original Data and Gamma Fit",
        x_axis_label="Turn-shift Speed",
        y_axis_label="Probability Density",
        width=800,
        height=400,
        x_range=(-4, 4),  # 横軸の範囲を指定
        y_range=(0, 1.2),  # 縦軸の範囲を指定
    )

    # ラベルと目盛りのフォントサイズを設定
    p.xaxis.axis_label_text_font_size = "16pt"
    p.yaxis.axis_label_text_font_size = "16pt"
    p.xaxis.major_label_text_font_size = "14pt"
    p.yaxis.major_label_text_font_size = "14pt"
    p.legend.label_text_font_size = "16pt"

    # 元のヒストグラムデータ
    p.quad(
        top=hist,
        bottom=0,
        left=edges[:-1],
        right=edges[1:],
        fill_color="navy",
        line_color="white",
        alpha=0.5,
        legend_label="Original Data",
    )

    # フィットしたガンマ分布
    p.line(x, pdf_fitted, line_color="orange", line_width=2, legend_label="Gamma Fit")

    # 凡例の配置
    p.legend.location = "top_right"

    # グラフを保存
    save(p)


def plot_gamma_fit_comparison_html(
    turn_shifts,
    base_params,
    adjusted_params,
    filename="output/gamma_fit_comparison.html",
):
    """
    元データのヒストグラム＋
    オリジナルフィット＋アジャスト後フィット を同一グラフにプロットして HTML 出力する。
    """
    turn_shifts = np.array(turn_shifts)
    hist, edges = np.histogram(turn_shifts, bins=50, density=True)
    x = np.linspace(min(turn_shifts), max(turn_shifts), 1000)

    # 各PDF を計算
    base_pdf = gamma.pdf(
        x,
        base_params["shape"],
        loc=base_params["loc"],
        scale=base_params["scale"],
    )
    adj_pdf = gamma.pdf(
        x,
        adjusted_params["shape"],
        loc=adjusted_params["loc"],
        scale=adjusted_params["scale"],
    )

    output_file(filename)
    p = figure(
        # title="Data vs. Gamma Fits (Original & Adjusted)",
        x_axis_label="Turn-shift Speed",
        y_axis_label="Probability Density",
        width=800,
        height=400,
    )

    # ラベルと目盛りのフォントサイズを設定
    p.xaxis.axis_label_text_font_size = "16pt"
    p.yaxis.axis_label_text_font_size = "16pt"
    p.xaxis.major_label_text_font_size = "14pt"
    p.yaxis.major_label_text_font_size = "14pt"
    p.legend.label_text_font_size = "16pt"

    # 元データ
    p.quad(
        top=hist,
        bottom=0,
        left=edges[:-1],
        right=edges[1:],
        fill_color="navy",
        line_color="white",
        alpha=0.5,
        legend_label="Original Data",
    )
    # オリジナルフィット
    p.line(
        x,
        base_pdf,
        line_width=2,
        line_dash="dashed",
        legend_label="Gamma Fit (Original)",
    )
    # アジャスト後フィット
    p.line(
        x,
        adj_pdf,
        line_width=2,
        legend_label="Gamma Fit (Adjusted)",
    )

    p.legend.location = "top_right"
    save(p)


def adjust_gamma_params(fitted_params, desired_mean, desired_variance):
    """
    Adjusts the parameters of a gamma distribution to match a desired mean and variance.

    Args:
    - fitted_params (dict): A dictionary containing the current 'shape', 'loc', and 'scale' of the distribution.
    - desired_mean (float): The target mean value.
    - desired_variance (float): The target variance value.

    Returns:
    - adjusted_params (dict): A dictionary containing the adjusted 'shape', 'loc', and 'scale' parameters.
    """
    shape = fitted_params["shape"]
    loc = fitted_params["loc"]
    scale = fitted_params["scale"]

    # Adjust scale and shape to match the desired mean and variance
    # adjusted_scale = np.sqrt(desired_variance / shape)
    # adjusted_shape = (desired_mean - loc) / adjusted_scale

    adjusted_scale = desired_variance / (desired_mean - loc)
    adjusted_shape = (desired_mean - loc) ** 2 / desired_variance

    # Return the new parameters
    adjusted_params = {
        "shape": adjusted_shape,
        "loc": loc,  # loc remains unchanged
        "scale": adjusted_scale,
    }

    return adjusted_params


def evaluate_gamma_fit_r2(turn_shifts, params, label="", bins=50):
    """
    - turn_shifts: データ配列
    - params: {'shape', 'loc', 'scale'}
    - bins: ヒストグラムのビン数
    → R² を返す (0–1 の値。*100 で％表示)
    """
    hist, edges = np.histogram(turn_shifts, bins=bins, density=True)
    mids = (edges[:-1] + edges[1:]) / 2
    mask = mids >= params["loc"]
    hist, mids = hist[mask], mids[mask]

    # 1) 実データの密度推定（ヒストグラムを確率密度化）
    hist, edges = np.histogram(turn_shifts, bins=bins, density=True)
    mids = (edges[:-1] + edges[1:]) / 2  # ビンの中心値

    # 2) モデルの PDF を同点で評価
    pdf_vals = gamma.pdf(
        mids,
        params["shape"],
        loc=params["loc"],
        scale=params["scale"],
    )

    # 3) 決定係数 R² の計算
    ss_res = np.sum((hist - pdf_vals) ** 2)
    ss_tot = np.sum((hist - np.mean(hist)) ** 2)
    r2 = 1 - ss_res / ss_tot
    logger.info(f"{label} gamma R²: {r2 * 100:.2f}%")
    return r2


def evaluate_gamma_fit_hellinger_overlap(data, params, n=1000, label=""):
    x = np.linspace(min(data), max(data), n)
    # KDE で経験分布をスムージング
    p_emp = gaussian_kde(data)(x)
    p_model = gamma.pdf(x, params["shape"], loc=params["loc"], scale=params["scale"])
    p_emp /= np.trapz(p_emp, x)  # 正規化
    p_model /= np.trapz(p_model, x)  # 正規化
    overlap = np.trapz(np.sqrt(p_emp * p_model), x)
    logger.info(f"{label} Hellinger一致率: {overlap * 100:.2f}%")
    return overlap * 100  # 0–100 %


# 2つのグループとそのガンマ分布を同一グラフに描画する関数
def plot_two_distributions_with_gamma(
    shift1,
    shift2,
    label1,
    label2,
    title,
    base_params,
    adjusted_params1,
    adjusted_params2,
    adjusted_label1,
    adjusted_label2,
):
    bins = np.arange(-4.0, 4.2, 0.2)
    hist1, edges = np.histogram(shift1, bins=bins, density=True)
    hist2, _ = np.histogram(shift2, bins=bins, density=True)

    p = figure(
        title=title,
        x_axis_label="Turn Shift Duration (s)",
        y_axis_label="Probability Density",
        width=800,
        height=600,
        x_range=(-4, 4),
    )

    # ラベルと目盛りのフォントサイズを設定
    p.xaxis.axis_label_text_font_size = "16pt"
    p.yaxis.axis_label_text_font_size = "16pt"
    p.xaxis.major_label_text_font_size = "14pt"
    p.yaxis.major_label_text_font_size = "14pt"
    p.legend.label_text_font_size = "16pt"

    # ヒストグラム
    p.quad(
        top=hist1,
        bottom=0,
        left=edges[:-1],
        right=edges[1:],
        fill_color="navy",
        alpha=0.5,
        legend_label=label1,
    )
    p.quad(
        top=hist2,
        bottom=0,
        left=edges[:-1],
        right=edges[1:],
        fill_color="orange",
        alpha=0.5,
        legend_label=label2,
    )

    # ガンマ分布のオーバーレイ
    x = np.linspace(-4, 4, 1000)
    base_pdf = gamma.pdf(
        x,
        base_params["shape"],
        loc=base_params["loc"],
        scale=base_params["scale"],
    )
    pdf_high = gamma.pdf(
        x,
        adjusted_params1["shape"],
        loc=adjusted_params1["loc"],
        scale=adjusted_params1["scale"],
    )
    pdf_low = gamma.pdf(
        x,
        adjusted_params2["shape"],
        loc=adjusted_params2["loc"],
        scale=adjusted_params2["scale"],
    )
    # ベース分布を破線、グループ分布を実線で
    p.line(x, base_pdf, line_dash="dashed", line_width=2, legend_label="Base Gamma")
    p.line(x, pdf_high, line_width=2, line_color="red", legend_label=adjusted_label1)
    p.line(x, pdf_low, line_width=2, line_color="green", legend_label=adjusted_label2)

    p.legend.location = "top_right"
    return p


# memo
filled_params = {
    "shape": np.float64(9.845759954342185),
    "loc": np.float64(-1.7494952473434728),
    "scale": np.float64(0.21313221400417537),
}

if __name__ == "__main__":
    """
    全体でのターンシフトの分布をガンマ分布にフィットさせる。
    NOTE: ここのコードは実装の音声が必要です。実際の音声が必要な方は著者に連絡してください。
    """
    # turn_shifts = analysis_all()
    # filled_params = fit_gamma_distribution(turn_shifts)
    # evaluate_gamma_fit(turn_shifts, filled_params, label="Based")
    # evaluate_gamma_fit_r2(turn_shifts, filled_params, label="Based")
    # evaluate_gamma_fit_hellinger_overlap(
    #     turn_shifts, filled_params, n=1000, label="Based"
    # )
    # logger.info(filled_params)
    # plot_gamma_fit_html(turn_shifts, filled_params)

    """
    ExpertとNoviceのターンシフトの分布を平均値と分散からガンマ分布にフィットさせる。
    """
    expert_turn_shifts, novice_turn_shifts = analysis_speaker()

    adjust_params_expert = adjust_gamma_params(
        filled_params,
        desired_mean=np.mean(expert_turn_shifts),
        desired_variance=np.var(expert_turn_shifts),
    )
    logger.info(adjust_params_expert)
    evaluate_gamma_fit(expert_turn_shifts, filled_params, label="Expert Based")
    evaluate_gamma_fit(
        expert_turn_shifts, adjust_params_expert, label="Expert Adjusted"
    )
    evaluate_gamma_fit_r2(expert_turn_shifts, filled_params, label="Expert Based")
    evaluate_gamma_fit_r2(
        expert_turn_shifts, adjust_params_expert, label="Expert Adjusted"
    )
    evaluate_gamma_fit_hellinger_overlap(
        expert_turn_shifts, filled_params, n=1000, label="Expert Based"
    )
    evaluate_gamma_fit_hellinger_overlap(
        expert_turn_shifts, adjust_params_expert, n=1000, label="Expert Adjusted"
    )
    plot_gamma_fit_comparison_html(
        expert_turn_shifts,
        filled_params,
        adjust_params_expert,
        filename="output/gamma_fit_expert.html",
    )

    adjust_params_novice = adjust_gamma_params(
        filled_params,
        desired_mean=np.mean(novice_turn_shifts),
        desired_variance=np.var(novice_turn_shifts),
    )
    logger.info(adjust_params_novice)
    evaluate_gamma_fit(novice_turn_shifts, filled_params, label="Novice Based")
    evaluate_gamma_fit(
        novice_turn_shifts, adjust_params_novice, label="Novice Adjusted"
    )
    evaluate_gamma_fit_r2(novice_turn_shifts, filled_params, label="Novice Based")
    evaluate_gamma_fit_r2(
        novice_turn_shifts, adjust_params_novice, label="Novice Adjusted"
    )
    evaluate_gamma_fit_hellinger_overlap(
        novice_turn_shifts, filled_params, n=1000, label="Novice Based"
    )
    evaluate_gamma_fit_hellinger_overlap(
        novice_turn_shifts, adjust_params_novice, n=1000, label="Novice Adjusted"
    )
    plot_gamma_fit_comparison_html(
        novice_turn_shifts,
        filled_params,
        adjust_params_novice,
        filename="output/gamma_fit_novice.html",
    )

    """
    friendとstrangerのターンシフトの分布を平均値と分散からガンマ分布にフィットさせる。
    """
    friend_turn_shifts, stranger_turn_shifts = analysis_relationship()

    adjust_params_friend = adjust_gamma_params(
        filled_params,
        desired_mean=np.mean(friend_turn_shifts),
        desired_variance=np.var(friend_turn_shifts),
    )
    logger.info(adjust_params_friend)
    evaluate_gamma_fit(friend_turn_shifts, filled_params, label="Friend Based")
    evaluate_gamma_fit(
        friend_turn_shifts, adjust_params_friend, label="Friend Adjusted"
    )
    evaluate_gamma_fit_r2(friend_turn_shifts, filled_params, label="Friend Based")
    evaluate_gamma_fit_r2(
        friend_turn_shifts, adjust_params_friend, label="Friend Adjusted"
    )
    evaluate_gamma_fit_hellinger_overlap(
        friend_turn_shifts, filled_params, n=1000, label="Friend Based"
    )
    evaluate_gamma_fit_hellinger_overlap(
        friend_turn_shifts, adjust_params_friend, n=1000, label="Friend Adjusted"
    )
    plot_gamma_fit_comparison_html(
        friend_turn_shifts,
        filled_params,
        adjust_params_friend,
        filename="output/gamma_fit_friend.html",
    )
    adjust_params_stranger = adjust_gamma_params(
        filled_params,
        desired_mean=np.mean(stranger_turn_shifts),
        desired_variance=np.var(stranger_turn_shifts),
    )
    logger.info(adjust_params_stranger)
    evaluate_gamma_fit(stranger_turn_shifts, filled_params, label="Stranger Based")
    evaluate_gamma_fit(
        stranger_turn_shifts, adjust_params_stranger, label="Stranger Adjusted"
    )
    evaluate_gamma_fit_r2(stranger_turn_shifts, filled_params, label="Stranger Based")
    evaluate_gamma_fit_r2(
        stranger_turn_shifts, adjust_params_stranger, label="Stranger Adjusted"
    )
    evaluate_gamma_fit_hellinger_overlap(
        stranger_turn_shifts, filled_params, n=1000, label="Stranger Based"
    )
    evaluate_gamma_fit_hellinger_overlap(
        stranger_turn_shifts, adjust_params_stranger, n=1000, label="Stranger Adjusted"
    )
    plot_gamma_fit_comparison_html(
        stranger_turn_shifts,
        filled_params,
        adjust_params_stranger,
        filename="output/gamma_fit_stranger.html",
    )

    """
    個人特性ごとのターンシフトの分布をガンマ分布にフィットさせる。
    """
    (
        high_extroversion,
        low_extroversion,
        high_openness,
        low_openness,
        high_agreeableness,
        low_agreeableness,
        high_conscientiousness,
        low_conscientiousness,
        high_neuroticism,
        low_neuroticism,
    ) = analysis_pesonal()

    # high_extroversion などのリストは既に定義済みとする
    output_file("output/turn_shift_statistics_BIG5.html")
    plots = []

    # 各特性ごとのデータ
    traits = [
        ("Extroversion", high_extroversion, low_extroversion),
        ("Openness", high_openness, low_openness),
        ("Agreeableness", high_agreeableness, low_agreeableness),
        ("Conscientiousness", high_conscientiousness, low_conscientiousness),
        ("Neuroticism", high_neuroticism, low_neuroticism),
    ]

    for name, high_group, low_group in traits:
        # 高群と低群それぞれの平均・分散に合わせてガンマ分布を調整
        adj_params_high = adjust_gamma_params(
            filled_params,
            desired_mean=np.mean(high_group),
            desired_variance=np.var(high_group),
        )
        adj_params_low = adjust_gamma_params(
            filled_params,
            desired_mean=np.mean(low_group),
            desired_variance=np.var(low_group),
        )
        logger.info(f"Adjusted params for {name} high: {adj_params_high}")
        logger.info(f"Adjusted params for {name} low: {adj_params_low}")

        plots.append(
            plot_two_distributions_with_gamma(
                high_group,
                low_group,
                f"High {name}",
                f"Low {name}",
                f"Turn Shift & Gamma - {name}",
                filled_params,
                adj_params_high,
                adj_params_low,
                f"High {name} Gamma",
                f"Low {name} Gamma",
            )
        )

        evaluate_gamma_fit(high_group, filled_params, label=f"High {name} Based")
        evaluate_gamma_fit(high_group, adj_params_high, label=f"High {name} Adjusted")
        evaluate_gamma_fit(low_group, filled_params, label=f"Low {name} Based")
        evaluate_gamma_fit(low_group, adj_params_low, label=f"Low {name} Adjusted")
        evaluate_gamma_fit_r2(high_group, filled_params, label=f"High {name} Based")
        evaluate_gamma_fit_r2(
            high_group, adj_params_high, label=f"High {name} Adjusted"
        )
        evaluate_gamma_fit_r2(low_group, filled_params, label=f"Low {name} Based")
        evaluate_gamma_fit_r2(low_group, adj_params_low, label=f"Low {name} Adjusted")
        evaluate_gamma_fit_hellinger_overlap(
            high_group, filled_params, n=1000, label=f"High {name} Based"
        )
        evaluate_gamma_fit_hellinger_overlap(
            high_group, adj_params_high, n=1000, label=f"High {name} Adjusted"
        )
        evaluate_gamma_fit_hellinger_overlap(
            low_group, filled_params, n=1000, label=f"Low {name} Based"
        )
        evaluate_gamma_fit_hellinger_overlap(
            low_group, adj_params_low, n=1000, label=f"Low {name} Adjusted"
        )
    save(column(*plots))