GAM部分¶
In [2]:
import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
from pygam import GAM, s, te, l
import matplotlib.pyplot as plt
import warnings
from sklearn.model_selection import KFold
warnings.filterwarnings("ignore")
# ====================================================================
# 1. 数据预处理与模型训练
# ====================================================================
# 设置中文显示
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 加载数据（确保 '附件_processed.csv' 文件存在）
try:
df = pd.read_csv('附件_processed.csv', encoding='gbk')
except FileNotFoundError:
print("Error: '附件_processed.csv' not found. Please check file path.")
exit()
# 数据清洗与筛选
num_cols = df.select_dtypes(include=[np.number]).columns.tolist()
# 需要排除的列
exclude_cols = ['身高', '体重', '年龄', '孕妇BMI']
# 实际参与筛选的数值型列
num_cols = [col for col in num_cols if col not in exclude_cols]
# 逐列筛选
for col in num_cols:
mean = df[col].mean()
std = df[col].std()
df = df[(df[col] >= mean - 3 * std) & (df[col] <= mean + 3 * std)]
df = df[df['Y染色体浓度'] <= 0.20].copy()
df = df[df['原始读段数'] >= 3000000].copy()
df = df.dropna().reset_index(drop=True)
# # 假设我们模拟一个均值为0，标准差为0.005的随机噪声
# # 您可以根据实际的仪器误差情况调整标准差
# noise_y = np.random.normal(0, 0.005, len(df))
# df['Y染色体浓度'] = df['Y染色体浓度'] + noise_y
# # 假设GC含量也有轻微误差
# noise_gc = np.random.normal(0, 0.01, len(df))
# df['GC含量'] = df['GC含量'] + noise_gc
# 连续压缩，以避免因 Y 染色体浓度为 0 或 1 而导致的模型问题
n = len(df)
epsilon = (np.log(n - 1) + 0.5) / n
df['y_it'] = df['Y染色体浓度'].clip(epsilon, 1 - epsilon)
# 孕周标准化：保存 StandardScaler 对象以供预测时复用
scaler_g = StandardScaler()
df['g_it_std'] = scaler_g.fit_transform(df[['检测孕周']])
# BMI 分位数归一化：保存用于排名的 Series
# 这将确保任何新的 BMI 输入都能根据训练集中的排名进行正确归一化
bmi_for_rank = df['孕妇BMI'].copy()
df['b_i_norm'] = bmi_for_rank.rank(pct=True)
# 交叉验证
def cross_validate_gam(df, n_splits=5):
kf = KFold(n_splits=n_splits, shuffle=True, random_state=42)
losses = []
for fold, (train_idx, test_idx) in enumerate(kf.split(df)):
df_train, df_test = df.iloc[train_idx].copy(), df.iloc[test_idx].copy()
# 标准化
scaler_g = StandardScaler()
df_train['g_it_std'] = scaler_g.fit_transform(df_train[['检测孕周']])
df_test['g_it_std'] = scaler_g.transform(df_test[['检测孕周']])
# BMI分位数
bmi_for_rank = df_train['孕妇BMI'].copy()
df_train['b_i_norm'] = bmi_for_rank.rank(pct=True)
combined_bmi = pd.concat([bmi_for_rank, df_test['孕妇BMI']])
df_test['b_i_norm'] = combined_bmi.rank(pct=True).iloc[-len(df_test):].values
# 类别变量编码
for col in ['IVF妊娠', '怀孕次数', '生产次数']:
df_train[col] = df_train[col].astype('category')
df_test[col] = df_test[col].astype('category')
X_train_smooth = df_train[['g_it_std', 'b_i_norm']]
X_train_linear = df_train[['年龄', '身高']].astype(float)
X_train_cat = pd.get_dummies(df_train[['IVF妊娠', '怀孕次数', '生产次数']], drop_first=True)
X_train = pd.concat([X_train_smooth, X_train_linear, X_train_cat], axis=1)
y_train = df_train['y_it']
# 测试集同样的列顺序
X_test_smooth = df_test[['g_it_std', 'b_i_norm']]
X_test_linear = df_test[['年龄', '身高']].astype(float)
X_test_cat = pd.get_dummies(df_test[['IVF妊娠', '怀孕次数', '生产次数']], drop_first=True)
# 保证测试集和训练集列一致
X_test_cat = X_test_cat.reindex(columns=X_train_cat.columns, fill_value=0)
X_test = pd.concat([X_test_smooth, X_test_linear, X_test_cat], axis=1)
y_test = df_test['y_it']
# 训练GAM
gam = GAM(s(0) + s(1) + te(0, 1) + l(2) + l(3) + l(4) + l(5) + l(6), distribution='gamma', link='log').fit(X_train, y_train)
y_pred = gam.predict(X_test)
mse = np.mean((y_pred - y_test) ** 2)
losses.append(mse)
print(f"Fold {fold+1}: MSE = {mse:.6f}")
print(f"平均MSE: {np.mean(losses):.6f}")
# 交叉验证评估 GAM
cross_validate_gam(df, n_splits=5)
# ...existing code...
# 假设df中 '身高' 为数值型，'IVF妊娠'、'怀孕次数'、'生产次数' 为类别型变量
# 先将类别变量转为category类型（如尚未转换）
df['IVF妊娠'] = df['IVF妊娠'].astype('category')
df['怀孕次数'] = df['怀孕次数'].astype('category')
df['生产次数'] = df['生产次数'].astype('category')
# GAM特征准备
X_smooth = df[['g_it_std', 'b_i_norm']]
X_linear = df[['年龄', '身高']].astype(float)
# 将类别变量用get_dummies编码
X_cat = pd.get_dummies(df[['IVF妊娠', '怀孕次数', '生产次数']], drop_first=True)
X = pd.concat([X_smooth, X_linear, X_cat], axis=1)
y = df['y_it']
# GAM公式：孕周、BMI为平滑项，年龄、身高、类别变量为线性项
gam = GAM(
s(0) + s(1) + te(0, 1) + l(2) + l(3) + l(4) + l(5) + l(6), # 线性项数量根据哑变量展开数自动调整
distribution='gamma',
link='log'
).fit(X, y)
print("GAM（含身高、IVF、怀孕次数、生产次数）模型已成功拟合。")
Fold 1: MSE = 0.000711 Fold 2: MSE = 0.000886 Fold 3: MSE = 0.000564 Fold 4: MSE = 0.000266 Fold 5: MSE = 0.000923 平均MSE: 0.000670 GAM（含身高、IVF、怀孕次数、生产次数）模型已成功拟合。
部分数据处理+代码版本1¶
In [ ]:
# 去掉同一个孕妇的重复记录，只保留第一次
df = df.drop_duplicates(subset=['孕妇代码'], keep='first')
import numpy as np
import pandas as pd
import random
import numpy as np
import pandas as pd
import random
# =========================
# 1. GA 参数
# =========================
SEED = 42
random.seed(SEED)
np.random.seed(SEED)
MAX_G = 6
MIN_G = 2
GROUP_PENALTY = -0.5
POP_SIZE = 200
N_GEN = 200
MUT_RATE = 0.7
CROSS_RATE = 0.7
#BMI_min, BMI_max = df['孕妇BMI'].min(), df['孕妇BMI'].max()
BMI_min, BMI_max = 26, 47
weeks_available = np.arange(14, 30)
N = df.shape[0]
# ...existing code...
# 先将类别变量转为category类型（如尚未转换）
df['IVF妊娠'] = df['IVF妊娠'].astype('category')
df['怀孕次数'] = df['怀孕次数'].astype('category')
df['生产次数'] = df['生产次数'].astype('category')
# 特征工程（与交叉验证一致）
scaler_g = StandardScaler()
df['g_it_std'] = scaler_g.fit_transform(df[['检测孕周']])
bmi_for_rank = df['孕妇BMI'].copy()
df['b_i_norm'] = bmi_for_rank.rank(pct=True)
X_smooth = df[['g_it_std', 'b_i_norm']]
X_linear = df[['年龄', '身高']].astype(float)
X_cat = pd.get_dummies(df[['IVF妊娠', '怀孕次数', '生产次数']], drop_first=True)
X = pd.concat([X_smooth, X_linear, X_cat], axis=1)
y = df['y_it']
X_columns = X.columns
X_cat_columns = X_cat.columns
gam = GAM(
s(0) + s(1) + te(0, 1) + l(2) + l(3) + l(4) + l(5) + l(6),
distribution='gamma',
link='log'
).fit(X, y)
print("GAM（含身高、IVF、怀孕次数、生产次数）模型已成功拟合。")
# ========== 1.5 预测每个孕妇Y染色体浓度首次达到0.04的孕周 ==========
weeks_range = np.arange(14, 31) # 合理孕周范围
threshold = 0.04
reach_week_list = []
for idx, row in df.iterrows():
found = False
for week in weeks_range:
# 构造单条预测数据
g_std = scaler_g.transform(np.array([[week]]))[0, 0]
bmi = row['孕妇BMI']
# 分位数归一化（用全体bmi_for_rank+当前bmi，取最后一个）
combined_bmi = pd.concat([bmi_for_rank, pd.Series([bmi])])
b_norm = combined_bmi.rank(pct=True).iloc[-1]
# 组装特征
x_pred = {
'g_it_std': g_std,
'b_i_norm': b_norm,
'年龄': row['年龄'],
'身高': row['身高']
}
# 类别变量one-hot
for col in X_cat.columns:
if col in row:
x_pred[col] = row[col]
else:
base_col = col.split('_')[0]
val = col[len(base_col)+1:]
x_pred[col] = int(str(row[base_col]) == val)
x_pred_df = pd.DataFrame([x_pred])
x_pred_df = x_pred_df[X.columns]
y_pred = gam.predict(x_pred_df)[0]
if y_pred >= threshold:
reach_week_list.append(week)
found = True
break
if not found:
reach_week_list.append(np.nan)
df['达标孕周'] = reach_week_list
# =========================
# 2. BMI 分组函数
# =========================
def assign_group(bmi_array, cuts):
"""bmi_array: array_like"""
return np.searchsorted(cuts, bmi_array)
# =========================
# 3. 初始化种群
# =========================
def init_population():
population = []
for _ in range(POP_SIZE):
while True:
G_var = np.random.randint(MIN_G, MAX_G+1)
# 生成严格递增的 cut
cuts = np.sort(np.random.uniform(BMI_min, BMI_max - 3*(G_var-1), G_var-1))
cuts = cuts + np.arange(G_var-1)*3 # 保证最小间距为3
if np.all(cuts[1:] - cuts[:-1] >= 3):
week_assign = np.random.choice(weeks_available, G_var)
population.append((cuts, week_assign))
break
return population
# =========================
# 4. 批量预测函数
# =========================
def predict_y_batch(gest_weeks, bmis, ages, heights, ivf_types, preg_times, birth_times, gam, scaler_g, bmi_for_rank, X_columns, X_cat_columns):
gest_weeks = np.array(gest_weeks)
g_it_std = scaler_g.transform(gest_weeks.reshape(-1,1)).flatten()
bmis = np.array(bmis)
combined_bmi = pd.concat([bmi_for_rank, pd.Series(bmis)])
b_i_norm_all = combined_bmi.rank(pct=True).iloc[-len(bmis):].values
# 构造DataFrame
X_pred = pd.DataFrame({
"g_it_std": g_it_std,
"b_i_norm": b_i_norm_all,
"年龄": ages,
"身高": heights
})
# 类别变量one-hot
cat_df = pd.DataFrame({
"IVF妊娠": ivf_types,
"怀孕次数": preg_times,
"生产次数": birth_times
})
X_cat_pred = pd.get_dummies(cat_df, drop_first=True)
X_cat_pred = X_cat_pred.reindex(columns=X_cat_columns, fill_value=0)
X_pred = pd.concat([X_pred, X_cat_pred], axis=1)
X_pred = X_pred[X_columns] # 保证顺序
y_preds = gam.predict(X_pred)
# 假阴性概率
fn_probs = np.ones_like(y_preds)
mask = y_preds > 0.04
fn_probs[mask] = np.exp(-50 * (y_preds[mask] - 0.04))
return y_preds, fn_probs
# =========================
# 5. GA 适应度函数
# =========================
def gc_score(gc):
# 38-42为合格，40-60为警示，其余为异常
if 38 <= gc <= 42:
return 1
elif 40 <= gc <= 60:
return 0.5
else:
return 0
def fitness(cuts, group_week_assign, df, gam, scaler_g, bmi_for_rank, X_columns, X_cat_columns):
bmis = df['孕妇BMI'].values
ages = df['年龄'].values
heights = df['身高'].values
ivf_types = df['IVF妊娠'].values
preg_times = df['怀孕次数'].values
birth_times = df['生产次数'].values
# 每个孕妇所属组
groups = np.array([assign_group(bmi, cuts) for bmi in bmis])
# 每个孕妇对应的检测周
check_weeks = np.array([group_week_assign[g] for g in groups])
# 批量预测 Y 浓度和假阴性概率
_, fn_probs = predict_y_batch(
check_weeks, bmis, ages, heights, ivf_types, preg_times, birth_times,
gam, scaler_g, bmi_for_rank, X_columns, X_cat_columns
)
# 总假阴性风险
total_risk = np.sum(fn_probs)
# 检测周评分
week_score = np.where(check_weeks <= 12, 1, np.where(check_weeks <= 27, 1.5, 5))
total_score = np.sum(week_score)
# BMI 区间惩罚
cut_diffs = np.diff(np.concatenate(([BMI_min], cuts, [BMI_max])))
small_interval_penalty = np.sum(cut_diffs < 3) * 100 # 区间小于3加惩罚
G_var = len(cuts) + 1
# GC含量合理度分数
gc_scores = np.array([gc_score(gc) for gc in df['GC含量']])
mean_gc_score = np.mean(gc_scores)
gc_penalty = 2 * (1 - mean_gc_score) # 你可以调整权重2为更大或更小
# 最终目标函数
obj = 10 * total_risk + 2.0 * total_score + GROUP_PENALTY * G_var + small_interval_penalty + gc_penalty
return obj
# =========================
# 6. 变异操作
# =========================
def mutate(indiv):
cuts, week_assign = indiv
cuts = cuts.copy()
week_assign = week_assign.copy()
for i in range(len(cuts)):
if np.random.rand() < MUT_RATE:
cuts[i] += np.random.uniform(-2, 2)
cuts = np.clip(cuts, BMI_min, BMI_max)
cuts = np.sort(cuts)
# 修正最小间距
for i in range(1, len(cuts)):
if cuts[i] - cuts[i-1] < 3:
cuts[i] = cuts[i-1] + 3
# 最后一组不超过 BMI_max
if cuts[-1] > BMI_max:
cuts[-1] = BMI_max
for i in range(len(week_assign)):
if np.random.rand() < MUT_RATE:
week_assign[i] = np.random.choice(weeks_available)
return (cuts, week_assign)
# =========================
# 6.5 交叉操作
# =========================
def crossover(parent1, parent2):
cuts1, weeks1 = parent1
cuts2, weeks2 = parent2
# 只允许同组数的个体交叉
if len(cuts1) != len(cuts2) or len(weeks1) != len(weeks2):
return parent1 # 不交叉，直接返回parent1
if len(cuts1) < 2 or len(weeks1) < 2:
return parent1 # 只有1个cut或1个week时无法交叉
# 随机选择一个交叉点
point = np.random.randint(1, len(cuts1))
new_cuts = np.concatenate([cuts1[:point], cuts2[point:]])
point_w = np.random.randint(1, len(weeks1))
new_weeks = np.concatenate([weeks1[:point_w], weeks2[point_w:]])
return (new_cuts, new_weeks)
# =========================
# 7. GA 主循环（向量化适应度）
# =========================
population = init_population()
best_obj = float('inf')
best_solution = None
for gen in range(N_GEN):
# 适应度计算
fitness_values = [
fitness(cuts, weeks, df, gam, scaler_g, bmi_for_rank, X_columns, X_cat_columns)
for cuts, weeks in population
]
# 更新最优解
idx = np.argmin(fitness_values)
if fitness_values[idx] < best_obj:
best_obj = fitness_values[idx]
best_solution = population[idx]
# 选择 + 变异
weights = np.exp(-np.array(fitness_values)) # 适应度越小概率越大
new_population = []
for _ in range(POP_SIZE):
if np.random.rand() < CROSS_RATE:
# 交叉
parent1 = random.choices(population, weights=weights, k=1)[0]
parent2 = random.choices(population, weights=weights, k=1)[0]
child = crossover(parent1, parent2)
else:
# 只变异
parent = random.choices(population, weights=weights, k=1)[0]
child = parent
# 变异
child = mutate(child)
new_population.append(child)
population = new_population
if gen % 10 == 0:
print(f"Generation {gen}: best objective = {best_obj:.4f}")
# =========================
# 8. 输出最优分组结果
# =========================
best_cuts, best_weeks = best_solution
print("\n最优 BMI 分组 cut:", best_cuts)
print("最优每组检测周:", best_weeks)
# 记录特征列顺序
X_columns = X.columns
X_cat_columns = X_cat.columns
# 计算每个孕妇的分组和分组检测周
bmis = df['孕妇BMI'].values
ages = df['年龄'].values
heights = df['身高'].values
ivf_types = df['IVF妊娠'].values
preg_times = df['怀孕次数'].values
birth_times = df['生产次数'].values
groups = np.array([assign_group(bmi, best_cuts) for bmi in bmis])
check_weeks = np.array([best_weeks[g] for g in groups])
y_preds, fn_probs = predict_y_batch(
check_weeks, bmis, ages, heights, ivf_types, preg_times, birth_times,
gam, scaler_g, bmi_for_rank, X_columns, X_cat_columns
)
# 输出示例
for code, bmi, week, y_pred, fn_prob in zip(df['孕妇代码'], bmis, check_weeks, y_preds, fn_probs):
print(f"孕妇 {code}: BMI={bmi:.2f}, 分组检测周={week}, 预测 Y={y_pred:.6f}, 假阴性={fn_prob:.4f}")
def cuts_to_intervals(cuts, BMI_min, BMI_max):
"""
cuts: array-like, BMI 分割点
BMI_min, BMI_max: BMI 范围
返回每个组的区间 [(下限, 上限), ...]
"""
intervals = []
prev = BMI_min
for c in cuts:
intervals.append((prev, c))
prev = c
intervals.append((prev, BMI_max)) # 最后一组
return intervals
best_intervals = cuts_to_intervals(best_cuts, BMI_min, BMI_max)
for i, (low, high) in enumerate(best_intervals):
print(f"组{i+1}: BMI {low:.1f}-{high:.1f}, 检测周: {best_weeks[i]}")
GAM（含身高、IVF、怀孕次数、生产次数）模型已成功拟合。 Generation 0: best objective = 333.2560 Generation 10: best objective = 332.1364 Generation 20: best objective = 332.1364 Generation 30: best objective = 332.1364 Generation 40: best objective = 332.1364 Generation 50: best objective = 332.1364 Generation 60: best objective = 332.1364 Generation 70: best objective = 332.1364 Generation 80: best objective = 332.1364 Generation 90: best objective = 332.1364 Generation 100: best objective = 332.1364 Generation 110: best objective = 332.1364 Generation 120: best objective = 332.1364 Generation 130: best objective = 332.1364 Generation 140: best objective = 332.1364 Generation 150: best objective = 332.1364 Generation 160: best objective = 332.1364 Generation 170: best objective = 332.1364 Generation 180: best objective = 332.1364 Generation 190: best objective = 332.1364 最优 BMI 分组 cut: [32.78289596 39.45134238 42.45134238] 最优每组检测周: [27 27 18 14] 孕妇 A003: BMI=30.74, 分组检测周=27, 预测 Y=0.065554, 假阴性=0.2787 孕妇 A008: BMI=29.55, 分组检测周=27, 预测 Y=0.082078, 假阴性=0.1220 孕妇 A010: BMI=35.56, 分组检测周=27, 预测 Y=0.088407, 假阴性=0.0889 孕妇 A012: BMI=30.80, 分组检测周=27, 预测 Y=0.060238, 假阴性=0.3635 孕妇 A018: BMI=28.13, 分组检测周=27, 预测 Y=0.068960, 假阴性=0.2350 孕妇 A023: BMI=30.04, 分组检测周=27, 预测 Y=0.073612, 假阴性=0.1863 孕妇 A026: BMI=28.04, 分组检测周=27, 预测 Y=0.060390, 假阴性=0.3608 孕妇 A029: BMI=33.11, 分组检测周=27, 预测 Y=0.083477, 假阴性=0.1137 孕妇 A041: BMI=36.81, 分组检测周=27, 预测 Y=0.065172, 假阴性=0.2840 孕妇 A048: BMI=35.23, 分组检测周=27, 预测 Y=0.069061, 假阴性=0.2339 孕妇 A055: BMI=35.00, 分组检测周=27, 预测 Y=0.058902, 假阴性=0.3886 孕妇 A061: BMI=36.79, 分组检测周=27, 预测 Y=0.076439, 假阴性=0.1617 孕妇 A065: BMI=32.05, 分组检测周=27, 预测 Y=0.064489, 假阴性=0.2939 孕妇 A066: BMI=31.22, 分组检测周=27, 预测 Y=0.081847, 假阴性=0.1234 孕妇 A069: BMI=34.18, 分组检测周=27, 预测 Y=0.081891, 假阴性=0.1231 孕妇 A071: BMI=34.72, 分组检测周=27, 预测 Y=0.058010, 假阴性=0.4064 孕妇 A073: BMI=34.55, 分组检测周=27, 预测 Y=0.072301, 假阴性=0.1989 孕妇 A092: BMI=30.48, 分组检测周=27, 预测 Y=0.052109, 假阴性=0.5458 孕妇 A095: BMI=39.14, 分组检测周=27, 预测 Y=0.063035, 假阴性=0.3161 孕妇 A096: BMI=27.90, 分组检测周=27, 预测 Y=0.073302, 假阴性=0.1892 孕妇 A105: BMI=32.41, 分组检测周=27, 预测 Y=0.059893, 假阴性=0.3698 孕妇 A109: BMI=29.41, 分组检测周=27, 预测 Y=0.064052, 假阴性=0.3004 孕妇 A130: BMI=30.10, 分组检测周=27, 预测 Y=0.069051, 假阴性=0.2340 孕妇 A133: BMI=28.27, 分组检测周=27, 预测 Y=0.050963, 假阴性=0.5780 孕妇 A146: BMI=31.64, 分组检测周=27, 预测 Y=0.078399, 假阴性=0.1466 孕妇 A147: BMI=31.41, 分组检测周=27, 预测 Y=0.076557, 假阴性=0.1608 孕妇 A155: BMI=30.47, 分组检测周=27, 预测 Y=0.067779, 假阴性=0.2493 孕妇 A163: BMI=33.31, 分组检测周=27, 预测 Y=0.057811, 假阴性=0.4104 孕妇 A169: BMI=31.35, 分组检测周=27, 预测 Y=0.080471, 假阴性=0.1322 孕妇 A170: BMI=29.44, 分组检测周=27, 预测 Y=0.066267, 假阴性=0.2689 孕妇 A171: BMI=32.70, 分组检测周=27, 预测 Y=0.075103, 假阴性=0.1729 孕妇 A172: BMI=33.74, 分组检测周=27, 预测 Y=0.088822, 假阴性=0.0871 孕妇 A175: BMI=31.83, 分组检测周=27, 预测 Y=0.074203, 假阴性=0.1808 孕妇 A176: BMI=30.50, 分组检测周=27, 预测 Y=0.081633, 假阴性=0.1247 孕妇 A178: BMI=34.57, 分组检测周=27, 预测 Y=0.062392, 假阴性=0.3264 孕妇 A192: BMI=33.29, 分组检测周=27, 预测 Y=0.083478, 假阴性=0.1137 孕妇 A194: BMI=34.65, 分组检测周=27, 预测 Y=0.067386, 假阴性=0.2543 孕妇 A196: BMI=31.36, 分组检测周=27, 预测 Y=0.081451, 假阴性=0.1259 孕妇 A198: BMI=30.65, 分组检测周=27, 预测 Y=0.076572, 假阴性=0.1606 孕妇 A200: BMI=32.09, 分组检测周=27, 预测 Y=0.078184, 假阴性=0.1482 孕妇 A201: BMI=31.26, 分组检测周=27, 预测 Y=0.072366, 假阴性=0.1982 孕妇 A202: BMI=31.68, 分组检测周=27, 预测 Y=0.075061, 假阴性=0.1732 孕妇 A203: BMI=28.17, 分组检测周=27, 预测 Y=0.075080, 假阴性=0.1731 孕妇 A204: BMI=33.73, 分组检测周=27, 预测 Y=0.077985, 假阴性=0.1497 孕妇 A205: BMI=29.58, 分组检测周=27, 预测 Y=0.072238, 假阴性=0.1995 孕妇 A215: BMI=32.72, 分组检测周=27, 预测 Y=0.078262, 假阴性=0.1476 孕妇 A216: BMI=30.28, 分组检测周=27, 预测 Y=0.071650, 假阴性=0.2055 孕妇 A219: BMI=29.01, 分组检测周=27, 预测 Y=0.080942, 假阴性=0.1291 孕妇 A224: BMI=32.19, 分组检测周=27, 预测 Y=0.085256, 假阴性=0.1041 孕妇 A227: BMI=33.01, 分组检测周=27, 预测 Y=0.071671, 假阴性=0.2052 孕妇 A230: BMI=34.21, 分组检测周=27, 预测 Y=0.074238, 假阴性=0.1805 孕妇 A233: BMI=32.17, 分组检测周=27, 预测 Y=0.074446, 假阴性=0.1787 孕妇 A234: BMI=34.75, 分组检测周=27, 预测 Y=0.075543, 假阴性=0.1691 孕妇 A235: BMI=34.27, 分组检测周=27, 预测 Y=0.079511, 假阴性=0.1387 孕妇 A236: BMI=32.76, 分组检测周=27, 预测 Y=0.070178, 假阴性=0.2212 孕妇 A242: BMI=30.03, 分组检测周=27, 预测 Y=0.075789, 假阴性=0.1671 孕妇 A245: BMI=30.76, 分组检测周=27, 预测 Y=0.081542, 假阴性=0.1253 孕妇 A247: BMI=33.21, 分组检测周=27, 预测 Y=0.084295, 假阴性=0.1092 孕妇 A248: BMI=33.10, 分组检测周=27, 预测 Y=0.073970, 假阴性=0.1830 孕妇 A250: BMI=33.57, 分组检测周=27, 预测 Y=0.077811, 假阴性=0.1510 孕妇 A255: BMI=32.14, 分组检测周=27, 预测 Y=0.069799, 假阴性=0.2254 孕妇 A256: BMI=33.34, 分组检测周=27, 预测 Y=0.079830, 假阴性=0.1365 孕妇 A263: BMI=29.35, 分组检测周=27, 预测 Y=0.075910, 假阴性=0.1660 孕妇 A266: BMI=32.97, 分组检测周=27, 预测 Y=0.073809, 假阴性=0.1844 孕妇 A267: BMI=30.70, 分组检测周=27, 预测 Y=0.080291, 假阴性=0.1334 组1: BMI 26.0-32.8, 检测周: 27 组2: BMI 32.8-39.5, 检测周: 27 组3: BMI 39.5-42.5, 检测周: 18 组4: BMI 42.5-47.0, 检测周: 14
代码版本2¶
In [4]:
import pandas as pd
import numpy as np
import random
from sklearn.preprocessing import StandardScaler
from pygam import GAM, s, l, te
# =========================
# 1. GA 参数
# =========================
SEED = 42
random.seed(SEED)
np.random.seed(SEED)
MAX_G = 6
MIN_G = 2
GROUP_PENALTY_WEIGHT = -0.2 # 惩罚项权重，可以调整
POP_SIZE = 200
N_GEN = 50
MUT_RATE = 0.7
CROSS_RATE = 0.7
# BMI 和孕周范围
BMI_min, BMI_max = 26, 47
weeks_available = np.arange(14, 30)
# =========================
# 2. 辅助函数
# =========================
def assign_group(bmi_array, cuts):
"""根据BMI分割点，为每个孕妇分配组别。"""
# np.searchsorted 返回插入点，用于确定元素在有序数组中的位置
# 0: BMI < cuts[0]
# 1: cuts[0] <= BMI < cuts[1]
# ...
# len(cuts): BMI >= cuts[-1]
return np.searchsorted(cuts, bmi_array)
def gc_score(gc):
"""根据GC含量合理度评分。"""
if 38 <= gc <= 42:
return 1
elif 40 <= gc <= 60:
return 0.5
else:
return 0
# =========================
# 3. 优化后的初始化种群
# =========================
def init_population():
"""初始化遗传算法种群，每个个体包含 BMI 分割点和对应的检测孕周。"""
population = []
for _ in range(POP_SIZE):
G_var = np.random.randint(MIN_G, MAX_G + 1)
# 优化：生成 G_var-1 个随机数，然后排序并强制添加间距
# 先生成 G_var-1 个点，确保它们在 BMI_min 到 BMI_max 范围内
cuts = np.sort(np.random.uniform(BMI_min, BMI_max, G_var - 1))
# 强制最小间距，确保组别划分有效
# 从第二个分割点开始，如果与前一个点的距离小于3，则将其调整到与前一个点相距3
for i in range(1, len(cuts)):
if cuts[i] - cuts[i-1] < 3:
cuts[i] = cuts[i-1] + 3
# 确保最后一个 cut 在合理范围内，如果调整后超过了 BMI_max，则将其设为 BMI_max
if cuts[-1] > BMI_max:
cuts[-1] = BMI_max
# 为每个分组随机选择一个检测孕周
week_assign = np.random.choice(weeks_available, G_var)
population.append({'cuts': cuts, 'weeks': week_assign})
return population
# =========================
# 4. 批量预测函数 (与原代码一致，但更清晰)
# =========================
def predict_y_batch(gest_weeks, bmis, ages, heights, ivf_types, preg_times, birth_times, gam, scaler_g, bmi_for_rank, X_columns, X_cat_columns):
"""
批量预测 Y 染色体浓度和假阴性概率。
假阴性概率公式已根据你的原代码进行保留，但请注意其经验性。
"""
gest_weeks = np.array(gest_weeks)
# 标准化检测孕周
g_it_std = scaler_g.transform(gest_weeks.reshape(-1, 1)).flatten()
bmis = np.array(bmis)
# 计算 BMI 的分位数排名，结合了训练集和当前预测集的数据
combined_bmi = pd.concat([bmi_for_rank, pd.Series(bmis)])
b_i_norm_all = combined_bmi.rank(pct=True).iloc[-len(bmis):].values
# 构建预测用的 DataFrame
X_pred = pd.DataFrame({
"g_it_std": g_it_std,
"b_i_norm": b_i_norm_all,
"年龄": ages,
"身高": heights
})
# 处理类别变量，并确保列顺序与训练时一致
cat_df = pd.DataFrame({
"IVF妊娠": ivf_types,
"怀孕次数": preg_times,
"生产次数": birth_times
})
X_cat_pred = pd.get_dummies(cat_df, drop_first=True)
# 确保列存在，不存在的列填充0，保持列结构一致
X_cat_pred = X_cat_pred.reindex(columns=X_cat_columns, fill_value=0)
# 合并数值型和类别型特征
X_pred = pd.concat([X_pred, X_cat_pred], axis=1)
# 确保特征列顺序与训练时GAM模型期望的一致
X_pred = X_pred[X_columns]
# 使用 GAM 模型进行预测
y_preds = gam.predict(X_pred)
# ==========================
# 添加高斯噪声
# ==========================
noise = np.random.normal(loc=0, scale=0.0005, size=y_preds.shape)
y_preds_noisy = y_preds # + noise
# ==========================
# 计算假阴性概率 (fn_probs)
fn_probs = np.ones_like(y_preds_noisy)
mask = y_preds_noisy > 0.04 # 只对预测 Y > 0.04 的情况计算指数衰减
# 指数衰减公式，预测 Y 越高，假阴性概率越低
fn_probs[mask] = np.exp(-50 * (y_preds_noisy[mask] - 0.04))
return y_preds_noisy, fn_probs
# =========================
# 5. 适应度函数
# =========================
def risk_score(weeks, bmis):
"""
将孕周分为三档，并根据BMI进行额外惩罚，以鼓励更早的孕周。
"""
weeks = np.array(weeks)
bmis = np.array(bmis)
score = np.zeros_like(weeks, dtype=float)
# 基础分
mask_early = weeks <= 12
score[mask_early] = 1 + 0.05 * weeks[mask_early]
mask_mid = (weeks >= 13) & (weeks <= 27)
score[mask_mid] = 2 + 0.02 * (weeks[mask_mid] - 13)
mask_late = weeks >= 28
score[mask_late] = 3 + 0.01 * (weeks[mask_late] - 28)
# 针对高BMI的额外惩罚（阈值可调）
mask_high_bmi = bmis > 35
score[mask_high_bmi] += (bmis[mask_high_bmi] - 35) * 0.1 # BMI越高，额外惩罚越大
return score
def fitness(cuts, group_week_assign, df, gam, scaler_g, bmi_for_rank, X_columns, X_cat_columns):
"""
计算个体（一组 BMI 分割点和孕周）的适应度。
目标是最小化此值。
"""
# 提取当前个体需要的所有孕妇数据
bmis = df['孕妇BMI'].values
ages = df['年龄'].values
heights = df['身高'].values
ivf_types = df['IVF妊娠'].values
preg_times = df['怀孕次数'].values
birth_times = df['生产次数'].values
# 根据个体的 BMI 分割点，为所有孕妇分配组别
groups = assign_group(bmis, cuts)
# 根据分配的组别，查找对应的检测孕周
check_weeks = np.array([group_week_assign[g] for g in groups])
# 使用批量预测函数计算假阴性概率
# 注意：这里使用的是添加了噪声的预测值
_, fn_probs = predict_y_batch(
gest_weeks=check_weeks,
bmis=bmis,
ages=ages,
heights=heights,
ivf_types=ivf_types,
preg_times=preg_times,
birth_times=birth_times,
gam=gam,
scaler_g=scaler_g,
bmi_for_rank=bmi_for_rank,
X_columns=X_columns,
X_cat_columns=X_cat_columns
)
# 计算总的假阴性风险（假阴性概率的总和）
total_risk = np.sum(fn_probs)
# ========================
# 检测周评分，使用上一题的 risk_score 函数（区间递增）
# ========================
# 调用修改后的 risk_score 函数，传入 bmis
week_score = risk_score(check_weeks, bmis)
total_score = np.sum(week_score)
# BMI 区间长度的惩罚项
cut_diffs = np.diff(np.concatenate(([BMI_min], cuts, [BMI_max])))
small_interval_penalty = np.sum(cut_diffs < 3) * 100
# 当前个体的分组数量
G_var = len(cuts) + 1
# GC 含量惩罚
gc_scores = np.array([gc_score(gc) for gc in df['GC含量']])
mean_gc_score = np.mean(gc_scores)
gc_penalty = 2 * (1 - mean_gc_score)
# 最终目标函数
obj = 10 * total_risk + 2 * total_score + GROUP_PENALTY_WEIGHT * G_var + small_interval_penalty + gc_penalty
return obj
# =========================
# 6. 变异操作
# =========================
def mutate(indiv):
"""
优化：变异操作增加随机性，不仅改变值，还可以改变组数。
"""
cuts = indiv['cuts'].copy()
weeks = indiv['weeks'].copy()
# 变异类型：
# 1. 变异 BMI 分割点的值
if np.random.rand() < MUT_RATE:
# 随机调整每个分割点的值，幅度在 -5 到 5 之间
cuts += np.random.uniform(-5, 5, len(cuts))
# 裁剪分割点到 BMI_min 和 BMI_max 范围内
cuts = np.clip(cuts, BMI_min, BMI_max)
# 重新排序分割点
cuts = np.sort(cuts)
# 修正最小间距（与初始化时类似）
for i in range(1, len(cuts)):
if cuts[i] - cuts[i-1] < 3:
cuts[i] = cuts[i-1] + 3
# 移除重复的分割点，确保唯一性
cuts = np.unique(cuts)
# 2. 变异孕周分配
if np.random.rand() < MUT_RATE:
for i in range(len(weeks)):
# 为每个分组随机选择一个新的检测孕周
weeks[i] = np.random.choice(weeks_available)
# 3. 变异组数（动态调整 G_var）
# 增加一个分组（添加一个 cut）
if np.random.rand() < MUT_RATE and len(cuts) + 1 < MAX_G:
# 在现有分割点之间随机选择一个位置插入新的分割点
new_cut = np.random.uniform(cuts.min(), cuts.max())
cuts = np.sort(np.append(cuts, new_cut))
# 减少一个分组（移除一个 cut）
if np.random.rand() < MUT_RATE and len(cuts) + 1 > MIN_G:
# 随机移除一个现有的分割点
cuts = np.delete(cuts, np.random.randint(len(cuts)))
# 确保 weeks 数组的长度与 cuts 数组的长度（加1）匹配
weeks_len = len(cuts) + 1
if len(weeks) != weeks_len:
# 如果长度不匹配，则重新随机生成 weeks
weeks = np.random.choice(weeks_available, weeks_len)
return {'cuts': cuts, 'weeks': weeks}
# =========================
# 7. 交叉操作
# =========================
def crossover(parent1, parent2):
"""
优化：允许不同组数的个体进行交叉。
"""
cuts1, weeks1 = parent1['cuts'], parent1['weeks']
cuts2, weeks2 = parent2['cuts'], parent2['weeks']
# 将两个父代的 cuts 和 weeks 合并，然后重新排序和选择
# 合并两个父代的分割点，并去重，得到一个更全面的分割点集合
all_cuts = np.unique(np.concatenate([cuts1, cuts2]))
# 随机确定子代的组数，允许在 MIN_G-1 到 MAX_G-1 之间随机选择 cut 的数量
# 注意：这里的 new_cuts_size 是 cut 的数量，实际组数是 new_cuts_size + 1
new_cuts_size = np.random.randint(MIN_G - 1, MAX_G)
# 从合并后的所有分割点中，随机选择 new_cuts_size 个作为子代的分割点
if len(all_cuts) > new_cuts_size:
new_cuts = np.sort(np.random.choice(all_cuts, new_cuts_size, replace=False))
else:
# 如果所有分割点少于或等于目标数量，则全部保留
new_cuts = all_cuts
# 为子代随机分配检测孕周，数量为 new_cuts_size + 1
new_weeks = np.random.choice(weeks_available, len(new_cuts) + 1)
return {'cuts': new_cuts, 'weeks': new_weeks}
# =========================
# 8. GA 主循环
# =========================
population = init_population()
best_obj = float('inf')
best_solution = None
for gen in range(N_GEN):
# 适应度计算：为种群中的每个个体计算适应度值
fitness_values = [
fitness(indiv['cuts'], indiv['weeks'], df, gam, scaler_g, bmi_for_rank, X_columns, X_cat_columns)
for indiv in population
]
# 更新最优解：记录当前最优的个体及其适应度值
idx = np.argmin(fitness_values)
if fitness_values[idx] < best_obj:
best_obj = fitness_values[idx]
best_solution = population[idx]
# 选择：使用轮盘赌选择法，适应度值越高的个体（目标函数值越小）被选中的概率越大
# 转换为概率，适应度值越小，权重越大
weights = np.exp(-np.array(fitness_values))
weights /= weights.sum() # 归一化权重
new_population = []
# 精英策略：将当前最优解直接加入下一代种群，保证最优解不丢失
new_population.append(best_solution)
# 繁殖下一代种群
for _ in range(POP_SIZE - 1): # 减去精英个体数量
# 交叉操作
if np.random.rand() < CROSS_RATE:
# 从当前种群中按权重随机选择两个父代
parent1 = random.choices(population, weights=weights, k=1)[0]
parent2 = random.choices(population, weights=weights, k=1)[0]
child = crossover(parent1, parent2)
else:
# 如果不进行交叉，则随机选择一个父代（相当于拷贝）
parent = random.choices(population, weights=weights, k=1)[0]
child = parent
# 变异操作
child = mutate(child)
new_population.append(child)
# 更新种群
population = new_population
# 每隔10代打印一次当前最优的适应度值
if gen % 10 == 0:
print(f"Generation {gen}: best objective = {best_obj:.4f}")
# =========================
# 9. 输出最优分组结果
# =========================
best_cuts, best_weeks = best_solution['cuts'], best_solution['weeks']
print("\n最优 BMI 分组 cut:", best_cuts)
print("最优每组检测周:", best_weeks)
# 重新计算并打印详细结果，用于验证
bmis = df['孕妇BMI'].values
# 为所有孕妇根据最优分割点分配组别
groups = assign_group(bmis, best_cuts)
# 获取每个孕妇对应的最优检测周
check_weeks = np.array([best_weeks[g] for g in groups])
# 使用最优分组方案进行预测
# 注意：这里的 predict_y_batch 函数已经包含了噪声
y_preds_noisy, fn_probs = predict_y_batch(
check_weeks, bmis, df['年龄'].values, df['身高'].values, df['IVF妊娠'].values,
df['怀孕次数'].values, df['生产次数'].values,
gam, scaler_g, bmi_for_rank, X_columns, X_cat_columns
)
# 打印每个孕妇的详细信息
for code, bmi, week, y_pred, fn_prob in zip(df['孕妇代码'], bmis, check_weeks, y_preds_noisy, fn_probs):
print(f"孕妇 {code}: BMI={bmi:.2f}, 分组检测周={week}, 预测 Y={y_pred:.6f}, 假阴性={fn_prob:.4f}")
# 辅助函数，将分割点转换为区间表示
def cuts_to_intervals(cuts, BMI_min, BMI_max):
"""
将 BMI 分割点转换为可读的区间表示。
"""
intervals = []
prev = BMI_min
for c in cuts:
intervals.append((prev, c))
prev = c
intervals.append((prev, BMI_max)) # 最后一组区间
return intervals
# 打印最终确定的 BMI 分组区间和对应的检测周
best_intervals = cuts_to_intervals(best_cuts, BMI_min, BMI_max)
for i, (low, high) in enumerate(best_intervals):
print(f"组{i+1}: BMI {low:.1f}-{high:.1f}, 检测周: {best_weeks[i]}")
Generation 0: best objective = 411.0337 Generation 10: best objective = 410.7823 Generation 20: best objective = 410.7823 Generation 30: best objective = 410.7823 Generation 40: best objective = 410.7823 最优 BMI 分组 cut: [36.72817988 39.72817988 42.72817988] 最优每组检测周: [14 16 26 23] 孕妇 A003: BMI=30.74, 分组检测周=14, 预测 Y=0.064980, 假阴性=0.2868 孕妇 A008: BMI=29.55, 分组检测周=14, 预测 Y=0.080113, 假阴性=0.1346 孕妇 A010: BMI=35.56, 分组检测周=14, 预测 Y=0.086613, 假阴性=0.0972 孕妇 A012: BMI=30.80, 分组检测周=14, 预测 Y=0.059882, 假阴性=0.3701 孕妇 A018: BMI=28.13, 分组检测周=14, 预测 Y=0.067807, 假阴性=0.2490 孕妇 A023: BMI=30.04, 分组检测周=14, 预测 Y=0.071955, 假阴性=0.2023 孕妇 A026: BMI=28.04, 分组检测周=14, 预测 Y=0.059429, 假阴性=0.3785 孕妇 A029: BMI=33.11, 分组检测周=14, 预测 Y=0.081796, 假阴性=0.1237 孕妇 A041: BMI=36.81, 分组检测周=16, 预测 Y=0.063958, 假阴性=0.3018 孕妇 A048: BMI=35.23, 分组检测周=14, 预测 Y=0.067729, 假阴性=0.2500 孕妇 A055: BMI=35.00, 分组检测周=14, 预测 Y=0.057825, 假阴性=0.4101 孕妇 A061: BMI=36.79, 分组检测周=16, 预测 Y=0.075065, 假阴性=0.1732 孕妇 A065: BMI=32.05, 分组检测周=14, 预测 Y=0.063831, 假阴性=0.3038 孕妇 A066: BMI=31.22, 分组检测周=14, 预测 Y=0.081451, 假阴性=0.1259 孕妇 A069: BMI=34.18, 分组检测周=14, 预测 Y=0.080705, 假阴性=0.1306 孕妇 A071: BMI=34.72, 分组检测周=14, 预测 Y=0.057060, 假阴性=0.4261 孕妇 A073: BMI=34.55, 分组检测周=14, 预测 Y=0.071248, 假阴性=0.2096 孕妇 A092: BMI=30.48, 分组检测周=14, 预测 Y=0.051266, 假阴性=0.5693 孕妇 A095: BMI=39.14, 分组检测周=16, 预测 Y=0.061827, 假阴性=0.3358 孕妇 A096: BMI=27.90, 分组检测周=14, 预测 Y=0.072168, 假阴性=0.2002 孕妇 A105: BMI=32.41, 分组检测周=14, 预测 Y=0.058695, 假阴性=0.3927 孕妇 A109: BMI=29.41, 分组检测周=14, 预测 Y=0.062601, 假阴性=0.3230 孕妇 A130: BMI=30.10, 分组检测周=14, 预测 Y=0.067583, 假阴性=0.2518 孕妇 A133: BMI=28.27, 分组检测周=14, 预测 Y=0.049993, 假阴性=0.6067 孕妇 A146: BMI=31.64, 分组检测周=14, 预测 Y=0.078001, 假阴性=0.1496 孕妇 A147: BMI=31.41, 分组检测周=14, 预测 Y=0.076240, 假阴性=0.1633 孕妇 A155: BMI=30.47, 分组检测周=14, 预测 Y=0.066557, 假阴性=0.2650 孕妇 A163: BMI=33.31, 分组检测周=14, 预测 Y=0.056805, 假阴性=0.4316 孕妇 A169: BMI=31.35, 分组检测周=14, 预测 Y=0.080174, 假阴性=0.1342 孕妇 A170: BMI=29.44, 分组检测周=14, 预测 Y=0.064711, 假阴性=0.2907 孕妇 A171: BMI=32.70, 分组检测周=14, 预测 Y=0.073515, 假阴性=0.1872 孕妇 A172: BMI=33.74, 分组检测周=14, 预测 Y=0.087507, 假阴性=0.0930 孕妇 A175: BMI=31.83, 分组检测周=14, 预测 Y=0.073594, 假阴性=0.1864 孕妇 A176: BMI=30.50, 分组检测周=14, 预测 Y=0.080470, 假阴性=0.1322 孕妇 A178: BMI=34.57, 分组检测周=14, 预测 Y=0.061457, 假阴性=0.3420 孕妇 A192: BMI=33.29, 分组检测周=14, 预测 Y=0.081949, 假阴性=0.1228 孕妇 A194: BMI=34.65, 分组检测周=14, 预测 Y=0.066335, 假阴性=0.2680 孕妇 A196: BMI=31.36, 分组检测周=14, 预测 Y=0.081149, 假阴性=0.1278 孕妇 A198: BMI=30.65, 分组检测周=14, 预测 Y=0.075629, 假阴性=0.1684 孕妇 A200: BMI=32.09, 分组检测周=14, 预测 Y=0.077220, 假阴性=0.1555 孕妇 A201: BMI=31.26, 分组检测周=14, 预测 Y=0.072069, 假阴性=0.2012 孕妇 A202: BMI=31.68, 分组检测周=14, 预测 Y=0.074576, 假阴性=0.1775 孕妇 A203: BMI=28.17, 分组检测周=14, 预测 Y=0.073745, 假阴性=0.1850 孕妇 A204: BMI=33.73, 分组检测周=14, 预测 Y=0.076794, 假阴性=0.1589 孕妇 A205: BMI=29.58, 分组检测周=14, 预测 Y=0.070511, 假阴性=0.2175 孕妇 A215: BMI=32.72, 分组检测周=14, 预测 Y=0.076553, 假阴性=0.1608 孕妇 A216: BMI=30.28, 分组检测周=14, 预测 Y=0.070235, 假阴性=0.2205 孕妇 A219: BMI=29.01, 分组检测周=14, 预测 Y=0.079296, 假阴性=0.1402 孕妇 A224: BMI=32.19, 分组检测周=14, 预测 Y=0.083685, 假阴性=0.1126 孕妇 A227: BMI=33.01, 分组检测周=14, 预测 Y=0.070125, 假阴性=0.2217 孕妇 A230: BMI=34.21, 分组检测周=14, 预测 Y=0.073174, 假阴性=0.1904 孕妇 A233: BMI=32.17, 分组检测周=14, 预测 Y=0.073213, 假阴性=0.1900 孕妇 A234: BMI=34.75, 分组检测周=14, 预测 Y=0.074236, 假阴性=0.1805 孕妇 A235: BMI=34.27, 分组检测周=14, 预测 Y=0.078369, 假阴性=0.1468 孕妇 A236: BMI=32.76, 分组检测周=14, 预测 Y=0.068628, 假阴性=0.2390 孕妇 A242: BMI=30.03, 分组检测周=14, 预测 Y=0.074014, 假阴性=0.1826 孕妇 A245: BMI=30.76, 分组检测周=14, 预测 Y=0.080954, 假阴性=0.1290 孕妇 A247: BMI=33.21, 分组检测周=14, 预测 Y=0.082673, 假阴性=0.1184 孕妇 A248: BMI=33.10, 分组检测周=14, 预测 Y=0.072421, 假阴性=0.1977 孕妇 A250: BMI=33.57, 分组检测周=14, 预测 Y=0.076577, 假阴性=0.1606 孕妇 A255: BMI=32.14, 分组检测周=14, 预测 Y=0.068788, 假阴性=0.2371 孕妇 A256: BMI=33.34, 分组检测周=14, 预测 Y=0.078507, 假阴性=0.1458 孕妇 A263: BMI=29.35, 分组检测周=14, 预测 Y=0.074273, 假阴性=0.1802 孕妇 A266: BMI=32.97, 分组检测周=14, 预测 Y=0.072187, 假阴性=0.2000 孕妇 A267: BMI=30.70, 分组检测周=14, 预测 Y=0.079450, 假阴性=0.1391 组1: BMI 26.0-36.7, 检测周: 14 组2: BMI 36.7-39.7, 检测周: 16 组3: BMI 39.7-42.7, 检测周: 26 组4: BMI 42.7-47.0, 检测周: 23
代码版本3¶
In [7]:
import pandas as pd
import numpy as np
import random
from sklearn.preprocessing import StandardScaler
from pygam import GAM, s, l, te
# =========================
# 1. GA 参数
# =========================
SEED = 42
random.seed(SEED)
np.random.seed(SEED)
MAX_G = 6
MIN_G = 2
GROUP_PENALTY_WEIGHT = -0.5 # 惩罚项权重，可以调整
POP_SIZE = 200
N_GEN = 50
MUT_RATE = 0.7
CROSS_RATE = 0.7
# BMI 和孕周范围
BMI_min, BMI_max = 26, 47
weeks_available = np.arange(14, 30)
# =========================
# 2. 辅助函数
# =========================
def assign_group(bmi_array, cuts):
"""根据BMI分割点，为每个孕妇分配组别。"""
# np.searchsorted 返回插入点，用于确定元素在有序数组中的位置
# 0: BMI < cuts[0]
# 1: cuts[0] <= BMI < cuts[1]
# ...
# len(cuts): BMI >= cuts[-1]
return np.searchsorted(cuts, bmi_array)
def gc_score(gc):
"""根据GC含量合理度评分。"""
if 38 <= gc <= 42:
return 1
elif 40 <= gc <= 60:
return 0.5
else:
return 0
# =========================
# 3. 优化后的初始化种群
# =========================
def init_population():
"""初始化遗传算法种群，每个个体包含 BMI 分割点和对应的检测孕周。"""
population = []
for _ in range(POP_SIZE):
G_var = np.random.randint(MIN_G, MAX_G + 1)
# 优化：生成 G_var-1 个随机数，然后排序并强制添加间距
# 先生成 G_var-1 个点，确保它们在 BMI_min 到 BMI_max 范围内
cuts = np.sort(np.random.uniform(BMI_min, BMI_max, G_var - 1))
# 强制最小间距，确保组别划分有效
# 从第二个分割点开始，如果与前一个点的距离小于3，则将其调整到与前一个点相距3
for i in range(1, len(cuts)):
if cuts[i] - cuts[i-1] < 3:
cuts[i] = cuts[i-1] + 3
# 确保最后一个 cut 在合理范围内，如果调整后超过了 BMI_max，则将其设为 BMI_max
if cuts[-1] > BMI_max:
cuts[-1] = BMI_max
# 为每个分组随机选择一个检测孕周
week_assign = np.random.choice(weeks_available, G_var)
population.append({'cuts': cuts, 'weeks': week_assign})
return population
# =========================
# 4. 批量预测函数 (与原代码一致，但更清晰)
# =========================
def predict_y_batch(gest_weeks, bmis, ages, heights, ivf_types, preg_times, birth_times, gam, scaler_g, bmi_for_rank, X_columns, X_cat_columns):
"""
批量预测 Y 染色体浓度和假阴性概率。
假阴性概率公式已根据你的原代码进行保留，但请注意其经验性。
"""
gest_weeks = np.array(gest_weeks)
# 标准化检测孕周
g_it_std = scaler_g.transform(gest_weeks.reshape(-1, 1)).flatten()
bmis = np.array(bmis)
# 计算 BMI 的分位数排名，结合了训练集和当前预测集的数据
combined_bmi = pd.concat([bmi_for_rank, pd.Series(bmis)])
b_i_norm_all = combined_bmi.rank(pct=True).iloc[-len(bmis):].values
# 构建预测用的 DataFrame
X_pred = pd.DataFrame({
"g_it_std": g_it_std,
"b_i_norm": b_i_norm_all,
"年龄": ages,
"身高": heights
})
# 处理类别变量，并确保列顺序与训练时一致
cat_df = pd.DataFrame({
"IVF妊娠": ivf_types,
"怀孕次数": preg_times,
"生产次数": birth_times
})
X_cat_pred = pd.get_dummies(cat_df, drop_first=True)
# 确保列存在，不存在的列填充0，保持列结构一致
X_cat_pred = X_cat_pred.reindex(columns=X_cat_columns, fill_value=0)
# 合并数值型和类别型特征
X_pred = pd.concat([X_pred, X_cat_pred], axis=1)
# 确保特征列顺序与训练时GAM模型期望的一致
X_pred = X_pred[X_columns]
# 使用 GAM 模型进行预测
y_preds = gam.predict(X_pred)
# ==========================
# 添加高斯噪声
# ==========================
noise = np.random.normal(loc=0, scale=0.0000005, size=y_preds.shape)
y_preds_noisy = y_preds + noise
# ==========================
# 计算假阴性概率 (fn_probs)
fn_probs = np.ones_like(y_preds_noisy)
mask = y_preds_noisy > 0.04 # 只对预测 Y > 0.04 的情况计算指数衰减
# 指数衰减公式，预测 Y 越高，假阴性概率越低
fn_probs[mask] = np.exp(-50 * (y_preds_noisy[mask] - 0.04))
return y_preds_noisy, fn_probs
# =========================
# 5. 适应度函数
# =========================
def risk_score(weeks, bmis):
"""
将孕周分为三档，并根据BMI进行额外惩罚，以鼓励更早的孕周。
"""
weeks = np.array(weeks)
bmis = np.array(bmis)
score = np.zeros_like(weeks, dtype=float)
# 基础分
mask_early = weeks <= 12
score[mask_early] = 1 + 0.05 * weeks[mask_early]
mask_mid = (weeks >= 13) & (weeks <= 27)
score[mask_mid] = 2 + 0.02 * (weeks[mask_mid] - 13)
mask_late = weeks >= 28
score[mask_late] = 3 + 0.01 * (weeks[mask_late] - 28)
# 针对高BMI的额外惩罚（阈值可调）
mask_high_bmi = bmis > 35
score[mask_high_bmi] += (bmis[mask_high_bmi] - 35) * 0.1 # BMI越高，额外惩罚越大
return score
def fitness(cuts, group_week_assign, df, gam, scaler_g, bmi_for_rank, X_columns, X_cat_columns):
"""
计算个体（一组 BMI 分割点和孕周）的适应度。
目标是最小化此值。
"""
# 提取当前个体需要的所有孕妇数据
bmis = df['孕妇BMI'].values
ages = df['年龄'].values
heights = df['身高'].values
ivf_types = df['IVF妊娠'].values
preg_times = df['怀孕次数'].values
birth_times = df['生产次数'].values
# 根据个体的 BMI 分割点，为所有孕妇分配组别
groups = assign_group(bmis, cuts)
# 根据分配的组别，查找对应的检测孕周
check_weeks = np.array([group_week_assign[g] for g in groups])
# 使用批量预测函数计算假阴性概率
# 注意：这里使用的是添加了噪声的预测值
_, fn_probs = predict_y_batch(
gest_weeks=check_weeks,
bmis=bmis,
ages=ages,
heights=heights,
ivf_types=ivf_types,
preg_times=preg_times,
birth_times=birth_times,
gam=gam,
scaler_g=scaler_g,
bmi_for_rank=bmi_for_rank,
X_columns=X_columns,
X_cat_columns=X_cat_columns
)
# 计算总的假阴性风险（假阴性概率的总和）
total_risk = np.sum(fn_probs)
# ========================
# 检测周评分，使用上一题的 risk_score 函数（区间递增）
# ========================
# 调用修改后的 risk_score 函数，传入 bmis
week_score = risk_score(check_weeks, bmis)
total_score = np.sum(week_score)
# BMI 区间长度的惩罚项
cut_diffs = np.diff(np.concatenate(([BMI_min], cuts, [BMI_max])))
small_interval_penalty = np.sum(cut_diffs < 3) * 100
# 当前个体的分组数量
G_var = len(cuts) + 1
# GC 含量惩罚
gc_scores = np.array([gc_score(gc) for gc in df['GC含量']])
mean_gc_score = np.mean(gc_scores)
gc_penalty = 2 * (1 - mean_gc_score)
# 最终目标函数
obj = 10 * total_risk + 2 * total_score + GROUP_PENALTY_WEIGHT * G_var + small_interval_penalty + gc_penalty
return obj
# =========================
# 6. 变异操作
# =========================
def mutate(indiv):
"""
优化：变异操作增加随机性，不仅改变值，还可以改变组数。
"""
cuts = indiv['cuts'].copy()
weeks = indiv['weeks'].copy()
# 变异类型：
# 1. 变异 BMI 分割点的值
if np.random.rand() < MUT_RATE:
# 随机调整每个分割点的值，幅度在 -5 到 5 之间
cuts += np.random.uniform(-5, 5, len(cuts))
# 裁剪分割点到 BMI_min 和 BMI_max 范围内
cuts = np.clip(cuts, BMI_min, BMI_max)
# 重新排序分割点
cuts = np.sort(cuts)
# 修正最小间距（与初始化时类似）
for i in range(1, len(cuts)):
if cuts[i] - cuts[i-1] < 3:
cuts[i] = cuts[i-1] + 3
# 移除重复的分割点，确保唯一性
cuts = np.unique(cuts)
# 2. 变异孕周分配
if np.random.rand() < MUT_RATE:
for i in range(len(weeks)):
# 为每个分组随机选择一个新的检测孕周
weeks[i] = np.random.choice(weeks_available)
# 3. 变异组数（动态调整 G_var）
# 增加一个分组（添加一个 cut）
if np.random.rand() < MUT_RATE and len(cuts) + 1 < MAX_G:
# 在现有分割点之间随机选择一个位置插入新的分割点
new_cut = np.random.uniform(cuts.min(), cuts.max())
cuts = np.sort(np.append(cuts, new_cut))
# 减少一个分组（移除一个 cut）
if np.random.rand() < MUT_RATE and len(cuts) + 1 > MIN_G:
# 随机移除一个现有的分割点
cuts = np.delete(cuts, np.random.randint(len(cuts)))
# 确保 weeks 数组的长度与 cuts 数组的长度（加1）匹配
weeks_len = len(cuts) + 1
if len(weeks) != weeks_len:
# 如果长度不匹配，则重新随机生成 weeks
weeks = np.random.choice(weeks_available, weeks_len)
return {'cuts': cuts, 'weeks': weeks}
# =========================
# 7. 交叉操作
# =========================
def crossover(parent1, parent2):
"""
优化：允许不同组数的个体进行交叉。
"""
cuts1, weeks1 = parent1['cuts'], parent1['weeks']
cuts2, weeks2 = parent2['cuts'], parent2['weeks']
# 将两个父代的 cuts 和 weeks 合并，然后重新排序和选择
# 合并两个父代的分割点，并去重，得到一个更全面的分割点集合
all_cuts = np.unique(np.concatenate([cuts1, cuts2]))
# 随机确定子代的组数，允许在 MIN_G-1 到 MAX_G-1 之间随机选择 cut 的数量
# 注意：这里的 new_cuts_size 是 cut 的数量，实际组数是 new_cuts_size + 1
new_cuts_size = np.random.randint(MIN_G - 1, MAX_G)
# 从合并后的所有分割点中，随机选择 new_cuts_size 个作为子代的分割点
if len(all_cuts) > new_cuts_size:
new_cuts = np.sort(np.random.choice(all_cuts, new_cuts_size, replace=False))
else:
# 如果所有分割点少于或等于目标数量，则全部保留
new_cuts = all_cuts
# 为子代随机分配检测孕周，数量为 new_cuts_size + 1
new_weeks = np.random.choice(weeks_available, len(new_cuts) + 1)
return {'cuts': new_cuts, 'weeks': new_weeks}
# =========================
# 8. GA 主循环
# =========================
population = init_population()
best_obj = float('inf')
best_solution = None
for gen in range(N_GEN):
# 适应度计算：为种群中的每个个体计算适应度值
fitness_values = [
fitness(indiv['cuts'], indiv['weeks'], df, gam, scaler_g, bmi_for_rank, X_columns, X_cat_columns)
for indiv in population
]
# 更新最优解：记录当前最优的个体及其适应度值
idx = np.argmin(fitness_values)
if fitness_values[idx] < best_obj:
best_obj = fitness_values[idx]
best_solution = population[idx]
# 选择：使用轮盘赌选择法，适应度值越高的个体（目标函数值越小）被选中的概率越大
# 转换为概率，适应度值越小，权重越大
weights = np.exp(-np.array(fitness_values))
weights /= weights.sum() # 归一化权重
new_population = []
# 精英策略：将当前最优解直接加入下一代种群，保证最优解不丢失
new_population.append(best_solution)
# 繁殖下一代种群
for _ in range(POP_SIZE - 1): # 减去精英个体数量
# 交叉操作
if np.random.rand() < CROSS_RATE:
# 从当前种群中按权重随机选择两个父代
parent1 = random.choices(population, weights=weights, k=1)[0]
parent2 = random.choices(population, weights=weights, k=1)[0]
child = crossover(parent1, parent2)
else:
# 如果不进行交叉，则随机选择一个父代（相当于拷贝）
parent = random.choices(population, weights=weights, k=1)[0]
child = parent
# 变异操作
child = mutate(child)
new_population.append(child)
# 更新种群
population = new_population
# 每隔10代打印一次当前最优的适应度值
if gen % 10 == 0:
print(f"Generation {gen}: best objective = {best_obj:.4f}")
# =========================
# 9. 输出最优分组结果
# =========================
best_cuts, best_weeks = best_solution['cuts'], best_solution['weeks']
print("\n最优 BMI 分组 cut:", best_cuts)
print("最优每组检测周:", best_weeks)
# 重新计算并打印详细结果，用于验证
bmis = df['孕妇BMI'].values
# 为所有孕妇根据最优分割点分配组别
groups = assign_group(bmis, best_cuts)
# 获取每个孕妇对应的最优检测周
check_weeks = np.array([best_weeks[g] for g in groups])
# 使用最优分组方案进行预测
# 注意：这里的 predict_y_batch 函数已经包含了噪声
y_preds_noisy, fn_probs = predict_y_batch(
check_weeks, bmis, df['年龄'].values, df['身高'].values, df['IVF妊娠'].values,
df['怀孕次数'].values, df['生产次数'].values,
gam, scaler_g, bmi_for_rank, X_columns, X_cat_columns
)
# 打印每个孕妇的详细信息
for code, bmi, week, y_pred, fn_prob in zip(df['孕妇代码'], bmis, check_weeks, y_preds_noisy, fn_probs):
print(f"孕妇 {code}: BMI={bmi:.2f}, 分组检测周={week}, 预测 Y={y_pred:.6f}, 假阴性={fn_prob:.4f}")
# 辅助函数，将分割点转换为区间表示
def cuts_to_intervals(cuts, BMI_min, BMI_max):
"""
将 BMI 分割点转换为可读的区间表示。
"""
intervals = []
prev = BMI_min
for c in cuts:
intervals.append((prev, c))
prev = c
intervals.append((prev, BMI_max)) # 最后一组区间
return intervals
# 打印最终确定的 BMI 分组区间和对应的检测周
best_intervals = cuts_to_intervals(best_cuts, BMI_min, BMI_max)
for i, (low, high) in enumerate(best_intervals):
print(f"组{i+1}: BMI {low:.1f}-{high:.1f}, 检测周: {best_weeks[i]}")
Generation 0: best objective = 410.4333 Generation 10: best objective = 409.6975 Generation 20: best objective = 409.4333 Generation 30: best objective = 409.4329 Generation 40: best objective = 409.4326 最优 BMI 分组 cut: [35.37965574 40.75138861 43.75138861] 最优每组检测周: [14 14 20 18] 孕妇 A003: BMI=30.74, 分组检测周=14, 预测 Y=0.064980, 假阴性=0.2868 孕妇 A008: BMI=29.55, 分组检测周=14, 预测 Y=0.080114, 假阴性=0.1346 孕妇 A010: BMI=35.56, 分组检测周=14, 预测 Y=0.086614, 假阴性=0.0972 孕妇 A012: BMI=30.80, 分组检测周=14, 预测 Y=0.059882, 假阴性=0.3701 孕妇 A018: BMI=28.13, 分组检测周=14, 预测 Y=0.067807, 假阴性=0.2490 孕妇 A023: BMI=30.04, 分组检测周=14, 预测 Y=0.071954, 假阴性=0.2024 孕妇 A026: BMI=28.04, 分组检测周=14, 预测 Y=0.059429, 假阴性=0.3785 孕妇 A029: BMI=33.11, 分组检测周=14, 预测 Y=0.081796, 假阴性=0.1237 孕妇 A041: BMI=36.81, 分组检测周=14, 预测 Y=0.063739, 假阴性=0.3052 孕妇 A048: BMI=35.23, 分组检测周=14, 预测 Y=0.067729, 假阴性=0.2500 孕妇 A055: BMI=35.00, 分组检测周=14, 预测 Y=0.057825, 假阴性=0.4101 孕妇 A061: BMI=36.79, 分组检测周=14, 预测 Y=0.074818, 假阴性=0.1754 孕妇 A065: BMI=32.05, 分组检测周=14, 预测 Y=0.063831, 假阴性=0.3038 孕妇 A066: BMI=31.22, 分组检测周=14, 预测 Y=0.081451, 假阴性=0.1259 孕妇 A069: BMI=34.18, 分组检测周=14, 预测 Y=0.080705, 假阴性=0.1306 孕妇 A071: BMI=34.72, 分组检测周=14, 预测 Y=0.057060, 假阴性=0.4261 孕妇 A073: BMI=34.55, 分组检测周=14, 预测 Y=0.071249, 假阴性=0.2096 孕妇 A092: BMI=30.48, 分组检测周=14, 预测 Y=0.051266, 假阴性=0.5693 孕妇 A095: BMI=39.14, 分组检测周=14, 预测 Y=0.061610, 假阴性=0.3394 孕妇 A096: BMI=27.90, 分组检测周=14, 预测 Y=0.072167, 假阴性=0.2002 孕妇 A105: BMI=32.41, 分组检测周=14, 预测 Y=0.058695, 假阴性=0.3927 孕妇 A109: BMI=29.41, 分组检测周=14, 预测 Y=0.062600, 假阴性=0.3230 孕妇 A130: BMI=30.10, 分组检测周=14, 预测 Y=0.067582, 假阴性=0.2518 孕妇 A133: BMI=28.27, 分组检测周=14, 预测 Y=0.049994, 假阴性=0.6067 孕妇 A146: BMI=31.64, 分组检测周=14, 预测 Y=0.078002, 假阴性=0.1496 孕妇 A147: BMI=31.41, 分组检测周=14, 预测 Y=0.076239, 假阴性=0.1633 孕妇 A155: BMI=30.47, 分组检测周=14, 预测 Y=0.066556, 假阴性=0.2651 孕妇 A163: BMI=33.31, 分组检测周=14, 预测 Y=0.056805, 假阴性=0.4316 孕妇 A169: BMI=31.35, 分组检测周=14, 预测 Y=0.080173, 假阴性=0.1342 孕妇 A170: BMI=29.44, 分组检测周=14, 预测 Y=0.064711, 假阴性=0.2907 孕妇 A171: BMI=32.70, 分组检测周=14, 预测 Y=0.073515, 假阴性=0.1872 孕妇 A172: BMI=33.74, 分组检测周=14, 预测 Y=0.087507, 假阴性=0.0930 孕妇 A175: BMI=31.83, 分组检测周=14, 预测 Y=0.073595, 假阴性=0.1864 孕妇 A176: BMI=30.50, 分组检测周=14, 预测 Y=0.080470, 假阴性=0.1322 孕妇 A178: BMI=34.57, 分组检测周=14, 预测 Y=0.061457, 假阴性=0.3420 孕妇 A192: BMI=33.29, 分组检测周=14, 预测 Y=0.081950, 假阴性=0.1228 孕妇 A194: BMI=34.65, 分组检测周=14, 预测 Y=0.066336, 假阴性=0.2680 孕妇 A196: BMI=31.36, 分组检测周=14, 预测 Y=0.081149, 假阴性=0.1278 孕妇 A198: BMI=30.65, 分组检测周=14, 预测 Y=0.075629, 假阴性=0.1684 孕妇 A200: BMI=32.09, 分组检测周=14, 预测 Y=0.077220, 假阴性=0.1555 孕妇 A201: BMI=31.26, 分组检测周=14, 预测 Y=0.072069, 假阴性=0.2012 孕妇 A202: BMI=31.68, 分组检测周=14, 预测 Y=0.074576, 假阴性=0.1775 孕妇 A203: BMI=28.17, 分组检测周=14, 预测 Y=0.073744, 假阴性=0.1850 孕妇 A204: BMI=33.73, 分组检测周=14, 预测 Y=0.076795, 假阴性=0.1589 孕妇 A205: BMI=29.58, 分组检测周=14, 预测 Y=0.070511, 假阴性=0.2175 孕妇 A215: BMI=32.72, 分组检测周=14, 预测 Y=0.076553, 假阴性=0.1608 孕妇 A216: BMI=30.28, 分组检测周=14, 预测 Y=0.070235, 假阴性=0.2205 孕妇 A219: BMI=29.01, 分组检测周=14, 预测 Y=0.079297, 假阴性=0.1402 孕妇 A224: BMI=32.19, 分组检测周=14, 预测 Y=0.083685, 假阴性=0.1126 孕妇 A227: BMI=33.01, 分组检测周=14, 预测 Y=0.070125, 假阴性=0.2217 孕妇 A230: BMI=34.21, 分组检测周=14, 预测 Y=0.073173, 假阴性=0.1904 孕妇 A233: BMI=32.17, 分组检测周=14, 预测 Y=0.073213, 假阴性=0.1900 孕妇 A234: BMI=34.75, 分组检测周=14, 预测 Y=0.074237, 假阴性=0.1805 孕妇 A235: BMI=34.27, 分组检测周=14, 预测 Y=0.078368, 假阴性=0.1468 孕妇 A236: BMI=32.76, 分组检测周=14, 预测 Y=0.068628, 假阴性=0.2390 孕妇 A242: BMI=30.03, 分组检测周=14, 预测 Y=0.074015, 假阴性=0.1826 孕妇 A245: BMI=30.76, 分组检测周=14, 预测 Y=0.080954, 假阴性=0.1290 孕妇 A247: BMI=33.21, 分组检测周=14, 预测 Y=0.082673, 假阴性=0.1184 孕妇 A248: BMI=33.10, 分组检测周=14, 预测 Y=0.072421, 假阴性=0.1977 孕妇 A250: BMI=33.57, 分组检测周=14, 预测 Y=0.076576, 假阴性=0.1606 孕妇 A255: BMI=32.14, 分组检测周=14, 预测 Y=0.068788, 假阴性=0.2371 孕妇 A256: BMI=33.34, 分组检测周=14, 预测 Y=0.078506, 假阴性=0.1458 孕妇 A263: BMI=29.35, 分组检测周=14, 预测 Y=0.074272, 假阴性=0.1802 孕妇 A266: BMI=32.97, 分组检测周=14, 预测 Y=0.072187, 假阴性=0.2000 孕妇 A267: BMI=30.70, 分组检测周=14, 预测 Y=0.079450, 假阴性=0.1391 组1: BMI 26.0-35.4, 检测周: 14 组2: BMI 35.4-40.8, 检测周: 14 组3: BMI 40.8-43.8, 检测周: 20 组4: BMI 43.8-47.0, 检测周: 18
代码版本4¶
In [8]:
import pandas as pd
import numpy as np
import random
from sklearn.preprocessing import StandardScaler
from pygam import GAM, s, l, te
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import itertools
# =========================
# 1. GA 参数
# =========================
SEED = 42
random.seed(SEED)
np.random.seed(SEED)
MAX_G = 6
MIN_G = 2
GROUP_PENALTY_WEIGHT = -0.5 # 惩罚项权重，可以调整
POP_SIZE = 200
N_GEN = 50
MUT_RATE = 0.7
CROSS_RATE = 0.7
# BMI 和孕周范围
BMI_min, BMI_max = 26, 47
weeks_available = np.arange(14, 30)
# =========================
# 2. 辅助函数
# =========================
def assign_group(bmi_array, cuts):
"""根据BMI分割点，为每个孕妇分配组别。"""
return np.searchsorted(cuts, bmi_array)
def gc_score(gc):
"""根据GC含量合理度评分。"""
if 38 <= gc <= 42:
return 1
elif 40 <= gc <= 60:
return 0.5
else:
return 0
# =========================
# 3. 优化后的初始化种群
# =========================
def init_population(pop_size, min_g, max_g):
"""初始化遗传算法种群，每个个体包含 BMI 分割点和对应的检测孕周。"""
population = []
for _ in range(pop_size):
G_var = np.random.randint(min_g, max_g + 1)
cuts = np.sort(np.random.uniform(BMI_min, BMI_max, G_var - 1))
for i in range(1, len(cuts)):
if cuts[i] - cuts[i-1] < 3:
cuts[i] = cuts[i-1] + 3
if cuts[-1] > BMI_max:
cuts[-1] = BMI_max
week_assign = np.random.choice(weeks_available, G_var)
population.append({'cuts': cuts, 'weeks': week_assign})
return population
# =========================
# 4. 批量预测函数 (与原代码一致，但更清晰)
# =========================
def predict_y_batch(gest_weeks, bmis, ages, heights, ivf_types, preg_times, birth_times, gam, scaler_g, bmi_for_rank, X_columns, X_cat_columns):
"""
批量预测 Y 染色体浓度和假阴性概率。
"""
gest_weeks = np.array(gest_weeks)
g_it_std = scaler_g.transform(gest_weeks.reshape(-1, 1)).flatten()
bmis = np.array(bmis)
combined_bmi = pd.concat([bmi_for_rank, pd.Series(bmis)])
b_i_norm_all = combined_bmi.rank(pct=True).iloc[-len(bmis):].values
X_pred = pd.DataFrame({
"g_it_std": g_it_std,
"b_i_norm": b_i_norm_all,
"年龄": ages,
"身高": heights
})
cat_df = pd.DataFrame({
"IVF妊娠": ivf_types,
"怀孕次数": preg_times,
"生产次数": birth_times
})
X_cat_pred = pd.get_dummies(cat_df, drop_first=True)
X_cat_pred = X_cat_pred.reindex(columns=X_cat_columns, fill_value=0)
X_pred = pd.concat([X_pred, X_cat_pred], axis=1)
X_pred = X_pred[X_columns]
y_preds = gam.predict(X_pred)
noise = np.random.normal(loc=0, scale=0.0005, size=y_preds.shape)
y_preds_noisy = y_preds #+ noise
fn_probs = np.ones_like(y_preds_noisy)
mask = y_preds_noisy > 0.04
fn_probs[mask] = np.exp(-50 * (y_preds_noisy[mask] - 0.04))
return y_preds_noisy, fn_probs
# =========================
# 5. 适应度函数
# =========================
def risk_score(weeks, bmis):
"""
将孕周分为三档，并根据BMI进行额外惩罚，以鼓励更早的孕周。
"""
weeks = np.array(weeks)
bmis = np.array(bmis)
score = np.zeros_like(weeks, dtype=float)
mask_early = weeks <= 12
score[mask_early] = 1 + 0.05 * weeks[mask_early]
mask_mid = (weeks >= 13) & (weeks <= 27)
score[mask_mid] = 2 + 0.02 * (weeks[mask_mid] - 13)
mask_late = weeks >= 28
score[mask_late] = 3 + 0.01 * (weeks[mask_late] - 28)
mask_high_bmi = bmis > 35
score[mask_high_bmi] += (bmis[mask_high_bmi] - 35) * 0.1
return score
def fitness(cuts, group_week_assign, df, gam, scaler_g, bmi_for_rank, X_columns, X_cat_columns, gp_weight, ts_weight):
"""
计算个体（一组 BMI 分割点和孕周）的适应度。
目标是最小化此值。
"""
bmis = df['孕妇BMI'].values
ages = df['年龄'].values
heights = df['身高'].values
ivf_types = df['IVF妊娠'].values
preg_times = df['怀孕次数'].values
birth_times = df['生产次数'].values
groups = assign_group(bmis, cuts)
check_weeks = np.array([group_week_assign[g] for g in groups])
_, fn_probs = predict_y_batch(
gest_weeks=check_weeks,
bmis=bmis,
ages=ages,
heights=heights,
ivf_types=ivf_types,
preg_times=preg_times,
birth_times=birth_times,
gam=gam,
scaler_g=scaler_g,
bmi_for_rank=bmi_for_rank,
X_columns=X_columns,
X_cat_columns=X_cat_columns
)
total_risk = np.sum(fn_probs)
week_score = risk_score(check_weeks, bmis)
total_score = np.sum(week_score)
cut_diffs = np.diff(np.concatenate(([BMI_min], cuts, [BMI_max])))
small_interval_penalty = np.sum(cut_diffs < 3) * 100
G_var = len(cuts) + 1
gc_scores = np.array([gc_score(gc) for gc in df['GC含量']])
mean_gc_score = np.mean(gc_scores)
gc_penalty = 2 * (1 - mean_gc_score)
obj = 10 * total_risk + ts_weight * total_score + gp_weight * G_var + small_interval_penalty + gc_penalty
return obj
# =========================
# 6. 变异操作
# =========================
def mutate(indiv):
"""
优化：变异操作增加随机性，不仅改变值，还可以改变组数。
"""
cuts = indiv['cuts'].copy()
weeks = indiv['weeks'].copy()
if np.random.rand() < MUT_RATE:
cuts += np.random.uniform(-5, 5, len(cuts))
cuts = np.clip(cuts, BMI_min, BMI_max)
cuts = np.sort(cuts)
for i in range(1, len(cuts)):
if cuts[i] - cuts[i-1] < 3:
cuts[i] = cuts[i-1] + 3
cuts = np.unique(cuts)
if np.random.rand() < MUT_RATE:
for i in range(len(weeks)):
weeks[i] = np.random.choice(weeks_available)
if np.random.rand() < MUT_RATE and len(cuts) + 1 < MAX_G:
new_cut = np.random.uniform(cuts.min(), cuts.max())
cuts = np.sort(np.append(cuts, new_cut))
if np.random.rand() < MUT_RATE and len(cuts) + 1 > MIN_G:
cuts = np.delete(cuts, np.random.randint(len(cuts)))
weeks_len = len(cuts) + 1
if len(weeks) != weeks_len:
weeks = np.random.choice(weeks_available, weeks_len)
return {'cuts': cuts, 'weeks': weeks}
# =========================
# 7. 交叉操作
# =========================
def crossover(parent1, parent2):
"""
优化：允许不同组数的个体进行交叉。
"""
cuts1, weeks1 = parent1['cuts'], parent1['weeks']
cuts2, weeks2 = parent2['cuts'], parent2['weeks']
all_cuts = np.unique(np.concatenate([cuts1, cuts2]))
new_cuts_size = np.random.randint(MIN_G - 1, MAX_G)
if len(all_cuts) > new_cuts_size:
new_cuts = np.sort(np.random.choice(all_cuts, new_cuts_size, replace=False))
else:
new_cuts = all_cuts
new_weeks = np.random.choice(weeks_available, len(new_cuts) + 1)
return {'cuts': new_cuts, 'weeks': new_weeks}
# =========================
# 8. GA 主循环
# =========================
def run_ga(pop_size, n_gen, gp_weight, ts_weight, **kwargs):
"""
执行遗传算法主循环并返回最优解。
"""
population = init_population(pop_size, MIN_G, MAX_G)
best_obj = float('inf')
best_solution = None
for gen in range(n_gen):
fitness_values = [
fitness(indiv['cuts'], indiv['weeks'], **kwargs, gp_weight=gp_weight, ts_weight=ts_weight)
for indiv in population
]
idx = np.argmin(fitness_values)
if fitness_values[idx] < best_obj:
best_obj = fitness_values[idx]
best_solution = population[idx]
weights = np.exp(-np.array(fitness_values))
weights /= weights.sum()
new_population = [best_solution] if best_solution else []
for _ in range(pop_size - 1):
if np.random.rand() < CROSS_RATE:
parent1 = random.choices(population, weights=weights, k=1)[0]
parent2 = random.choices(population, weights=weights, k=1)[0]
child = crossover(parent1, parent2)
else:
parent = random.choices(population, weights=weights, k=1)[0]
child = parent
child = mutate(child)
new_population.append(child)
population = new_population
if gen % 10 == 0:
print(f"Generation {gen}: best objective = {best_obj:.4f}")
return best_solution, best_obj
# =========================
# 9. 网格化调参
# =========================
param_grid = {
'GROUP_PENALTY_WEIGHT': [-0.2, -0.5, -0.8],
'TOTAL_SCORE_WEIGHT': [0.5, 1.0, 2.0]
}
param_combinations = list(itertools.product(
param_grid['GROUP_PENALTY_WEIGHT'],
param_grid['TOTAL_SCORE_WEIGHT']
))
N_GEN_GRID = 30
POP_SIZE_GRID = 100
print("\n--- 开始网格化调参 ---")
for gp_weight, ts_weight in param_combinations:
print(f"\n正在尝试参数组合: GROUP_PENALTY_WEIGHT={gp_weight}, TOTAL_SCORE_WEIGHT={ts_weight}")
best_solution, best_obj = run_ga(
pop_size=POP_SIZE_GRID,
n_gen=N_GEN_GRID,
gp_weight=gp_weight,
ts_weight=ts_weight,
df=df,
gam=gam,
scaler_g=scaler_g,
bmi_for_rank=bmi_for_rank,
X_columns=X_columns,
X_cat_columns=X_cat_columns
)
if best_solution:
num_groups = len(best_solution['cuts']) + 1
is_monotonic = np.all(np.diff(best_solution['weeks']) >= 0)
is_early_week = np.all(best_solution['weeks'] < 21)
if 3 <= num_groups <= 5 and is_monotonic and is_early_week:
print("\n🎉 找到一个符合条件的解！")
print(f"参数: GROUP_PENALTY_WEIGHT={gp_weight}, TOTAL_SCORE_WEIGHT={ts_weight}")
print(f"最优解: BMI 分组={best_solution['cuts']}, 推荐孕周={best_solution['weeks']}")
print(f"目标函数值: {best_obj:.4f}")
else:
print("\n该参数组合在30代内未找到符合条件的解。")
print("\n--- 网格化调参结束 ---")
--- 开始网格化调参 --- 正在尝试参数组合: GROUP_PENALTY_WEIGHT=-0.2, TOTAL_SCORE_WEIGHT=0.5 Generation 0: best objective = 212.2737 Generation 10: best objective = 212.1625 Generation 20: best objective = 212.0298 该参数组合在30代内未找到符合条件的解。 正在尝试参数组合: GROUP_PENALTY_WEIGHT=-0.2, TOTAL_SCORE_WEIGHT=1.0 Generation 0: best objective = 278.7910 Generation 10: best objective = 278.6810 Generation 20: best objective = 278.4810 该参数组合在30代内未找到符合条件的解。 正在尝试参数组合: GROUP_PENALTY_WEIGHT=-0.2, TOTAL_SCORE_WEIGHT=2.0 Generation 0: best objective = 411.0337 Generation 10: best objective = 411.0337 Generation 20: best objective = 410.9918 🎉 找到一个符合条件的解！ 参数: GROUP_PENALTY_WEIGHT=-0.2, TOTAL_SCORE_WEIGHT=2.0 最优解: BMI 分组=[39.53081909 43.76147835], 推荐孕周=[14 15 19] 目标函数值: 410.8337 正在尝试参数组合: GROUP_PENALTY_WEIGHT=-0.5, TOTAL_SCORE_WEIGHT=0.5 Generation 0: best objective = 210.2540 Generation 10: best objective = 210.0002 Generation 20: best objective = 210.0002 该参数组合在30代内未找到符合条件的解。 正在尝试参数组合: GROUP_PENALTY_WEIGHT=-0.5, TOTAL_SCORE_WEIGHT=1.0 Generation 0: best objective = 278.2861 Generation 10: best objective = 277.7810 Generation 20: best objective = 277.7810 该参数组合在30代内未找到符合条件的解。 正在尝试参数组合: GROUP_PENALTY_WEIGHT=-0.5, TOTAL_SCORE_WEIGHT=2.0 Generation 0: best objective = 410.4337 Generation 10: best objective = 409.9337 Generation 20: best objective = 409.9337 该参数组合在30代内未找到符合条件的解。 正在尝试参数组合: GROUP_PENALTY_WEIGHT=-0.8, TOTAL_SCORE_WEIGHT=0.5 Generation 0: best objective = 208.1489 Generation 10: best objective = 208.0642 Generation 20: best objective = 208.0642 该参数组合在30代内未找到符合条件的解。 正在尝试参数组合: GROUP_PENALTY_WEIGHT=-0.8, TOTAL_SCORE_WEIGHT=1.0 Generation 0: best objective = 278.1566 Generation 10: best objective = 276.0953 Generation 20: best objective = 276.0810 该参数组合在30代内未找到符合条件的解。 正在尝试参数组合: GROUP_PENALTY_WEIGHT=-0.8, TOTAL_SCORE_WEIGHT=2.0 Generation 0: best objective = 409.8337 Generation 10: best objective = 409.0337 Generation 20: best objective = 409.0337 该参数组合在30代内未找到符合条件的解。 --- 网格化调参结束 ---