预期表型值（EPD）在动物育种中的核心应用与计算详解

你好！很高兴和你聊聊动物育种中这个非常重要的概念。作为一个长期关注畜牧业遗传改良的”老伙计”，我发现很多人对EPD（Expected Progeny Difference）既熟悉又陌生——知道它很重要，但真正理解其计算原理和应用方法的人并不多。今天，我们就来把这件事讲透彻，让无论是学生还是从业者都能看懂、会用。

一、什么是预期表型值（EPD）？

1.1 一个通俗易懂的定义

想象一下，你有一头特别优秀的种猪，它的后代普遍长得又快又好。那么问题来了：这种”优秀”是因为它本身的遗传基因好，还是因为饲养管理得当呢？

EPD就是用来回答这个问题的科学工具。它不是一个动物自身表现的简单记录，而是预测这头动物传给后代的遗传能力。说得再直白一点：EPD告诉你，如果你用这头公猪配种，它的后代平均会比群体平均水平好多少（或多差多少）。

1.2 为什么EPD如此重要？

在传统育种中，人们往往根据动物自身的外貌或生产记录来选种。但这样做有几个明显的问题：

• 自身表现好，不代表遗传给后代的基因也优秀
• 环境因素（饲料、温度、管理）会干扰判断
• 不同动物在不同时间、不同条件下生产，难以直接比较

EPD的出现，正是为了解决这些痛点。它通过统计学方法，把”遗传”和”环境”分开，让你能够：

• 公平地比较不同年龄、不同饲养条件下的动物
• 预测后代的潜在表现，而不是只看父母的表现
• 在多个性状之间找到最优平衡点

二、EPD的计算原理：从数据到决策

2.1 数据的标准化处理

在计算EPD之前，所有原始数据都必须经过严格的标准化处理。这一步非常关键，因为如果不消除环境因素的影响，EPD的结果就会失真。

举个实际例子：一头母猪在冬天生产的仔猪断奶重可能比夏天生产的低1-2公斤，这主要是温度差异造成的，而不是遗传差异。在计算EPD时，需要通过统计模型把这种”季节效应”剔除掉。

标准化处理通常包括以下几个步骤：

# 假设我们有以下原始性能记录数据
import pandas as pd
import numpy as np

# 原始数据示例
data = {
    '动物ID': ['P001', 'P002', 'P003', 'P004', 'P005'],
    '性状': ['断奶重', '背膘厚', '日增重', '产仔数', '断奶重'],
    '表型值': [12.5, 14.2, 850, 12, 13.1],  # 原始记录
    '出生日期': ['2023-01-15', '2023-03-20', '2023-06-10', '2023-02-28', '2023-04-05'],
    '出生季节': ['冬季', '春季', '夏季', '冬季', '春季'],
    '性别': ['公', '公', '公', '母', '母'],
    '窝重': [10, 12, 11, 8, 10]
}

df = pd.DataFrame(data)
print("原始数据：")
print(df)

# 进行环境因子校正
# 1. 季节效应校正
season_effects = {'冬季': 0.5, '春季': 0.0, '夏季': -0.3, '秋季': 0.2}

# 2. 性别效应校正
sex_effects = {'公': 0.0, '母': -0.2}

# 3. 窝重效应（窝越大，个体断奶重可能越低）
litter_effect = -0.05  # 每多一头仔猪，个体断奶重减少0.05kg

def adjust_phenotype(row):
    """对环境因素进行校正"""
    adjusted = row['表型值']
    adjusted -= season_effects.get(row['出生季节'], 0)
    adjusted -= sex_effects.get(row['性别'], 0)
    adjusted -= litter_effect * row['窝重']
    return adjusted

df['校正后表型值'] = df.apply(adjust_phytype, axis=1)
print("\n校正后数据：")
print(df[['动物ID', '性状', '表型值', '校正后表型值']])

2.2 遗传参数的估计

校正后的数据只是第一步，接下来需要估计遗传参数。其中最核心的两个参数是：

遗传力（Heritability, h²）： 表示性状受遗传控制的程度。遗传力越高，选种效果越好。

• 高遗传力性状（h² > 0.3）：背膘厚、日增重
• 中遗传力性状（h² = 0.1-0.3）：断奶重、窝重
• 低遗传力性状（h² < 0.1）：繁殖性能、抗病力

遗传相关（Genetic Correlation, r_g）： 表示不同性状之间的遗传联系。

# 遗传参数估计示例
heredity_params = {
    '断奶重': {'遗传力': 0.25, '表型标准差': 1.8},
    '背膘厚': {'遗传力': 0.45, '表型标准差': 2.1},
    '日增重': {'遗传力': 0.35, '表型标准差': 50},
    '产仔数': {'遗传力': 0.10, '表型标准差': 2.5},
    '饲料转化率': {'遗传力': 0.30, '表型标准差': 0.15}
}

# 遗传相关矩阵示例（断奶重与背膘厚负相关）
genetic_correlation = {
    '断奶重-背膘厚': -0.35,  # 负相关：长得快的猪通常背膘厚
    '断奶重-日增重': 0.75,   # 强正相关
    '背膘厚-饲料转化率': 0.40, # 正相关：背膘厚的猪饲料转化差
    '产仔数-断奶重': -0.20    # 弱负相关：产仔多的母猪，仔猪平均断奶重略低
}

print("各性状遗传力：")
for trait, params in heredity_params.items():
    print(f"  {trait}: 遗传力 = {params['遗传力']:.2f}")

print("\n主要遗传相关：")
for pair, corr in genetic_correlation.items():
    direction = "正相关" if corr > 0 else "负相关"
    print(f"  {pair}: {corr:.2f} ({direction})")

2.3 BLUP方法：现代EPD计算的核心

BLUP（Best Linear Unbiased Prediction，最佳线性无偏预测）是目前最主流的EPD计算方法。它的核心思想是通过求解混合模型方程组，同时利用个体自身的记录、亲属记录以及系谱信息，来获得最准确的育种值估计。

混合模型的一般形式为：

\[y = Xb + Za + e\]

其中：

• y：表型值向量
• b：固定效应向量（如季节、性别、批别等）
• a：随机遗传效应向量（育种值）
• e：随机残差向量
• X、Z：设计矩阵

# 简化的BLUP计算示例
import numpy as np
from scipy.linalg import solve

def calculate_BLUP_APU(Y, A, H2, sigma_e2, sigma_a2):
    """
    简化版的BLUP计算
    Y: 表型值向量
    A: 关系矩阵（加性遗传关系矩阵）
    H2: 遗传力
    sigma_e2: 残差方差
    sigma_a2: 遗传方差
    """
    n = len(Y)
    
    # 构建K矩阵（A的加权形式）
    K = sigma_a2 * A
    
    # 构建R矩阵（对角矩阵，残差方差）
    R = sigma_e2 * np.eye(n)
    
    # 求解混合模型方程组
    # 简化版：只考虑个体自身记录
    # EBV = (A * inv(A + R/sigma_a2)) * (Y - mean(Y))
    
    # 计算系数矩阵
    C = np.linalg.inv(A + R/sigma_a2)
    
    # 计算育种值估计
    EBV = A @ C @ (Y - np.mean(Y))
    
    # EPD = 0.5 * EBV
    EPD = 0.5 * EBV
    
    return EBV, EPD

# 示例数据
# 假设我们有5头猪的断奶重记录
phenotypes = np.array([12.5, 13.2, 11.8, 14.0, 12.8])  # 断奶重(kg)

# 简化的关系矩阵（假设P001和P002是兄弟）
A_matrix = np.array([
    [1.0, 0.5, 0.0, 0.25, 0.0],
    [0.5, 1.0, 0.0, 0.25, 0.0],
    [0.0, 0.0, 1.0, 0.0, 0.5],
    [0.25, 0.25, 0.0, 1.0, 0.0],
    [0.0, 0.0, 0.5, 0.0, 1.0]
])

# 参数设定
h2 = 0.25  # 断奶重的遗传力
phenotypic_std = 1.8  # 表型标准差
sigma_a2 = h2 * phenotypic_std**2  # 遗传方差
sigma_e2 = (1 - h2) * phenotypic_std**2  # 残差方差

# 计算BLUP
EBV, EPD = calculate_BLUP_APU(phenotypes, A_matrix, h2, sigma_e2, sigma_a2)

print("各动物的育种值（EBV）和预期表型值（EPD）：")
for i, (pb, ebv, epd) in enumerate(zip(['P001', 'P002', 'P003', 'P004', 'P005'], EBV, EPD)):
    print(f"  {pb}: EBV = {ebv:.3f}, EPD = {epd:.3f}")

print(f"\n群体平均断奶重: {np.mean(phenotypes):.2f} kg")
print("EPD > 0 表示遗传潜力高于群体平均")
print("EPD < 0 表示遗传潜力低于群体平均")

三、EPD的实际计算公式

3.1 基本计算公式

EPD的基本计算公式非常简洁：

EPD = 0.5 × EBV

其中EBV是育种值（Estimated Breeding Value）。系数0.5的原因是：每个亲本只传递一半的基因给后代。

但EBV的计算要复杂得多。在实际应用中，EBV的估计需要考虑多种信息源：

def calculate_EBV_simple(self_record, relative_records, h2, reliability):
    """
    简化的EBV计算
    self_record: 个体自身记录
    relative_records: 亲属记录
    h2: 遗传力
    reliability: 准确度（0-1之间）
    """
    # 基于个体自身记录的EBV
    ebv_self = h2 * (self_record - population_mean)
    
    # 基于亲属记录的EBV（简化版）
    # 假设只有一个亲属记录
    r = 0.5  # 亲子相关系数
    ebv_relative = r * h2 * (relative_records - population_mean) / h2
    
    # 综合EBV
    ebv_combined = (reliability * ebv_self + (1 - reliability) * ebv_relative)
    
    return ebv_combined

# 实际应用示例
population_mean = 12.0  # 群体平均断奶重

# 假设P001的记录
self_record = 13.5  # P001的断奶重
relative_record = 12.8  # P001父亲的断奶重
reliability = 0.75  # 基于10条记录的准确度

ebv = calculate_EBV_simple(self_record, relative_record, 0.25, reliability)
epd = 0.5 * ebv

print(f"P001的EBV = {ebv:.3f}")
print(f"P001的EPD = {epd:.3f}")
print(f"预测：P001的后代平均断奶重 = {population_mean + epd:.2f} kg")

3.2 多性状EPD计算

在实际育种中，我们通常同时关注多个性状。这时候就需要考虑性状的遗传相关，使用多性状模型：

\[ \begin{bmatrix} y_1 \\ y_2 \\ ... \\ y_n \end{bmatrix} = \begin{bmatrix} X_1 & 0 & ... & 0 \\ 0 & X_2 & ... & 0 \\ ... & ... & ... & ... \\ 0 & 0 & ... & X_n \end{bmatrix} \begin{bmatrix} b_1 \\ b_2 \\ ... \\ b_n \end{bmatrix} + \begin{bmatrix} Z_1 & 0 & ... & 0 \\ 0 & Z_2 & ... & 0 \\ ... & ... & ... & ... \\ 0 & 0 & ... & Z_n \end{bmatrix} \begin{bmatrix} a_1 \\ a_2 \\ ... \\ a_n \end{bmatrix} + \begin{bmatrix} e_1 \\ e_2 \\ ... \\ e_n \end{bmatrix} \]

# 多性状BLUP计算示例
class MultiTrait_BLUP:
    def __init__(self, G_matrix, R_matrix, A_matrix, Y, X, Z):
        """
        G_matrix: 加性遗传协方差矩阵
        R_matrix: 残差协方差矩阵
        A_matrix: 关系矩阵
        Y: 表型值向量
        X, Z: 设计矩阵
        """
        self.G = G_matrix
        self.R = R_matrix
        self.A = A_matrix
        self.Y = Y
        self.X = X
        self.Z = Z
        
    def solve(self):
        """求解混合模型方程组"""
        n_traits = self.G.shape[0] // 2
        n_animals = self.A.shape[0]
        
        # 构建K矩阵
        K = np.kron(self.A, self.G)
        
        # 构建R矩阵
        R = np.kron(np.eye(n_animals), self.R)
        
        # 求解
        C = np.linalg.inv(K + R)
        
        # 计算育种值
        a_hat = K @ C @ self.Y
        
        # 计算EPD
        EPD = 0.5 * a_hat
        
        return EPD

# 示例：断奶重和背膘厚的双性状计算
# 遗传协方差矩阵
G = np.array([
    [1.8**2 * 0.25, -0.35 * 1.8 * 2.1 * np.sqrt(0.25 * 0.45)],
    [-0.35 * 1.8 * 2.1 * np.sqrt(0.25 * 0.45), 2.1**2 * 0.45]
])

# 残差协方差矩阵
R = np.array([
    [1.8**2 * 0.75, 0],
    [0, 2.1**2 * 0.55]
])

print("遗传协方差矩阵（断奶重-背膘厚）：")
print(G)
print("\n残差协方差矩阵：")
print(R)

3.3 不同物种的EPD应用

虽然基本原理相同，但不同物种的EPD计算有其特殊性：

猪的EPD计算特点：

• 主要性状：产仔数、断奶重、背膘厚、日增重、饲料转化率
• 数据特点：每胎产仔数多，个体记录容易获取
• 特殊考虑：窝效应需要特别处理

# 生猪EPD计算完整示例
class Pig_EPD_Calculator:
    def __init__(self):
        # 遗传参数（生猪典型值）
        self.h2 = {
            '断奶重': 0.25,
            '背膘厚': 0.45,
            '日增重': 0.35,
            '产仔数': 0.10,
            '饲料转化率': 0.30
        }
        
        self.phenotypic_std = {
            '断奶重': 1.8,
            '背膘厚': 2.1,
            '日增重': 50,
            '产仔数': 2.5,
            '饲料转化率': 0.15
        }
        
        # 遗传相关矩阵
        self.genetic_corr = {
            '断奶重-背膘厚': -0.35,
            '断奶重-日增重': 0.75,
            '背膘厚-饲料转化率': 0.40,
            '产仔数-断奶重': -0.20
        }
    
    def calculate_epd(self, animal_id, records, pedigree):
        """
        计算单头猪的EPD
        records: 性能记录字典
        pedigree: 系谱信息
        """
        epd_results = {}
        
        for trait in self.h2.keys():
            if trait in records:
                # 计算该性状的EPD
                phenotypic_value = records[trait]
                h2 = self.h2[trait]
                std = self.phenotypic_std[trait]
                
                # 简化计算：基于个体记录和系谱信息
                # 实际应用中需要使用BLUP方法
                pb_mean = 12.0  # 群体平均值（断奶重）
                
                # 计算相对选择指数
                if trait == '断奶重':
                    # 考虑多个相关性状
                    bbepd = epd_results.get('背膘厚', 0)
                    epd_value = h2 * (phenotypic_value - pb_mean)
                    # 调整背膘厚的影响
                    epd_value += self.genetic_corr.get('断奶重-背膘厚', 0) * bbepd
                else:
                    epd_value = h2 * (phenotypic_value - pb_mean)
                
                epd_results[trait] = epd_value
        
        return epd_results
    
    def calculate_selection_index(self, epd_values, weights):
        """
        计算综合选择指数
        epd_values: 各性状的EPD值
        weights: 各性状的经济权重
        """
        index_value = 0
        for trait, epd in epd_values.items():
            if trait in weights:
                index_value += weights[trait] * epd
        
        return index_value

# 应用示例
calculator = Pig_EPD_Calculator()

# 假设P001的性能记录
records_p001 = {
    '断奶重': 13.5,
    '背膘厚': 13.0,
    '日增重': 850,
    '产仔数': 12,
    '饲料转化率': 2.8
}

# 计算EPD
epd_p001 = calculator.calculate_epd('P001', records_p001, None)

print("P001各性状的EPD值：")
for trait, epd in epd_p001.items():
    print(f"  {trait}: EPD = {epd:.3f}")

# 计算选择指数（假设权重）
weights = {
    '断奶重': 1.0,
    '背膘厚': -0.5,  # 负权重，因为希望背膘薄
    '日增重': 0.8,
    '产仔数': 0.3,
    '饲料转化率': -0.6  # 负权重，希望饲料转化好
}

index = calculator.calculate_selection_index(epd_p001, weights)
print(f"\nP001的综合选择指数 = {index:.3f}")

四、EPD在实际育种中的应用

4.1 种猪选择决策

在实际生产中，EPD最直接的应用就是辅助种猪选择。一个好的种猪应该是：

• 关键性状EPD值高
• 各性状之间协调性好
• 遗传疾病风险低

def evaluate_boar_selection(boar_data, target_traits, weights):
    """
    评估公猪的选种价值
    """
    print("=" * 50)
    print(f"公猪ID: {boar_data['id']}")
    print("=" * 50)
    
    # 显示各性状EPD
    print("\n各性状EPD值：")
    for trait in target_traits:
        epd = boar_data['epd'].get(trait, 0)
        std = boar_data['std'].get(trait, 1)
        reliability = boar_data['reliability'].get(trait, 0.5)
        
        # 计算EPD的置信区间
        ci_lower = epd - 1.96 * std * np.sqrt(1 - reliability)
        ci_upper = epd + 1.96 * std * np.sqrt(1 - reliability)
        
        print(f"  {trait}: EPD = {epd:+.3f} (95% CI: [{ci_lower:+.3f}, {ci_upper:+.3f}])")
    
    # 计算综合指数
    total_index = 0
    for trait, weight in weights.items():
        epd = boar_data['epd'].get(trait, 0)
        total_index += weight * epd
    
    print(f"\n综合选择指数: {total_index:+.3f}")
    
    # 风险评估
    risk_traits = []
    for trait in target_traits:
        reliability = boar_data['reliability'].get(trait, 0.5)
        if reliability < 0.6:
            risk_traits.append(trait)
    
    if risk_traits:
        print(f"\n⚠️ 警告：以下性状准确度较低，需谨慎参考：")
        for trait in risk_traits:
            print(f"   - {trait} (准确度: {boar_data['reliability'][trait]:.1%})")
    
    # 给出选择建议
    print("\n📋 选择建议：")
    if total_index > 1.5:
        print("   ✅ 推荐：综合遗传潜力优秀，可作为核心群选种")
    elif total_index > 0.5:
        print("   ⚠️ 可用：遗传潜力良好，建议配合后裔测定验证")
    else:
        print("   ❌ 不建议：遗传潜力一般，建议淘汰")
    
    return total_index

# 实际应用示例
boar_data = {
    'id': 'XY-2024-001',
    'epd': {
        '断奶重': 0.45,
        '背膘厚': -0.32,
        '日增重': 12.5,
        '产仔数': 0.8,
        '饲料转化率': -0.08
    },
    'std': {
        '断奶重': 0.18,
        '背膘厚': 0.21,
        '日增重': 5.0,
        '产仔数': 0.25,
        '饲料转化率': 0.03
    },
    'reliability': {
        '断奶重': 0.85,
        '背膘厚': 0.78,
        '日增重': 0.92,
        '产仔数': 0.45,  # 准确度较低
        '饲料转化率': 0.88
    }
}

target_traits = ['断奶重', '背膘厚', '日增重', '产仔数', '饲料转化率']
weights = {
    '断奶重': 1.0,
    '背膘厚': -0.4,
    '日增重': 0.6,
    '产仔数': 0.3,
    '饲料转化率': -0.8
}

evaluate_boar_selection(boar_data, target_traits, weights)

4.2 杂交配套体系设计

EPD不仅在纯种选育中有用，在杂交配套体系的设计中同样关键。一个好的杂交组合应该：

• 父母本各有优势性状
• 杂种优势得到充分利用
• 后代生产性能稳定

class Crossbreeding_Strategy:
    def __init__(self):
        self.breeds = ['长白', '大白', '杜洛克']
        self.mating_systems = [
            '长白 × 大白 → 杜洛克杂交',
            '大白 × 长白 → 杜洛克杂交',
            '杜洛克 × 长白 → 大白杂交'
        ]
    
    def calculate_hybrid_vigor(self, sire_epd, dam_epd, target_trait):
        """
        计算杂种优势
        """
        # 杂种优势 = 杂交后代表现 - 纯种平均表现
        purebred_mean = (sire_epd + dam_epd) / 2
        hybrid_expected = purebred_mean * 1.15  # 假设15%的杂种优势
        
        return hybrid_expected - purebred_mean
    
    def evaluate_crossbreeding(self, sire_line, dam_line, traits):
        """
        评估杂交组合
        """
        results = {}
        
        for trait in traits:
            sire_epd = sire_line['epd'].get(trait, 0)
            dam_epd = dam_line['epd'].get(trait, 0)
            
            # 计算杂种优势
            vigor = self.calculate_hybrid_vigor(sire_epd, dam_epd, trait)
            
            # 计算杂交后代预期表现
            expected_performance = (sire_epd + dam_epd) / 2 + vigor
            
            results[trait] = {
                '父本EPD': sire_epd,
                '母本EPD': dam_epd,
                '杂种优势': vigor,
                '预期表现': expected_performance
            }
        
        return results
    
    def recommend_crossbreeding(self, sire_lines, dam_lines, traits, weights):
        """
        推荐最佳杂交组合
        """
        best_combination = None
        best_score = -float('inf')
        
        for sire in sire_lines:
            for dam in dam_lines:
                results = self.evaluate_crossbreeding(sire, dam, traits)
                
                # 计算综合评分
                total_score = 0
                for trait in traits:
                    weight = weights.get(trait, 0)
                    performance = results[trait]['预期表现']
                    total_score += weight * performance
                
                if total_score > best_score:
                    best_score = total_score
                    best_combination = {
                        '父本': sire['id'],
                        '母本': dam['id'],
                        '评分': total_score,
                        '各性状预期表现': results
                    }
        
        return best_combination

# 实际应用
crossbreeder = Crossbreeding_Strategy()

# 假设两个品系的EPD数据
sire_line_a = {
    'id': '品系A-杜洛克',
    'epd': {
        '日增重': 18.5,
        '饲料转化率': -0.12,
        '背膘厚': -0.45
    }
}

sire_line_b = {
    'id': '品系B-杜洛克',
    'epd': {
        '日增重': 15.2,
        '饲料转化率': -0.08,
        '背膘厚': -0.30
    }
}

dam_line = {
    'id': '母本品系-长白×大白',
    'epd': {
        '产仔数': 1.5,
        '断奶重': 0.35,
        '泌乳力': 0.8
    }
}

traits = ['日增重', '饲料转化率', '背膘厚', '产仔数', '断奶重', '泌乳力']
weights = {
    '日增重': 0.3,
    '饲料转化率': 0.25,
    '背膘厚': 0.15,
    '产仔数': 0.15,
    '断奶重': 0.1,
    '泌乳力': 0.05
}

best = crossbreeder.recommend_crossbreeding(
    [sire_line_a, sire_line_b], 
    [dam_line], 
    traits, 
    weights
)

print("最佳杂交组合推荐：")
print(f"父本: {best['父本']}")
print(f"母本: {best['母本']}")
print(f"综合评分: {best['评分']:.3f}")

print("\n各性状预期表现：")
for trait, data in best['各性状预期表现'].items():
    print(f"  {trait}: 预期表现 = {data['预期表现']:.3f}")

4.3 遗传进展评估

通过EPD的长期追踪，可以评估育种工作的效果。一个成功的育种计划应该实现：

• 目标性状EPD值持续提高
• 遗传进展速度符合预期
• 遗传多样性保持在合理水平

class Genetic_Progress_Tracker:
    def __init__(self, base_year, current_year):
        self.base_year = base_year
        self.current_year = current_year
        self.trait_history = {}
    
    def record_yearly_epd(self, year, traits_epd):
        """记录每年的EPD数据"""
        if year not in self.trait_history:
            self.trait_history[year] = {}
        
        for trait, epd_value in traits_epd.items():
            if trait not in self.trait_history[year]:
                self.trait_history[year][trait] = []
            self.trait_history[year][trait].append(epd_value)
    
    def calculate_genetic_progress(self, trait):
        """计算某性状的遗传进展"""
        if trait not in self.trait_history:
            return None
        
        progress_data = []
        for year in sorted(self.trait_history.keys()):
            avg_epd = np.mean(self.trait_history[year][trait])
            progress_data.append({
                'year': year,
                'avg_epd': avg_epd
            })
        
        if len(progress_data) < 2:
            return None
        
        # 计算遗传进展速率
        first_year = progress_data[0]
        last_year = progress_data[-1]
        
        years_passed = last_year['year'] - first_year['year']
        epd_change = last_year['avg_epd'] - first_year['avg_epd']
        
        annual_progress = epd_change / years_passed if years_passed > 0 else 0
        
        return {
            'trait': trait,
            'start_year': first_year['year'],
            'end_year': last_year['year'],
            'start_epd': first_year['avg_epd'],
            'end_epd': last_year['avg_epd'],
            'total_change': epd_change,
            'annual_progress': annual_progress,
            'progress_per_year_pct': (annual_progress / first_year['avg_epd'] * 100) if first_year['avg_epd'] != 0 else 0
        }
    
    def generate_progress_report(self):
        """生成遗传进展报告"""
        report = {
            'base_year': self.base_year,
            'current_year': self.current_year,
            'traits_progress': {}
        }
        
        for trait in self.trait_history:
            progress = self.calculate_genetic_progress(trait)
            if progress:
                report['traits_progress'][trait] = progress
        
        return report
    
    def visualize_progress(self):
        """可视化遗传进展"""
        import matplotlib.pyplot as plt
        import matplotlib
        matplotlib.use('Agg')
        
        fig, axes = plt.subplots(2, 2, figsize=(12, 10))
        axes = axes.flatten()
        
        for idx, (trait, progress) in enumerate(self.trait_history.items()):
            if idx >= 4:
                break
            
            years = [data['year'] for data in progress]
            epd_values = [data['avg_epd'] for data in progress]
            
            axes[idx].plot(years, epd_values, 'o-', linewidth=2, markersize=8)
            axes[idx].set_title(f'{trait} EPD趋势', fontsize=12, fontweight='bold')
            axes[idx].set_xlabel('年份')
            axes[idx].set_ylabel('EPD值')
            axes[idx].grid(True, alpha=0.3)
            
            # 添加趋势线
            if len(years) >= 2:
                z = np.polyfit(years, epd_values, 1)
                p = np.poly1d(z)
                axes[idx].plot(years, p(years), "r--", alpha=0.6, label='趋势线')
                axes[idx].legend()
        
        plt.tight_layout()
        plt.savefig('genetic_progress.png', dpi=150, bbox_inches='tight')
        print("图表已保存为 'genetic_progress.png'")

# 实际应用示例
tracker = Genetic_Progress_Tracker(2015, 2024)

# 模拟过去10年的EPD数据
np.random.seed(42)
for year in range(2015, 2025):
    # 模拟各性状EPD的逐年进步
    yearly_epd = {
        '断奶重': 0.15 + (year - 2015) * 0.03 + np.random.normal(0, 0.02),
        '背膘厚': -0.08 - (year - 2015) * 0.02 + np.random.normal(0, 0.01),
        '日增重': 1.5 + (year - 2015) * 0.8 + np.random.normal(0, 0.3),
        '产仔数': 0.05 + (year - 2015) * 0.02 + np.random.normal(0, 0.01)
    }
    tracker.record_yearly_epd(year, yearly_epd)

# 生成报告
report = tracker.generate_progress_report()

print("\n" + "="*60)
print("遗传进展评估报告")
print("="*60)
print(f"评估期间: {report['base_year']} - {report['current_year']}年")
print()

for trait, progress in report['traits_progress'].items():
    print(f"【{trait}】")
    print(f"  起始EPD: {progress['start_epd']:+.3f}")
    print(f"  结束EPD: {progress['end_epd']:+.3f}")
    print(f"  总变化: {progress['total_change']:+.3f}")
    print(f"  年进展: {progress['annual_progress']:+.3f}/年")
    print(f"  进展率: {progress['progress_per_year_pct']:+.1f}%/年")
    print()

# 生成可视化图表
tracker.visualize_progress()

五、EPD应用的常见问题与解决方案

5.1 数据质量问题

在实际应用中，数据质量是影响EPD准确性的最大因素。常见问题包括：

def data_quality_check(performance_data):
    """
    数据质量检查函数
    """
    quality_issues = {
        'missing_values': [],
        'outliers': [],
        'inconsistent_records': [],
        'missing_pedigree': []
    }
    
    for record in performance_data:
        # 检查缺失值
        if pd.isna(record.get('表型值')):
            quality_issues['missing_values'].append(record['动物ID'])
        
        # 检查异常值（使用IQR方法）
        if '表型值' in record:
            q1 = np.percentile([r['表型值'] for r in performance_data if not pd.isna(r.get('表型值'))], 25)
            q3 = np.percentile([r['表型值'] for r in performance_data if not pd.isna(r.get('表型值'))], 75)
            iqr = q3 - q1
            
            if record['表型值'] < q1 - 1.5 * iqr or record['表型值'] > q3 + 1.5 * iqr:
                quality_issues['outliers'].append(record['动物ID'])
        
        # 检查系谱信息
        if not record.get('父本ID') or not record.get('母本ID'):
            quality_issues['missing_pedigree'].append(record['动物ID'])
    
    return quality_issues

def suggest_data_correction(issues):
    """根据问题类型提供修正建议"""
    suggestions = []
    
    if issues['missing_values']:
        suggestions.append(f"发现 {len(issues['missing_values'])} 条缺失表型值记录，建议：")
        suggestions.append("   1. 重新核实原始记录")
        suggestions.append("   2. 使用KNN插补或均值插补")
        suggestions.append("   3. 删除无法恢复的记录")
    
    if issues['outliers']:
        suggestions.append(f"发现 {len(issues['outliers'])} 条异常值记录，建议：")
        suggestions.append("   1. 核实是否为记录错误")
        suggestions.append("   2. 使用winsorization方法处理")
        suggestions.append("   3. 或者使用稳健统计方法")
    
    if issues['missing_pedigree']:
        suggestions.append(f"发现 {len(issues['missing_pedigree'])} 条缺失系谱记录，建议：")
        suggestions.append("   1. 补充系谱信息")
        suggestions.append("   2. 使用基因组信息替代系谱")
        suggestions.append("   3. 将这些个体单独分析")
    
    return suggestions

# 示例应用
sample_data = [
    {'动物ID': 'P001', '表型值': 12.5, '父本ID': 'S001', '母本ID': 'D001'},
    {'动物ID': 'P002', '表型值': None, '父本ID': 'S002', '母本ID': 'D002'},
    {'动物ID': 'P003', '表型值': 50.0, '父本ID': None, '母本ID': 'D003'},
    {'动物ID': 'P004', '表型值': 13.2, '父本ID': 'S004', '母本ID': None},
]

issues = data_quality_check(sample_data)
suggestions = suggest_data_correction(issues)

print("数据质量检查结果：")
for category, items in issues.items():
    if items:
        print(f"  {category}: {len(items)}条问题记录")

print("\n修正建议：")
for suggestion in suggestions:
    print(suggestion)

5.2 遗传多样性保护

在追求遗传进展的同时，不能忽视遗传多样性的保护。过度选种可能导致近交系数上升，遗传基础变窄。

class Genetic_Diversity_Monitor:
    def __init__(self, population_size, generation_interval):
        self.population_size = population_size
        self.generation_interval = generation_interval  # 代间隔（年）
        self.inbreeding_coefficients = {}
        self.effective_population_size = population_size
    
    def calculate_inbreeding(self, pedigree):
        """计算个体近交系数"""
        inbreeding = {}
        
        for animal_id in pedigree:
            parents = pedigree[animal_id]
            
            if not parents['father'] or not parents['mother']:
                inbreeding[animal_id] = 0
                continue
            
            # 简化计算：如果父母有共同祖先，则近交系数>0
            # 实际应用中需要使用更复杂的方法
            common_ancestors = self.find_common_ancestors(pedigree, parents['father'], parents['mother'])
            
            if common_ancestors:
                # 近交系数 = Σ(1/2)^(n1+n2+1) * (1+FA)
                inbreeding_coeff = 0
                for ancestor in common_ancestors:
                    n1 = self.get_generation_distance(pedigree, parents['father'], ancestor)
                    n2 = self.get_generation_distance(pedigree, parents['mother'], ancestor)
                    inbreeding_coeff += (0.5 ** (n1 + n2 + 1))
                
                inbreeding[animal_id] = inbreeding_coeff
            else:
                inbreeding[animal_id] = 0
        
        return inbreeding
    
    def find_common_ancestors(self, pedigree, animal1, animal2, max_depth=10):
        """查找共同祖先"""
        ancestors1 = self.get_all_ancestors(pedigree, animal1, max_depth)
        ancestors2 = self.get_all_ancestors(pedigree, animal2, max_depth)
        
        return list(set(ancestors1) & set(ancestors2))
    
    def get_all_ancestors(self, pedigree, animal_id, max_depth):
        """获取所有祖先"""
        ancestors = []
        
        if animal_id not in pedigree:
            return ancestors
        
        parents = pedigree[animal_id]
        
        if parents['father'] and parents['father'] not in ancestors:
            ancestors.append(parents['father'])
            ancestors.extend(self.get_all_ancestors(pedigree, parents['father'], max_depth - 1))
        
        if parents['mother'] and parents['mother'] not in ancestors:
            ancestors.append(parents['mother'])
            ancestors.extend(self.get_all_ancestors(pedigree, parents['mother'], max_depth - 1))
        
        return ancestors
    
    def get_generation_distance(self, pedigree, animal_id, ancestor_id, depth=0):
        """获取代距离"""
        if animal_id == ancestor_id:
            return depth
        
        if animal_id not in pedigree:
            return float('inf')
        
        parents = pedigree[animal_id]
        
        min_distance = float('inf')
        
        if parents['father']:
            dist = self.get_generation_distance(pedigree, parents['father'], ancestor_id, depth + 1)
            min_distance = min(min_distance, dist)
        
        if parents['mother']:
            dist = self.get_generation_distance(pedigree, parents['mother'], ancestor_id, depth + 1)
            min_distance = min(min_distance, dist)
        
        return min_distance
    
    def monitor_diversity(self, pedigree, current_generation):
        """监测遗传多样性"""
        inbreeding = self.calculate_inbreeding(pedigree)
        
        avg_inbreeding = np.mean(list(inbreeding.values()))
        max_inbreeding = np.max(list(inbreeding.values()))
        
        # 计算有效群体大小
        # Ne = 4 / (Ft + Ft-1)，其中Ft是当前近交系数
        if avg_inbreeding > 0:
            effective_size = 1 / avg_inbreeding if avg_inbreeding > 0 else self.population_size
        else:
            effective_size = self.population_size
        
        return {
            'generation': current_generation,
            'avg_inbreeding': avg_inbreeding,
            'max_inbreeding': max_inbreeding,
            'effective_population_size': effective_size,
            'recommendations': []
        }
    
    def generate_diversity_report(self, pedigree, generations_data):
        """生成遗传多样性报告"""
        report = {
            'generations': [],
            'warnings': [],
            'recommendations': []
        }
        
        for gen_data in generations_data:
            monitoring_result = self.monitor_diversity(pedigree, gen_data['generation'])
            report['generations'].append(monitoring_result)
            
            # 生成警告
            if monitoring_result['avg_inbreeding'] > 0.0625:  # 6.25%
                report['warnings'].append(
                    f"第{monitoring_result['generation']}代：平均近交系数过高（{monitoring_result['avg_inbreeding']:.2%}）"
                )
            
            if monitoring_result['effective_population_size'] < 50:
                report['warnings'].append(
                    f"第{monitoring_result['generation']}代：有效群体大小过小（{monitoring_result['effective_population_size']:.0f}）"
                )
        
        # 生成建议
        if report['warnings']:
            report['recommendations'].append("建议采取以下措施保护遗传多样性：")
            report['recommendations'].append("1. 控制选种强度，适当扩大种公猪使用数量")
            report['recommendations'].append("2. 引入外血，增加遗传多样性")
            report['recommendations'].append("3. 建立最小近交配对制度")
            report['recommendations'].append("4. 定期评估遗传进展与多样性平衡")
        
        return report

# 实际应用
diversity_monitor = Genetic_Diversity_Monitor(population_size=200, generation_interval=3)

# 模拟系谱数据
simulated_pedigree = {
    'P001': {'father': 'S001', 'mother': 'D001'},
    'P002': {'father': 'S001', 'mother': 'D002'},
    'P003': {'father': 'S002', 'mother': 'D001'},
    'P004': {'father': 'S002', 'mother': 'D003'},
    'P005': {'father': 'S003', 'mother': 'D002'},
    # ... 更多个体
}

generations_data = [
    {'generation': 1, 'avg_inbreeding': 0.02},
    {'generation': 2, 'avg_inbreeding': 0.04},
    {'generation': 3, 'avg_inbreeding': 0.06},
    {'generation': 4, 'avg_inbreeding': 0.08},
]

report = diversity_monitor.generate_diversity_report(simulated_pedigree, generations_data)

print("遗传多样性监测报告：")
print("=" * 60)

for gen in report['generations']:
    print(f"\n第{gen['generation']}代：")
    print(f"  平均近交系数: {gen['avg_inbreeding']:.2%}")
    print(f"  最大近交系数: {gen['max_inbreeding']:.2%}")
    print(f"  有效群体大小: {gen['effective_population_size']:.0f}")

if report['warnings']:
    print("\n⚠️ 警告：")
    for warning in report['warnings']:
        print(f"  {warning}")

if report['recommendations']:
    print("\n📋 建议：")
    for rec in report['recommendations']:
        print(f"  {rec}")

六、EPD的未来发展趋势

6.1 基因组选择的应用

随着基因组测序技术的普及，基因组选择正在改变EPD的计算方式。基因组选择可以在动物幼年时期就准确预测其育种值，大大缩短了世代间隔。

class Genomic_EPD_Calculator:
    def __init__(self, snp_matrix, phenotype_data, pedigree):
        """
        snp_matrix: SNP基因型矩阵 (n_individuals x n_snps)
        phenotype_data: 表型数据
        pedigree: 系谱信息
        """
        self.snp_matrix = snp_matrix
        self.phenotype_data = phenotype_data
        self.pedigree = pedigree
        self.gwas_results = None
        self.marker_effects = None
    
    def perform_gwas(self, trait, min_maf=0.05, p_threshold=1e-5):
        """
        进行全基因组关联分析（GWAS）
        """
        print(f"正在进行{trait}的全基因组关联分析...")
        
        gwas_results = []
        n_snps = self.snp_matrix.shape[1]
        
        for i in range(n_snps):
            # 简化版GWAS：线性回归
            x = self.snp_matrix[:, i].astype(float)
            y = self.phenotype_data[trait]
            
            # 计算MAF
            allele_freq = np.mean(x) / 2
            if allele_freq < min_maf or allele_freq > (1 - min_maf):
                continue
            
            # 简单线性回归
            if np.std(x) > 0 and np.std(y) > 0:
                correlation = np.corrcoef(x, y)[0, 1]
                # 简化p值计算
                t_stat = correlation * np.sqrt(len(y) - 2) / np.sqrt(1 - correlation**2)
                p_value = 2 * (1 - self._normal_cdf(abs(t_stat)))
            else:
                p_value = 1.0
            
            gwas_results.append({
                'snp_id': f'SNP_{i+1}',
                'chromosome': (i // (n_snps // 18)) + 1,  # 假设18条染色体
                'position': i % 100000000,
                'maf': allele_freq,
                'r2': correlation**2 if np.std(x) > 0 and np.std(y) > 0 else 0,
                'p_value': p_value
            })
        
        # 筛选显著SNP
        significant_snps = [r for r in gwas_results if r['p_value'] < p_threshold]
        self.gwas_results = significant_snps
        
        print(f"发现 {len(significant_snps)} 个显著关联SNP")
        
        return significant_snps
    
    def calculate_gibf(self, gwas_results, heritability):
        """
        计算基因信息育种值（GBLUP）
        """
        print("\n计算基因组育种值...")
        
        # 构建G矩阵（基因组关系矩阵）
        n_individuals = self.snp_matrix.shape[0]
        
        # 简化版G矩阵计算
        M = self.snp_matrix - 2 * np.mean(self.snp_matrix, axis=0)
        G = (M @ M.T) / np.sum(np.var(self.snp_matrix, axis=0) * 2)
        
        # GBLUP计算
        # 简化版
        H2 = heritability
        sigma_a2 = H2 * np.var(self.phenotype_data['表型值'])
        sigma_e2 = (1 - H2) * np.var(self.phenotype_data['表型值'])
        
        # 求解混合模型
        K = sigma_a2 * G
        R = sigma_e2 * np.eye(n_individuals)
        
        try:
            C = np.linalg.inv(K + R)
            EBV = K @ C @ (self.phenotype_data['表型值'] - np.mean(self.phenotype_data['表型值']))
            GEBV = EBV
        except np.linalg.LinAlgError:
            print("矩阵奇异，使用简化方法")
            GEBV = 0.5 * np.dot(G, self.phenotype_data['表型值'])
        
        self.marker_effects = GEBV
        
        return GEBV
    
    def predict_genomic_epd(self, individual_id, trait, reliability=0.8):
        """
        预测个体的基因组EPD
        """
        if self.marker_effects is None:
            print("请先进行基因组育种值计算")
            return None
        
        if individual_id not in self.phenotype_data['索引']:
            print(f"未找到个体 {individual_id}")
            return None
        
        idx = self.phenotype_data['索引'].index(individual_id)
        
        # 基因组EBV
        genomic_ebv = self.marker_effects[idx]
        
        # 考虑准确度的EPD
        epd = 0.5 * genomic_ebv * np.sqrt(reliability)
        
        return {
            'individual_id': individual_id,
            'trait': trait,
            'gebv': genomic_ebv,
            'epd': epd,
            'reliability': reliability,
            'conf_interval': [
                epd - 1.96 * abs(epd) * np.sqrt(1 - reliability),
                epd + 1.96 * abs(epd) * np.sqrt(1 - reliability)
            ]
        }
    
    def _normal_cdf(self, x):
        """标准正态分布累积函数"""
        return 0.5 * (1 + np.erf(x / np.sqrt(2)))

# 应用示例
np.random.seed(42)

# 模拟数据
n_individuals = 500
n_snps = 50000

# SNP矩阵（0, 1, 2表示基因型）
snp_matrix = np.random.randint(0, 3, (n_individuals, n_snps))

# 表型数据（假设受100个SNP位点影响）
true_effects = np.zeros(n_snps)
true_effects[:100] = np.random.normal(0, 0.5, 100)  # 100个QTL
genetic_value = snp_matrix @ true_effects
phenotype = genetic_value + np.random.normal(0, 10, n_individuals)

phenotype_data = {
    '索引': [f'P{i+1}' for i in range(n_individuals)],
    '表型值': phenotype
}

# 创建基因组EPD计算器
genomic_calculator = Genomic_EPD_Calculator(snp_matrix, phenotype_data, None)

# 执行GWAS
significant_snps = genomic_calculator.perform_gwas('表型值', min_maf=0.05, p_threshold=1e-5)

# 计算GBLUP
heritability = 0.4
gebv = genomic_calculator.calculate_gibf(significant_snps, heritability)

# 预测基因组EPD
prediction = genomic_calculator.predict_genomic_epd('P1', '表型值', reliability=0.85)

if prediction:
    print(f"\n基因组EPD预测结果：")
    print(f"  个体: {prediction['individual_id']}")
    print(f"  性状: {prediction['trait']}")
    print(f"  基因组EBV: {prediction['gebv']:.3f}")
    print(f"  基因组EPD: {prediction['epd']:.3f}")
    print(f"  准确度: {prediction['reliability']:.0%}")
    print(f"  95%置信区间: [{prediction['conf_interval'][0]:.3f}, {prediction['conf_interval'][1]:.3f}]")

6.2 人工智能与大数据应用

AI和大数据技术正在让EPD计算更加智能化：

class AI_EPD_Predictor:
    def __init__(self):
        self.model = None
        self.feature_importance = None
    
    def prepare_training_data(self, performance_data, pedigree_data, env_data):
        """准备AI训练数据"""
        # 特征工程
        features = []
        
        for record in performance_data:
            feature_vector = [
                record.get('断奶重', 0),
                record.get('背膘厚', 0),
                record.get('日增重', 0),
                record.get('饲料转化率', 0),
                record.get('产仔数', 0),
                record.get('窝重', 0),
                record.get('胎次', 1),
                self._encode_season(record.get('出生季节', '春季')),
                self._encode_sex(record.get('性别', '公'))
            ]
            
            # 添加系谱特征
            if '父本EPD' in record:
                feature_vector.extend([
                    record.get('父本EPD', {}).get('断奶重', 0),
                    record.get('父本EPD', {}).get('背膘厚', 0),
                    record.get('母本EPD', {}).get('断奶重', 0),
                    record.get('母本EPD', {}).get('背膘厚', 0)
                ])
            
            # 添加环境特征
            if '饲养条件' in record:
                feature_vector.extend([
                    record.get('饲养条件', {}).get('温度', 20),
                    record.get('饲养条件', {}).get('湿度', 60),
                    record.get('饲养条件', {}).get('饲养密度', 1.5)
                ])
            
            features.append(feature_vector)
        
        return np.array(features)
    
    def _encode_season(self, season):
        """编码季节"""
        encoding = {'春季': 0, '夏季': 1, '秋季': 2, '冬季': 3}
        return encoding.get(season, 0)
    
    def _encode_sex(self, sex):
        """编码性别"""
        encoding = {'公': 0, '母': 1}
        return encoding.get(sex, 0)
    
    def train_ml_model(self, X_train, y_train):
        """训练机器学习模型"""
        from sklearn.ensemble import RandomForestRegressor
        from sklearn.model_selection import cross_val_score
        
        # 使用随机森林模型
        model = RandomForestRegressor(
            n_estimators=100,
            max_depth=15,
            min_samples_split=5,
            random_state=42
        )
        
        # 训练模型
        model.fit(X_train, y_train)
        
        # 交叉验证
        scores = cross_val_score(model, X_train, y_train, cv=5, scoring='r2')
        
        self.model = model
        self.feature_importance = model.feature_importances_
        
        print(f"模型训练完成！")
        print(f"交叉验证R²: {scores.mean():.3f} ± {scores.std():.3f}")
        
        return scores.mean()
    
    def predict_epd(self, new_record):
        """预测新个体的EPD"""
        if self.model is None:
            print("请先训练模型")
            return None
        
        # 准备特征
        feature_vector = self.prepare_training_data([new_record], None, None)[0]
        
        # 预测
        prediction = self.model.predict([feature_vector])[0]
        
        # 获取特征重要性
        importance_dict = dict(zip(
            ['断奶重', '背膘厚', '日增重', '饲料转化率', '产仔数', '窝重', '胎次', '季节', '性别'],
            self.feature_importance[:9]
        ))
        
        return {
            'predicted_epd': prediction,
            'feature_importance': importance_dict,
            'top_factors': sorted(importance_dict.items(), key=lambda x: x[1], reverse=True)[:3]
        }
    
    def analyze_selection_strategy(self, population_data, target_traits):
        """分析最佳选种策略"""
        if self.model is None:
            print("请先训练模型")
            return None
        
        predictions = []
        
        for animal in population_data:
            pred = self.predict_epd(animal)
            if pred:
                predictions.append({
                    'animal_id': animal.get('动物ID'),
                    'predicted_epd': pred['predicted_epd'],
                    'top_factors': pred['top_factors']
                })
        
        # 按预测EPD排序
        predictions.sort(key=lambda x: x['predicted_epd'], reverse=True)
        
        # 推荐选种
        top_candidates = predictions[:10]
        
        return {
            'recommendations': top_candidates,
            'selection_strategy': self._generate_strategy(top_candidates, target_traits)
        }
    
    def _generate_strategy(self, candidates, target_traits):
        """生成选种策略建议"""
        strategy = {
            'primary_selection': [],
            'backup_selection': [],
            'genetic_diversity_considerations': []
        }
        
        # 主要选种标准
        for candidate in candidates[:3]:
            strategy['primary_selection'].append({
                'animal_id': candidate['animal_id'],
                'predicted_epd': candidate['predicted_epd'],
                'key_factors': candidate['top_factors']
            })
        
        # 备选方案
        for candidate in candidates[3:6]:
            strategy['backup_selection'].append({
                'animal_id': candidate['animal_id'],
                'predicted_epd': candidate['predicted_epd']
            })
        
        # 遗传多样性考虑
        if len(candidates) > 0:
            unique_parents = set()
            for candidate in candidates[:10]:
                # 假设从动物ID中提取父本信息
                parent_id = candidate['animal_id'][:3]
                unique_parents.add(parent_id)
            
            if len(unique_parents) < 5:
                strategy['genetic_diversity_considerations'].append(
                    "注意：候选种猪遗传背景相似度高，建议扩大选种范围"
                )
        
        return strategy

# 应用示例
ai_predictor = AI_EPD_Predictor()

# 模拟训练数据
np.random.seed(42)
n_samples = 1000

training_data = []
for i in range(n_samples):
    record = {
        '断奶重': np.random.normal(12, 1.5),
        '背膘厚': np.random.normal(14, 2),
        '日增重': np.random.normal(800, 80),
        '饲料转化率': np.random.normal(2.8, 0.2),
        '产仔数': np.random.normal(12, 2),
        '窝重': np.random.normal(18, 2),
        '胎次': np.random.randint(1, 8),
        '出生季节': np.random.choice(['春季', '夏季', '秋季', '冬季']),
        '性别': np.random.choice(['公', '母']),
        '父本EPD': {
            '断奶重': np.random.normal(0.3, 0.2),
            '背膘厚': np.random.normal(-0.2, 0.15)
        },
        '母本EPD': {
            '断奶重': np.random.normal(0.2, 0.2),
            '背膘厚': np.random.normal(-0.15, 0.15)
        }
    }
    training_data.append(record)

# 准备训练数据
X_train = ai_predictor.prepare_training_data(training_data, None, None)
y_train = np.array([r['断奶重'] for r in training_data])

# 训练模型
r2_score = ai_predictor.train_ml_model(X_train, y_train)

# 预测新个体
new_animal = {
    '动物ID': 'NEW001',
    '断奶重': 13.5,
    '背膘厚': 13.0,
    '日增重': 850,
    '饲料转化率': 2.7,
    '产仔数': 12,
    '窝重': 19,
    '胎次': 2,
    '出生季节': '春季',
    '性别': '公',
    '父本EPD': {'断奶重': 0.5, '背膘厚': -0.3},
    '母本EPD': {'断奶重': 0.4, '背膘厚': -0.2}
}

prediction = ai_predictor.predict_epd(new_animal)

if prediction:
    print(f"\nAI预测结果：")
    print(f"  预测EPD: {prediction['predicted_epd']:.3f}")
    print(f"\n  最重要影响因素：")
    for factor in prediction['top_factors']:
        print(f"    - {factor[0]}: 重要性 {factor[1]:.3f}")

# 分析选种策略
population_sample = training_data[:50]
strategy = ai_predictor.analyze_selection_strategy(population_sample, ['断奶重', '背膘厚'])

if strategy:
    print(f"\nAI选种策略建议：")
    print(f"  主要候选种猪：")
    for i, candidate in enumerate(strategy['recommendations']['primary_selection'], 1):
        print(f"    {i}. {candidate['animal_id']}: EPD = {candidate['predicted_epd']:.3f}")
        print(f"       关键因素: {candidate['key_factors']}")
    
    if strategy['recommendations']['genetic_diversity_considerations']:
        print(f"\n  遗传多样性警告：")
        for warning in strategy['recommendations']['genetic_diversity_considerations']:
            print(f"    ⚠️ {warning}")

七、结语：EPD让育种更科学

聊了这么多，我想用一个简单的比喻来总结EPD的意义：

EPD就像是一只”透视眼”——它让我们能够透过表型的外在表现，看到动物真正的遗传潜力。过去，我们只能根据猪长得好不好来选种；现在，我们可以预测这只猪的后代会有多好。

对于养殖户来说，EPD意味着更准确的选种决策，更低的试错成本，更快的遗传进展。对于育种企业来说，EPD是连接科研与生产的桥梁，让遗传改良的成果能够真正转化为经济效益。

当然，EPD也不是万能的。它依赖于高质量的数据、科学的计算方法和合理的应用策略。但无论如何，EPD的出现让动物育种从”经验主义”走向”数据驱动”，这是整个行业的巨大进步。

希望这篇文章能帮助你更好地理解EPD的原理和应用。如果你在实际工作中遇到具体问题，欢迎继续交流！记住，好的育种不仅是科学，也是一门艺术——需要数据，也需要经验，更需要对动物的热爱。