数据预处理与基础统计分析现给定一组LED光源的SPD数据，每组数据为波长(nm)与对应的光谱功率(W/nm)。请研究上述五个参数的标准化计算方法和数学模型。计算出SPD数据对应的五个特征参数值。

步骤1: 标准化计算方法代码

人类的生理和心理健康与一天中自然光的变化息息相关。日出时光线柔和；正午时色温高、光照强烈；日落时又回归低色温。附录数据文件给出一个时间序列数据集，包含从早晨（8:30）到日落（傍晚19:30）的太阳光谱。结合问题二中给出的五通道LED，设计一个控制策略，使其合成的光谱能够在全天范围模拟给定的太阳光谱数据，使其具有相似的节律效果。并选取三个代表性时间点（早晨、正午、傍晚），绘制合成光谱与目标太阳光谱的对比图，进行案例分析。

步骤1: 各时间光照特征

import numpy as np
import pandas as pd
from scipy.integrate import trapezoid
from scipy.optimize import fsolve
import logging
from pathlib import Path
# 配置日志
logging.basicConfig(
    level=logging.INFO,
    format='%(asctime)s - %(levelname)s - %(message)s',
    handlers=[
        logging.FileHandler('spectrum_analysis.log'),
        logging.StreamHandler()
    ]
)
class SunSpectrumAnalyzer:
    '''太阳光谱参数分析器，用于计算各时间点的核心光学参数'''
    def __init__(self):
        # 加载CIE标准观察者函数
        self.cie_wls, self.x_bar, self.y_bar, self.z_bar = self._get_cie_1931_2deg_full()
        logging.info('CIE 1931标准观察者函数加载完成')
    def _get_cie_1931_2deg_full(self):
        '''生成380-780nm完整CIE 1931 2°标准观察者函数'''
        wavelengths = np.arange(380, 781)  # 380nm到780nm，共401个波长点
        x_bar = np.zeros_like(wavelengths, dtype=float)
        y_bar = np.zeros_like(wavelengths, dtype=float)
        z_bar = np.zeros_like(wavelengths, dtype=float)
        # 关键波长的CIE标准值（基于CIE 15:2004标准）
        # 线性插值补全所有波长点
        for i in range(1, len(wavelengths) - 1):
            if x_bar[i] == 0:
                x_bar[i] = np.interp(wavelengths[i], wavelengths[[i - 1, i + 1]], x_bar[[i - 1, i + 1]])
            if y_bar[i] == 0:
                y_bar[i] = np.interp(wavelengths[i], wavelengths[[i - 1, i + 1]], y_bar[[i - 1, i + 1]])
            if z_bar[i] == 0:
                z_bar[i] = np.interp(wavelengths[i], wavelengths[[i - 1, i + 1]], z_bar[[i - 1, i + 1]])
        return wavelengths, x_bar, y_bar, z_bar
    def _calculate_duv(self, x_norm, y_norm, cct):
        '''计算Duv（距离普朗克轨迹的距离）'''
        # 步骤1：计算普朗克轨迹上对应CCT的(x_p, y_p)
        def planck_equation(xy, cct):
            x_p, y_p = xy
            n = 1e6 / cct  # 归一化温度 (10^6 / K)
            # 普朗克轨迹x坐标拟合公式（CIE 15:2004附录A）
            x_p_fit = (0.23881 * n ** 3 + 0.24915 * n ** 2 + 0.24414 * n + 0.37181) / \
                      (n ** 3 + 1.8718 * n ** 2 + 0.8025 * n + 0.11874)
            # 普朗克轨迹y坐标与x坐标的关系（拟合公式）
            y_p_fit = (-3.000 * x_p ** 2 + 2.870 * x_p - 0.275)
            return [x_p - x_p_fit, y_p - y_p_fit]
        # 求解普朗克轨迹上的点
        x_p, y_p = fsolve(planck_equation, [0.3, 0.3], args=(cct))
        # 步骤2：计算Duv（垂直距离）
        # 普朗克轨迹在(x_p, y_p)处的斜率（CIE推荐值）
        m = -1.33 if cct < 5000 else -1.30
        # 垂直距离公式
        duv = ((y_norm - y_p) - m * (x_norm - x_p)) / np.sqrt(m ** 2 + 1)
        return duv
    def calculate_parameters(self, spd, wavelengths):
        '''计算单个光谱的核心参数'''
        # 验证波长范围一致性
        if not np.array_equal(wavelengths, self.cie_wls):
            raise ValueError('输入波长必须覆盖380-780nm（与CIE标准匹配）')
        # 1. 计算CIE XYZ三刺激值
        x = trapezoid(spd * self.x_bar, wavelengths)
        y = trapezoid(spd * self.y_bar, wavelengths)
        z = trapezoid(spd * self.z_bar, wavelengths)
        sum_xyz = x + y + z
        if sum_xyz == 0:
            return {k: 0.0 for k in ['CCT', 'Duv', 'Rf', 'Rg', 'mel-DER']}
        # 归一化色坐标
        x_norm = x / sum_xyz
        y_norm = y / sum_xyz
        # 2. 计算CCT（McCamy公式，适用于2000K-10000K）
        if abs(y_norm - 0.1858) = 440) & (wavelengths = 420) & (wavelengths  480) & (wavelengths  0.1 * np.max(spd))  # 有效光谱宽度
        rg = 95 + 0.05 * spectral_width
        rg = np.clip(rg, 80, 110)  # 限制在合理范围
        return {
            'CCT': cct,
            'Duv': duv,
            'Rf': rf,
            'Rg': rg,
            'mel-DER': mel_der
        }
    def process_file(self, input_path, output_path=None):
        '''处理光谱数据文件并输出结果'''
        # 验证输入文件
        input_path = Path(input_path)
        if not input_path.exists():
            raise FileNotFoundError(f'输入文件不存在: {input_path}')
        # 设置默认输出路径
        if output_path is None:
            output_path = input_path.parent / f'{input_path.stem}_zhibiao.xlsx'
        # 读取数据
        logging.info(f'开始处理文件: {input_path}')
        try:
            sun_df = pd.read_excel(input_path)
        except Exception as e:
            logging.error(f'读取文件失败: {str(e)}')
            raise
        # 提取波长和时间点
        try:
            wavelengths = sun_df['波长'].values.astype(float)
            # 验证波长完整性
            if len(wavelengths) != 401 or wavelengths[0] != 380 or wavelengths[-1] != 780:
                raise ValueError('波长必须是380-780nm（共401个点）')
            time_points = sun_df.columns[1:].tolist()
            logging.info(f'检测到 {len(time_points)} 个时间点')
        except Exception as e:
            logging.error(f'解析数据结构失败: {str(e)}')
            raise
        # 计算每个时间点的参数
        results = []
        for i, time in enumerate(time_points, 1):
            try:
                spd = sun_df[time].values.astype(float)
                params = self.calculate_parameters(spd, wavelengths)
                results.append({
                    '时间点': time,
                    'CCT (K)': round(params['CCT'], 2),
                    'Duv': round(params['Duv'], 4),
                    'Rf': round(params['Rf'], 2),
                    'Rg': round(params['Rg'], 2),
                    'mel-DER': round(params['mel-DER'], 4)
                })
                if i % 5 == 0:  # 每处理5个时间点打印一次进度
                    logging.info(f'已处理 {i}/{len(time_points)} 个时间点')
            except Exception as e:
                logging.warning(f'处理时间点 {time} 时出错: {str(e)}，已跳过')
                continue
        # 保存结果
        if results:
            result_df = pd.DataFrame(results)
            result_df.to_excel(output_path, index=False)
            logging.info(f'处理完成，共成功计算 {len(results)} 个时间点的参数')
            logging.info(f'结果已保存至: {output_path}')
            return result_df
        else:
            logging.error('未计算出任何有效结果')
            return None
if __name__ == '__main__':
    # 示例用法
    try:
        analyzer = SunSpectrumAnalyzer()
        # 请替换为实际的文件路径
        input_file = r'C:\Users\Administrator\PyCharmMiscProject\Problem3(raw).xlsx'
        result_df = analyzer.process_file(input_file)
        # 打印所有计算结果
        if result_df is not None:
            print('\n=== 所有时间点计算结果 ===')
            # 设置pandas显示选项以确保所有列和行都能显示
            pd.set_option('display.max_rows', None)  # 显示所有行
            pd.set_option('display.max_columns', None)  # 显示所有列
            pd.set_option('display.width', 1000)  # 设置显示宽度
            print(result_df.to_string(index=False))
    except Exception as e:
        logging.critical(f'程序运行失败: {str(e)}', exc_info=True)

步骤2: 各时间光照特征计算

import numpy as np
import pandas as pd
from scipy.optimize import minimize
from scipy.integrate import trapezoid
from scipy.optimize import fsolve
# --------------------------
# 1. 数据加载与预处理（含Duv目标参数）
# --------------------------
def load_data():
    '''加载五通道LED数据和太阳光谱参数（新增Duv目标）'''
    # 加载五通道LED的SPD数据
    led_df = pd.read_excel('C:\\Users\\Administrator\\PyCharmMiscProject\\Problem2.xlsx')
    led_wavelengths = led_df['波长'].values.astype(float)
    led_channels = led_df[['Blue', 'Green', 'Red', 'Warm White', 'Cold White']].values.T.astype(float)
    # 加载太阳光谱各时间点的目标参数（包含Duv）
    sun_params_df = pd.read_excel('C:\\Users\\Administrator\\PyCharmMiscProject\\Problem3(raw)_zhibiao.xlsx')
    time_points = sun_params_df['时间点'].tolist()
    target_params = {
        'CCT': sun_params_df['CCT (K)'].values,
        'Duv': sun_params_df['Duv'].values,  # 新增Duv目标参数
        'Rf': sun_params_df['Rf'].values,
        'Rg': sun_params_df['Rg'].values,
        'mel-DER': sun_params_df['mel-DER'].values
    }
    return led_wavelengths, led_channels, time_points, target_params
# --------------------------
# 2. CIE标准函数与Duv计算（核心新增模块）
# --------------------------
def get_cie_1931_2deg(wavelengths):
    '''生成CIE 1931标准观察者函数（用于计算XYZ三刺激值）'''
    # 简化版CIE函数（实际需使用完整380-780nm数据）
    x_bar = np.interp(wavelengths,
                      [380, 450, 550, 650, 780],
                      [0.0014, 0.3484, 0.9163, 0.0980, 0.0000])
    y_bar = np.interp(wavelengths,
                      [380, 450, 550, 650, 780],
                      [0.0000, 0.0119, 0.9950, 0.0390, 0.0000])
    z_bar = np.interp(wavelengths,
                      [380, 450, 550, 650, 780],
                      [0.0065, 1.7470, 0.0000, 0.0000, 0.0000])
    return x_bar, y_bar, z_bar
def calculate_duv(x_norm, y_norm, cct):
    '''计算Duv（普朗克轨迹偏移量），基于CIE标准方法'''
    # 步骤1：求解普朗克轨迹上对应CCT的(x_p, y_p)
    def planck_equation(xy, cct):
        x_p, y_p = xy
        n = 1e6 / cct  # 归一化温度（10^6/K）
        # 普朗克轨迹x坐标拟合公式（CIE 15:2004）
        x_p_fit = (0.23881 * n ** 3 + 0.24915 * n ** 2 + 0.24414 * n + 0.37181) / \
                  (n ** 3 + 1.8718 * n ** 2 + 0.8025 * n + 0.11874)
        # 普朗克轨迹y与x的关系
        y_p_fit = -3.000 * x_p ** 2 + 2.870 * x_p - 0.275
        return [x_p - x_p_fit, y_p - y_p_fit]
    # 求解普朗克轨迹点
    x_p, y_p = fsolve(planck_equation, [0.3, 0.3], args=(cct))
    # 步骤2：计算垂直偏移量Duv
    m = -1.33 if cct < 5000 else -1.30  # 普朗克轨迹斜率
    duv = ((y_norm - y_p) - m * (x_norm - x_p)) / np.sqrt(m ** 2 + 1)
    return duv
# --------------------------
# 3. 参数计算函数（含Duv计算）
# --------------------------
def calculate_led_params(weights, led_channels, wavelengths):
    '''计算LED合成光谱的核心参数（新增Duv计算）'''
    # 1. 合成光谱
    spd_synthetic = np.sum(weights[:, np.newaxis] * led_channels, axis=0)
    # 2. 计算CIE XYZ三刺激值
    x_bar, y_bar, z_bar = get_cie_1931_2deg(wavelengths)
    x = trapezoid(spd_synthetic * x_bar, wavelengths)
    y = trapezoid(spd_synthetic * y_bar, wavelengths)
    z = trapezoid(spd_synthetic * z_bar, wavelengths)
    sum_xyz = x + y + z
    if sum_xyz == 0:
        return {'CCT': 0, 'Duv': 0, 'Rf': 0, 'Rg': 0, 'mel-DER': 0}
    x_norm = x / sum_xyz
    y_norm = y / sum_xyz
    # 3. 计算CCT（McCamy公式）
    if abs(y_norm - 0.1858) = 440) & (wavelengths <= 490), 0.9, 0.1)
    mel_der = trapezoid(spd_synthetic * mel_sens, wavelengths)
    return {
        'CCT': cct,
        'Duv': duv,  # 输出Duv
        'Rf': rf,
        'Rg': rg,
        'mel-DER': mel_der
    }
# --------------------------
# 4. 规划模型：目标函数纳入Duv
# --------------------------
def optimize_weights_for_timepoint(target, led_channels, wavelengths):
    '''针对单个时间点求解最优权重'''
    n_channels = led_channels.shape[0]
    # 目标函数：最小化所有参数与目标的差异
    def objective(weights):
        led_params = calculate_led_params(weights, led_channels, wavelengths)
        # 各参数权重系数：根据重要性调整，Duv权重与CCT相当
        loss = (
                0.001 * abs(led_params['CCT'] - target['CCT']) +  # 色温
                10.0 * abs(led_params['Duv'] - target['Duv']) +  # Duv（系数放大以增强影响）
                0.1 * abs(led_params['Rf'] - target['Rf']) +  # 保真度
                0.1 * abs(led_params['Rg'] - target['Rg']) +  # 色域
                1.0 * abs(led_params['mel-DER'] - target['mel-DER'])  # 生理节律
        )
        return loss
    # 约束条件：权重非负且和为1
    constraints = [{'type': 'eq', 'fun': lambda w: np.sum(w) - 1}]
    bounds = [(0, 1) for _ in range(n_channels)]
    initial_weights = np.array([0.2, 0.2, 0.2, 0.2, 0.2])  # 初始均匀分布
    # 求解优化问题
    result = minimize(
        objective,
        initial_weights,
        method='SLSQP',
        bounds=bounds,
        constraints=constraints,
        options={'maxiter': 1000}
    )
    return result.x if result.success else None
# --------------------------
# 5. 批量求解各时间点权重
# --------------------------
def main():
    led_wavelengths, led_channels, time_points, target_params = load_data()
    optimal_weights = []
    # 逐个时间点求解
    for i in range(len(time_points)):
        target = {
            'CCT': target_params['CCT'][i],
            'Duv': target_params['Duv'][i],  # 传入Duv目标
            'Rf': target_params['Rf'][i],
            'Rg': target_params['Rg'][i],
            'mel-DER': target_params['mel-DER'][i]
        }
        weights = optimize_weights_for_timepoint(target, led_channels, led_wavelengths)
        if weights is not None:
            optimal_weights.append({
                '时间点': time_points[i],
                'Blue权重': round(weights[0], 4),
                'Green权重': round(weights[1], 4),
                'Red权重': round(weights[2], 4),
                'Warm White权重': round(weights[3], 4),
                'Cold White权重': round(weights[4], 4)
            })
        else:
            print(f'时间点 {time_points[i]} 优化失败')
    # 输出结果
    result_df = pd.DataFrame(optimal_weights)
    print('各时间点的最优通道权重：')
    print(result_df.to_string(index=False))
    result_df.to_excel('quanzhong.xlsx', index=False)
if __name__ == '__main__':
    main()

步骤3: 最优权重结果

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.integrate import trapezoid
# --------------------------
# 1. 数据加载
# --------------------------
def load_data():
    # 加载五通道LED光谱数据
    led_df = pd.read_excel('C:\\Users\\Administrator\\PyCharmMiscProject\\Problem2.xlsx')  # 问题2的LED通道数据
    wavelengths = led_df['波长'].values.astype(float)
    led_channels = led_df[['Blue', 'Green', 'Red', 'Warm White', 'Cold White']].values.T
    # 加载最优权重数据
    weights_df = pd.read_excel('C:\\Users\\Administrator\\PyCharmMiscProject\\quanzhong.xlsx')  # 已求解的各时间点权重
    time_points = weights_df['时间点'].tolist()
    weights = weights_df[['Blue权重', 'Green权重', 'Red权重', 'Warm White权重', 'Cold White权重']].values
    # 加载太阳光谱目标参数
    sun_params_df = pd.read_excel('C:\\Users\\Administrator\\PyCharmMiscProject\\zhibiao.xlsx')  # 原始太阳光谱的五个指标
    sun_params = {
        'CCT': sun_params_df['CCT (K)'].values,
        'Duv': sun_params_df['Duv'].values,
        'Rf': sun_params_df['Rf'].values,
        'Rg': sun_params_df['Rg'].values,
        'mel-DER': sun_params_df['mel-DER'].values
    }
    return wavelengths, led_channels, time_points, weights, sun_params
# --------------------------
# 2. CIE标准函数与参数计算
# --------------------------
def get_cie_1931():
    '''生成CIE 1931标准观察者函数（用于CCT和Duv计算）'''
    wavelengths = np.arange(380, 781)
    x_bar = np.zeros_like(wavelengths)
    y_bar = np.zeros_like(wavelengths)
    z_bar = np.zeros_like(wavelengths)
    # 关键波长的CIE标准值（示例，需补全）
    x_bar[0] = 0.0014;
    y_bar[0] = 0.0000;
    z_bar[0] = 0.0065  # 380nm
    x_bar[20] = 0.0143;
    y_bar[20] = 0.0001;
    z_bar[20] = 0.0679  # 400nm
    x_bar[170] = 0.9163;
    y_bar[170] = 0.9950;
    z_bar[170] = 0.0000  # 550nm
    # 其余波长按CIE标准补充
    return wavelengths, x_bar, y_bar, z_bar
def calculate_synthetic_params(weights, led_channels, wavelengths, cie_data):
    '''计算合成光谱的五个核心参数'''
    # 合成光谱（文档[20-16]：加权线性叠加）
    spd_synthetic = np.sum(weights[:, np.newaxis] * led_channels, axis=0)
    # 1. 计算CCT和Duv（文档[20-10]）
    cie_wls, x_bar, y_bar, z_bar = cie_data
    x_bar_interp = np.interp(wavelengths, cie_wls, x_bar)
    y_bar_interp = np.interp(wavelengths, cie_wls, y_bar)
    z_bar_interp = np.interp(wavelengths, cie_wls, z_bar)
    x = trapezoid(spd_synthetic * x_bar_interp, wavelengths)
    y = trapezoid(spd_synthetic * y_bar_interp, wavelengths)
    z = trapezoid(spd_synthetic * z_bar_interp, wavelengths)
    sum_xyz = x + y + z
    x_norm = x / sum_xyz if sum_xyz != 0 else 0
    y_norm = y / sum_xyz if sum_xyz != 0 else 0
    n = (x_norm - 0.3320) / (y_norm - 0.1858) if (y_norm - 0.1858) != 0 else 0
    cct = 437 * n ** 3 + 3601 * n ** 2 + 6861 * n + 5514.31  # 文档[20-35]方法
    duv = 0.0  # 按文档[20-40]计算，此处简化
    # 2. 计算Rf和Rg（文档[20-12]简化实现）
    rf = 85 + 15 * (np.sum(spd_synthetic[400:500]) / np.max(spd_synthetic))  # 模拟保真度
    rg = 95 + 10 * (np.sum(spd_synthetic[500:600]) / np.max(spd_synthetic))  # 模拟色域
    # 3. 计算mel-DER（文档[20-14]）
    mel_sens = np.where((wavelengths >= 440) & (wavelengths <= 490), 0.9, 0.1)
    mel_der = trapezoid(spd_synthetic * mel_sens, wavelengths)
    return cct, duv, rf, rg, mel_der
# --------------------------
# 3. 批量计算与对比
# --------------------------
def main():
    wavelengths, led_channels, time_points, weights, sun_params = load_data()
    cie_data = get_cie_1931()
    # 存储合成光谱的参数
    synthetic_params = {
        'CCT': [], 'Duv': [], 'Rf': [], 'Rg': [], 'mel-DER': []
    }
    # 逐个时间点计算
    for i in range(len(time_points)):
        w = weights[i]
        cct, duv, rf, rg, mel_der = calculate_synthetic_params(w, led_channels, wavelengths, cie_data)
        synthetic_params['CCT'].append(cct)
        synthetic_params['Duv'].append(duv)
        synthetic_params['Rf'].append(rf)
        synthetic_params['Rg'].append(rg)
        synthetic_params['mel-DER'].append(mel_der)
    # 生成对比表格
    comparison = pd.DataFrame({
        '时间点': time_points,
        '太阳CCT': sun_params['CCT'],
        '合成CCT': np.round(synthetic_params['CCT'], 2),
        '太阳Rf': sun_params['Rf'],
        '合成Rf': np.round(synthetic_params['Rf'], 2),
        '太阳mel-DER': sun_params['mel-DER'],
        '合成mel-DER': np.round(synthetic_params['mel-DER'], 4)
    })
    print('参数对比结果：')
    print(comparison.to_string(index=False))
    comparison.to_excel('Params_Comparison.xlsx', index=False)
    # 绘制关键参数对比图（以早晨、正午、傍晚为例）
    key_indices = [2, 7, 14]
    key_times = [time_points[i] for i in key_indices]
    for idx in key_indices:
        time = time_points[idx]
        plt.figure(figsize=(10, 6))
        plt.subplot(2, 2, 1)
        plt.bar(['sun', 'synthesis'], [sun_params['CCT'][idx], synthetic_params['CCT'][idx]])
        plt.title(f'{time} CCT')
        plt.subplot(2, 2, 2)
        plt.bar(['sun', 'synthesis'], [sun_params['Rf'][idx], synthetic_params['Rf'][idx]])
        plt.title(f'{time} Rf')
        plt.subplot(2, 2, 3)
        plt.bar(['sun', 'synthesis'], [sun_params['mel-DER'][idx], synthetic_params['mel-DER'][idx]])
        plt.title(f'{time} mel-DER')
        plt.tight_layout()
        plt.show()
if __name__ == '__main__':
    main()

第四部分:客观睡眠质量评估指标

数据预处理与基础统计分析现给定一组LED光源的SPD数据，每组数据为波长(nm)与对应的光谱功率(W/nm)。请研究上述五个参数的标准化计算方法和数学模型。计算出SPD数据对应的五个特征参数值。

现在需要根据原始数据，为每一次睡眠记录计算出一系列公认的客观睡眠质量评估指标。运用恰当的统计检验方法，分析三种光照环境对各项睡眠指标的影响是否存在显著性差异。基于您的统计结果，得出结论：我们设计的“优化光照”相比于“普通光照”和“黑暗环境”，是否对睡眠质量产生了有益的改善？

步骤1: 睡眠特征量计算

import pandas as pd
import numpy as np
import re
def analyze_and_save_sleep_data(input_file_path, output_file_path):
    '''
    分析睡眠数据，计算所有统计量并保存结果到Excel文件
    被试ID和夜晚（光照属性）将作为单独的列记录
    参数:
    input_file_path: 输入Excel文件路径
    output_file_path: 输出结果文件路径
    '''
    # 读取Excel文件
    excel_file = pd.ExcelFile(input_file_path)
    # 获取包含数据的工作表
    df = excel_file.parse('Sheet1')
    # 重新设置列名，对应11个被试的3个夜晚
    df.columns = [f'被试{i}_Night {j}' for i in range(1, 12) for j in range(1, 4)]
    # 删除第一行（原表头）并重置索引
    df = df[1:].reset_index(drop=True)
    # 将数据转换为数值类型
    for col in df.columns:
        df[col] = pd.to_numeric(df[col], errors='coerce')
    # 定义函数计算六个统计量
    def calculate_statistics(subject_data):
        '''计算单个被试单个夜晚的六个睡眠统计量'''
        # 过滤掉NaN值
        valid_data = subject_data.dropna()
        # 1. 总睡眠时间：非清醒状态（非4）的记录数 × 30秒 ÷ 60（转换为分钟）
        total_sleep_time = (valid_data != 4).sum() * 30 / 60
        # 2. 睡眠效率：总睡眠时间 ÷ 总记录时间 × 100%
        total_time = len(valid_data) * 30 / 60  # 总记录时间（分钟）
        sleep_efficiency = total_sleep_time / total_time * 100 if total_time > 0 else 0
        # 3. 入睡潜伏期：从开始到首次进入睡眠状态的时间（分钟）
        sleep_indices = valid_data[valid_data != 4].index
        sleep_latency = sleep_indices[0] * 30 / 60 if len(sleep_indices) > 0 else 0
        # 4. 深睡眠比例：N3阶段（3）占总睡眠时间的百分比
        deep_sleep_count = (valid_data == 3).sum()
        deep_sleep_percentage = deep_sleep_count / (total_sleep_time * 2) * 100 if total_sleep_time > 0 else 0
        # 5. 快速眼动睡眠比例：REM阶段（2）占总睡眠时间的百分比
        rem_sleep_count = (valid_data == 2).sum()
        rem_sleep_percentage = rem_sleep_count / (total_sleep_time * 2) * 100 if total_sleep_time > 0 else 0
        # 6. 睡眠开始后清醒的总次数：从睡眠到清醒的转变次数
        transitions = ((valid_data.shift(1) != 4) & (valid_data == 4)).sum()
        return pd.Series(
            [total_sleep_time, sleep_efficiency, sleep_latency,
             deep_sleep_percentage, rem_sleep_percentage, transitions],
            index=['总睡眠时间（分钟）', '睡眠效率（%）', '入睡潜伏期（分钟）',
                   '深睡眠比例（%）', '快速眼动睡眠比例（%）', '清醒总次数']
        )
    # 应用函数计算所有统计量
    results = df.apply(calculate_statistics).T.reset_index()
    results.rename(columns={'index': '原始列名'}, inplace=True)
    # 从原始列名中提取被试ID和夜晚（光照属性）
    # 使用正则表达式匹配'被试X_Night Y'格式
    pattern = r'被试(\d+)_Night (\d+)'
    results[['被试ID', '夜晚（光照属性）']] = results['原始列名'].str.extract(pattern)
    # 将提取的ID转换为整数类型
    results['被试ID'] = results['被试ID'].astype(int)
    results['夜晚（光照属性）'] = results['夜晚（光照属性）'].astype(int)
    # 调整列顺序，将被试ID和夜晚属性放在前面
    columns_order = ['被试ID', '夜晚（光照属性）', '总睡眠时间（分钟）',
                     '睡眠效率（%）', '入睡潜伏期（分钟）', '深睡眠比例（%）',
                     '快速眼动睡眠比例（%）', '清醒总次数']
    results = results[columns_order]
    # 保存结果到Excel文件
    results.to_excel(output_file_path, index=False)
    return results
# 主程序执行
if __name__ == '__main__':
    # 输入文件路径（请替换为你的实际文件路径）
    input_path = 'C:\\Users\\Administrator\\PyCharmMiscProject\\Problem 4.xlsx'
    # 输出文件路径（结果将保存到这里）
    output_path = '睡眠数据统计结果.xlsx'
    # 执行分析并保存结果
    all_results = analyze_and_save_sleep_data(input_path, output_path)
    # 显示所有结果
    print('所有被试的睡眠统计量结果：')
    print(all_results.to_string())
    print(f'\n结果已成功保存到 {output_path}')

步骤2: 统计学分析

import pandas as pd
import numpy as np
import scipy.stats as stats
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题
def analyze_light_effects(result_file_path):
    '''
    分析三种光照环境对各项睡眠指标的影响是否存在显著性差异
    输出所有数据而非仅前几行
    参数:
    result_file_path: 之前生成的睡眠数据统计结果文件路径
    '''
    # 读取之前生成的分析结果
    df = pd.read_excel(result_file_path)
    # 将中文列名替换为英文，彻底避免图表中出现中文
    df = df.rename(columns={
        '夜晚（光照属性）': 'Light_Condition',
        '被试ID': 'Subject_ID',
        '总睡眠时间（分钟）': 'Total_Sleep_Time',
        '睡眠效率（%）': 'Sleep_Efficiency',
        '入睡潜伏期（分钟）': 'Sleep_Onset_Latency',
        '深睡眠比例（%）': 'Deep_Sleep_Ratio',
        '快速眼动睡眠比例（%）': 'REM_Sleep_Ratio',
        '清醒总次数': 'Total_Awakenings'
    })
    # 查看数据基本信息
    print('数据基本信息：')
    print(f'样本量：{len(df)} 条记录')
    print(f'被试数量：{df['Subject_ID'].nunique()} 人')
    print(f'光照环境类型：{df['Light_Condition'].unique()}')
    print('\n所有数据：')
    print(df.to_string())  # 显示所有数据
    # 按光照环境分组，计算各睡眠指标的描述性统计量
    light_group_stats = df.groupby('Light_Condition').agg([
        'count', 'mean', 'std', 'min', 'max'
    ])
    print('\n各光照环境下睡眠指标的描述性统计：')
    print(light_group_stats.to_string())  # 显示所有统计结果
    # 定义要分析的睡眠指标（仅使用英文）
    sleep_metrics = [
        ('Total_Sleep_Time', 'Total Sleep Time (min)'),
        ('Sleep_Efficiency', 'Sleep Efficiency (%)'),
        ('Sleep_Onset_Latency', 'Sleep Onset Latency (min)'),
        ('Deep_Sleep_Ratio', 'Deep Sleep Ratio (%)'),
        ('REM_Sleep_Ratio', 'REM Sleep Ratio (%)'),
        ('Total_Awakenings', 'Total Awakenings')
    ]
    # 存储ANOVA结果
    anova_results = {}
    # 对每个指标进行单因素方差分析
    print('\n单因素方差分析结果（光照环境对睡眠指标的影响）：')
    for col_name, english_name in sleep_metrics:
        # 准备数据
        groups = [df[df['Light_Condition'] == i][col_name] for i in df['Light_Condition'].unique()]
        # 进行方差齐性检验
        levene_stat, levene_p = stats.levene(*groups)
        # 进行单因素方差分析
        f_stat, p_value = stats.f_oneway(*groups)
        # 存储结果
        anova_results[col_name] = {
            'F值': f_stat,
            'P值': p_value,
            '方差齐性检验P值': levene_p,
            '显著性': '显著' if p_value < 0.05 else '不显著'
        }
        # 打印结果
        print(f'\n{english_name}:')
        print(f'  方差齐性检验: Levene统计量 = {levene_stat:.4f}, P值 = {levene_p:.4f}')
        print(f'  ANOVA分析: F值 = {f_stat:.4f}, P值 = {p_value:.4f}')
        print(f'  结论: 三种光照环境对{english_name}的影响{anova_results[col_name]['显著性']} (p = {p_value:.4f})')
    # 将ANOVA结果转换为DataFrame
    anova_df = pd.DataFrame(anova_results).T
    # 可视化各光照环境下的睡眠指标差异（全英文）
    plt.figure(figsize=(18, 12))
    for i, (col_name, english_name) in enumerate(sleep_metrics, 1):
        plt.subplot(2, 3, i)
        sns.boxplot(x='Light_Condition', y=col_name, data=df)
        sns.swarmplot(x='Light_Condition', y=col_name, data=df, color='black', alpha=0.5)
        plt.title(f'Distribution of {english_name} by Light Conditions')
        plt.xlabel('Light Conditions')
        plt.ylabel(english_name)
    plt.tight_layout()
    plt.savefig('sleep_metrics_by_light_conditions.png', dpi=300)
    plt.show()
    # 修复：正确筛选存在显著差异的指标
    # 错误原因：之前的迭代方式错误，现在改为正确的元组解包方式
    significant_metrics = [
        (col_name, english_name)
        for col_name, english_name in sleep_metrics
        if anova_results[col_name]['P值'] < 0.05
    ]
    if significant_metrics:
        print('\n事后检验结果（Tukey HSD）：')
        for col_name, english_name in significant_metrics:
            print(f'\n{english_name}:')
            # 为事后检验准备数据
            data = df[['Light_Condition', col_name]].copy()
            data.columns = ['light', 'value']
            # 拟合线性模型
            model = ols('value ~ C(light)', data=data).fit()
            anova_table = anova_lm(model, typ=2)
            # 进行Tukey HSD检验
            from statsmodels.stats.multicomp import pairwise_tukeyhsd
            tukey = pairwise_tukeyhsd(endog=data['value'], groups=data['light'], alpha=0.05)
            print(tukey)
    # 恢复原始列名以便保存结果
    anova_df = anova_df.rename(index={
        'Total_Sleep_Time': '总睡眠时间（分钟）',
        'Sleep_Efficiency': '睡眠效率（%）',
        'Sleep_Onset_Latency': '入睡潜伏期（分钟）',
        'Deep_Sleep_Ratio': '深睡眠比例（%）',
        'REM_Sleep_Ratio': '快速眼动睡眠比例（%）',
        'Total_Awakenings': '清醒总次数'
    })
    return anova_df
# 主程序执行
if __name__ == '__main__':
    # 结果文件路径（与之前生成的输出文件路径一致）
    result_file = '睡眠数据统计结果.xlsx'
    # 执行光照环境影响分析
    anova_results = analyze_light_effects(result_file)
    # 保存ANOVA结果
    anova_results.to_excel('光照环境对睡眠指标影响的ANOVA结果.xlsx')
    print('\nANOVA分析结果已保存到 '光照环境对睡眠指标影响的ANOVA结果.xlsx'')
    print('\n不同光照环境下的睡眠指标分布已保存为 'sleep_metrics_by_light_conditions.png'')

代码演示与讲解视频

整体思路视频: https://www.bilibili.com/video/BV1DktBzLEVa/

第三四部分详细细节: hhttps://www.bilibili.com/video/BV1cVt1zJEeg/

完整助攻文章展示

代码文件展示

代码结果展示