# -*- coding: utf-8 -*- """ Created on Wed Nov 26 21:07:08 2025 @author: ALY """ # -*- coding: utf-8 -*- ''' #顶部注释Format: File-》Setting-》File and Code Templates-》Python Script @project: Python_CodesSample @author: Jason.Fan @datetime: 2025/10/13 9:52 @IDE: PyCharm @file: demo3.py ''' import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns from sklearn.model_selection import train_test_split, cross_val_score from sklearn.preprocessing import StandardScaler, LabelEncoder from sklearn.impute import SimpleImputer from sklearn.linear_model import LogisticRegression from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier from sklearn.tree import DecisionTreeClassifier, plot_tree from sklearn.naive_bayes import GaussianNB from sklearn.neural_network import MLPClassifier #用于分类任务的多层感知器（Multilayer Perceptron）模型，也被称为人工神经网络（ANN） from sklearn.metrics import (accuracy_score, precision_score, recall_score, f1_score, roc_auc_score, confusion_matrix, classification_report, roc_curve) from sklearn.preprocessing import LabelEncoder import warnings warnings.filterwarnings('ignore') # 设置中文字体 plt.rcParams['font.sans-serif'] = ['SimHei'] plt.rcParams['axes.unicode_minus'] = False def create_wealth_adequacy_target(df): """ 根据文档定义创建二分类目标变量 """ def calculate_annual_consumption(df): """ 计算年消费需求 - 加总各项消费支出 """ # 基础生活消费 basic_consumption = df['clothing_food_housing_etc'] # 品质消费（教育、医疗、通讯等） quality_consumption = (df['education_fee'] + df['medic_care'] + df['communi_fees'] + df['online_store']) # 住房相关消费 housing_consumption = df['hous_mainte'] # 绿色消费（可持续交通、环保家居） green_consumption = df['susta_transp'] + df['env_fri_hom_furn'] # 年消费总额 annual_consumption = (basic_consumption + quality_consumption + housing_consumption + green_consumption) return annual_consumption # 计算所需退休财富（基于文档公式） def calculate_required_wealth(row): annual_consumption = row['annual_consumption'] # 年消费需求 current_age = row['age'] retirement_age = 60 max_age = 99 discount_rate = 0.025 years_remaining = max_age - max(current_age, retirement_age) if years_remaining <= 0: return 0 # 年金现值公式 required_wealth = annual_consumption * ( (1 - (1 + discount_rate) ** -years_remaining) / discount_rate ) return required_wealth # 计算总财富（基于心理账户理论） def calculate_total_wealth(row): # 当前收入账户 current_income = (row.get('wage', 0) + row.get('oper_inc', 0) + row.get('prop_inc', 0) + row.get('transf_inc', 0) + row.get('other_inc', 0)) # 当前资产账户 current_asset = (row.get('deposits_stocks_bonds_funds_etc', 0) + row.get('real_esta_hous_prop_etc', 0) + row.get('car', 0) + row.get('other_nonfin_ass', 0) + row.get('digital_asset', 0)) # 未来收入账户 future_income = (row.get('hou_pro_fund1', 0) + row.get('soci_sec', 0) + row.get('corp_hous_fund', 0) + row.get('medi_ins', 0)) return current_income + current_asset + future_income # 应用计算 df['annual_consumption'] = calculate_annual_consumption(df) df['required_wealth'] = df.apply(calculate_required_wealth, axis=1) df['total_wealth'] = df.apply(calculate_total_wealth, axis=1) # 创建二分类目标变量 df['wealth_adequacy'] = (df['total_wealth'] >= df['required_wealth']).astype(int) # 输出统计信息 print(f"总样本数: {len(df)}") print(f"财富充足样本: {df['wealth_adequacy'].sum()}") print(f"财富不足样本: {len(df) - df['wealth_adequacy'].sum()}") print(f"充足率: {df['wealth_adequacy'].mean():.2%}") print(f"平均所需财富: {df['required_wealth'].mean():.2f}") print(f"平均实际财富: {df['total_wealth'].mean():.2f}") # return df['wealth_adequacy'] return df class PensionWealthPredictor: def __init__(self): self.models = {} self.scaler = StandardScaler() self.imputer = SimpleImputer(strategy='median') self.results = {} self.feature_importance = None def generate_synthetic_data(self, n_samples=5000): """ 生成合成数据模拟CHFS调查数据 基于文档中描述的特征变量 """ np.random.seed(42) data = { # 当前收入相关特征 'wage_income': np.random.lognormal(10, 0.8, n_samples), 'operating_income': np.random.lognormal(8, 1.2, n_samples), 'property_income': np.random.lognormal(7, 1.5, n_samples), 'transfer_income': np.random.lognormal(9, 0.5, n_samples), 'other_income': np.random.lognormal(6, 1.0, n_samples), # 当前资产相关特征 'land_assets': np.random.lognormal(11, 1.0, n_samples), 'car_assets': np.random.lognormal(9, 0.8, n_samples), 'other_properties': np.random.lognormal(10, 1.2, n_samples), 'dividend_income': np.random.lognormal(7, 1.3, n_samples), 'bank_savings': np.random.lognormal(12, 0.7, n_samples), 'digital_assets': np.random.lognormal(6, 2.0, n_samples), # 未来收入相关特征 'home_equity': np.random.lognormal(13, 0.9, n_samples), 'provident_fund': np.random.lognormal(11, 0.6, n_samples), 'insurance_assets': np.random.lognormal(9, 0.8, n_samples), 'private_pension': np.random.lognormal(10, 0.7, n_samples), # 支出相关特征 'necessity_expenditure': np.random.lognormal(10, 0.5, n_samples), 'quality_expenditure': np.random.lognormal(8, 0.9, n_samples), 'green_expenditure': np.random.lognormal(7, 1.1, n_samples), # 人口统计特征 'age': np.random.randint(40, 85, n_samples), 'gender': np.random.choice([0, 1], n_samples, p=[0.48, 0.52]), 'education': np.random.choice([0, 1, 2, 3], n_samples, p=[0.1, 0.3, 0.4, 0.2]), 'marital_status': np.random.choice([0, 1], n_samples, p=[0.2, 0.8]), 'employment_status': np.random.choice([0, 1, 2], n_samples, p=[0.3, 0.5, 0.2]), 'household_size': np.random.randint(1, 6, n_samples), 'health_status': np.random.choice([1, 2, 3, 4, 5], n_samples, p=[0.1, 0.2, 0.4, 0.2, 0.1]) } df = pd.DataFrame(data) # 计算总财富（基于行为生命周期理论的心理账户） df['current_income_total'] = df[['wage_income', 'operating_income', 'property_income', 'transfer_income', 'other_income']].sum(axis=1) df['current_asset_total'] = df[['land_assets', 'car_assets', 'other_properties', 'dividend_income', 'bank_savings', 'digital_assets']].sum(axis=1) df['future_income_total'] = df[['home_equity', 'provident_fund', 'insurance_assets', 'private_pension']].sum(axis=1) df['total_expenditure'] = df[['necessity_expenditure', 'quality_expenditure', 'green_expenditure']].sum(axis=1) # 计算总财富 df['total_wealth'] = df['current_income_total'] + df['current_asset_total'] + df['future_income_total'] # 生成目标变量：财富充足性（基于文档中的计算方法） # 简化的计算：如果总财富超过年龄调整的阈值，则为充足 age_adjusted_threshold = df['age'] * 10000 + df['education'] * 50000 + df['employment_status'] * 30000 df['wealth_adequacy'] = (df['total_wealth'] > age_adjusted_threshold).astype(int) # 添加一些噪声使数据更真实 noise = np.random.normal(0, 0.1, n_samples) df['wealth_adequacy'] = (df['total_wealth'] * (1 + noise) > age_adjusted_threshold).astype(int) return df def load_data(self): df = pd.read_stata('5ys.dta') return create_wealth_adequacy_target(df) def preprocess_data(self, df): """数据预处理""" # 分离特征和目标变量 # categorical_cols = df.select_dtypes(include=['object', 'category']).columns # df = pd.get_dummies(df, columns=categorical_cols, drop_first=True) # print(df['province'][:10]) non_numeric_cols = df.select_dtypes(exclude=[np.number]).columns ## 识别非数值列 le = LabelEncoder() # df[non_numeric_cols] = df[non_numeric_cols].apply(le.fit_transform) #同下功能 # 对每个非数值列进行标签编码 for col in non_numeric_cols: # 使用fit_transform进行拟合和转换 df[col] = le.fit_transform(df[col]) # # print(df['province'][:10]) # input('=====') X = df.drop('wealth_adequacy', axis=1) y = df['wealth_adequacy'] # 处理缺失值 X_imputed = self.imputer.fit_transform(X) X_imputed = pd.DataFrame(X_imputed, columns=X.columns) # 标准化数值特征 numerical_cols = X.select_dtypes(include=[np.number]).columns X_imputed[numerical_cols] = self.scaler.fit_transform(X_imputed[numerical_cols]) return X_imputed, y def initialize_models(self): """初始化所有机器学习模型""" self.models = { 'Logistic Regression': LogisticRegression(random_state=42, max_iter=1000), 'Decision Tree': DecisionTreeClassifier(random_state=42, max_depth=10), 'Random Forest': RandomForestClassifier(random_state=42, n_estimators=100), 'Gradient Boosting': GradientBoostingClassifier(random_state=42, n_estimators=100), 'Naive Bayes': GaussianNB(), 'Neural Network': MLPClassifier(random_state=42, hidden_layer_sizes=(100, 50), max_iter=1000) } def train_and_evaluate(self, X_train, X_test, y_train, y_test): """训练和评估所有模型""" results = {} for name, model in self.models.items(): print(f"训练 {name}...") # 训练模型 model.fit(X_train, y_train) # 预测 y_pred = model.predict(X_test) y_pred_proba = model.predict_proba(X_test)[:, 1] if hasattr(model, 'predict_proba') else model.decision_function( X_test) # 计算评估指标 results[name] = { 'accuracy': accuracy_score(y_test, y_pred), 'precision': precision_score(y_test, y_pred), 'recall': recall_score(y_test, y_pred), 'f1_score': f1_score(y_test, y_pred), 'roc_auc': roc_auc_score(y_test, y_pred_proba) if len(np.unique(y_test)) > 1 else 0.5, 'model': model } if name == 'Decision Tree': # 绘制决策树 plt.figure(figsize=(20, 12)) # plot_tree(model, feature_names=['annual_consumption','required_wealth', 'total_wealth'], class_names=['0', '1'], filled=True,max_depth=3) plot_tree(model, feature_names=X_train.columns.tolist(), class_names=['0', '1'], filled=True, fontsize=12, max_depth=5, rounded=True, proportion=True) plt.savefig('tree.png', dpi =300, bbox_inches = 'tight', facecolor = 'white', transparent = False, pad_inches = 0.1, orientation = 'landscape' ) plt.show() print(f"{name} 完成 - 准确率: {results[name]['accuracy']:.4f}") self.results = results return results def plot_feature_importance(self, feature_names, top_n=15): # """绘制特征重要性（基于随机森林）""" # if 'Random Forest' in self.models: # rf_model = self.models['Random Forest'] # importance = rf_model.feature_importances_ # # # 创建特征重要性DataFrame # feature_imp_df = pd.DataFrame({ # 'feature': feature_names, # 'importance': importance # }).sort_values('importance', ascending=False) # # self.feature_importance = feature_imp_df # # # 绘制图表 # plt.figure(figsize=(12, 8)) # sns.barplot(data=feature_imp_df.head(top_n), x='importance', y='feature') # plt.title('特征重要性排名 - 随机森林模型') # plt.xlabel('重要性分数') # plt.tight_layout() # plt.show() # # return feature_imp_df if 'Decision Tree' in self.models: rf_model = self.models['Decision Tree'] importance = rf_model.feature_importances_ # 创建特征重要性DataFrame feature_imp_df = pd.DataFrame({ 'feature': feature_names, 'importance': importance }).sort_values('importance', ascending=False) self.feature_importance = feature_imp_df # 绘制图表 plt.figure(figsize=(12, 8)) sns.barplot(data=feature_imp_df.head(top_n), x='importance', y='feature') plt.title('特征重要性排名 - 决策树模型') plt.xlabel('Importance score') plt.tight_layout() plt.show() return feature_imp_df def plot_model_comparison(self): """比较所有模型的性能""" if not self.results: print("请先训练模型") return metrics = ['accuracy', 'precision', 'recall', 'f1_score', 'roc_auc'] fig, axes = plt.subplots(2, 3, figsize=(18, 12)) axes = axes.ravel() for i, metric in enumerate(metrics): model_names = list(self.results.keys()) scores = [self.results[name][metric] for name in model_names] bars = axes[i].bar(model_names, scores, color=sns.color_palette('viridis', len(model_names))) axes[i].set_title(f'{metric.upper()}') axes[i].set_ylabel(metric) axes[i].tick_params(axis='x', rotation=45) # 在柱子上显示数值 for bar, score in zip(bars, scores): height = bar.get_height() axes[i].text(bar.get_x() + bar.get_width() / 2., height, f'{score:.3f}', ha='center', va='bottom') # 移除多余的子图 for i in range(len(metrics), len(axes)): fig.delaxes(axes[i]) plt.tight_layout() plt.show() def plot_roc_curves(self, X_test, y_test): """绘制所有模型的ROC曲线""" plt.figure(figsize=(10, 8)) for name, result in self.results.items(): if 'model' in result: model = result['model'] if hasattr(model, 'predict_proba'): y_pred_proba = model.predict_proba(X_test)[:, 1] fpr, tpr, _ = roc_curve(y_test, y_pred_proba) auc_score = roc_auc_score(y_test, y_pred_proba) plt.plot(fpr, tpr, label=f'{name} (AUC = {auc_score:.3f})') plt.plot([0, 1], [0, 1], 'k--', label='Random Classifier') plt.xlabel('False positive rate') plt.ylabel('True positive rate') plt.title('ROC曲线比较') plt.legend() plt.grid(True) plt.show() def get_best_model(self): """获取最佳模型（基于F1分数）""" if not self.results: print("请先训练模型") return None best_model_name = max(self.results.keys(), key=lambda x: self.results[x]['f1_score']) best_model = self.results[best_model_name]['model'] print(f"最佳模型: {best_model_name}") print(f"F1分数: {self.results[best_model_name]['f1_score']:.4f}") print(f"准确率: {self.results[best_model_name]['accuracy']:.4f}") return best_model_name, best_model def detailed_model_report(self): """生成详细的模型评估报告""" if not self.results: print("请先训练模型") return report_data = [] for name, metrics in self.results.items(): report_data.append({ 'Model': name, 'Accuracy': f"{metrics['accuracy']:.4f}", 'Precision': f"{metrics['precision']:.4f}", 'Recall': f"{metrics['recall']:.4f}", 'F1-Score': f"{metrics['f1_score']:.4f}", 'AUC-ROC': f"{metrics['roc_auc']:.4f}" }) report_df = pd.DataFrame(report_data) print("模型性能详细报告:") print("=" * 80) print(report_df.to_string(index=False)) return report_df # 主执行函数 def main(): print("养老财富充足性预测模型") print("=" * 50) # 初始化预测器 predictor = PensionWealthPredictor() # 1. 生成模拟数据 # print("1. 生成模拟数据...") # df = predictor.generate_synthetic_data(n_samples=8000) # 1.1 读取原始数据 df = predictor.load_data() print(f"数据形状: {df.shape}") print(f"目标变量分布:\n{df['wealth_adequacy'].value_counts()}") print(f"财富充足比例: {df['wealth_adequacy'].mean():.2%}") # 2. 数据预处理 print("\n2. 数据预处理...") X, y = predictor.preprocess_data(df) # 3. 划分训练集和测试集 X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4, random_state=42, stratify=y) print(f"训练集大小: {X_train.shape}") print(f"测试集大小: {X_test.shape}") # 4. 初始化模型 print("\n3. 初始化机器学习模型...") predictor.initialize_models() # 5. 训练和评估模型 print("\n4. 训练和评估模型...") results = predictor.train_and_evaluate(X_train, X_test, y_train, y_test) # 6. 生成报告和可视化 print("\n5. 生成分析报告...") # 详细报告 report_df = predictor.detailed_model_report() # 特征重要性 print("\n6. 特征重要性分析...") feature_imp_df = predictor.plot_feature_importance(X.columns) print("前10个重要特征:") print(feature_imp_df.head(10)) # 模型比较 print("\n7. 模型性能比较...") predictor.plot_model_comparison() # ROC曲线 print("\n8. ROC曲线分析...") predictor.plot_roc_curves(X_test, y_test) # 最佳模型 print("\n9. 最佳模型识别...") best_model_name, best_model = predictor.get_best_model() # 10. 政策建议（基于特征重要性） print("\n10. 基于分析结果的政策建议:") print("=" * 50) top_features = feature_imp_df.head(5)['feature'].tolist() print("影响养老财富充足性的关键因素:") for i, feature in enumerate(top_features, 1): print(f"{i}. {feature}") # 影响养老财富充足性的关键因素: # 1. # required_wealth # 2. # annual_consumption # 3. # total_wealth # 4. # real_esta_hous_prop_etc ## 房地产资产 # 5. # Hous_net_worth ## 家庭净资产 if __name__ == "__main__": main()