3. 北京市住房价格预测

# 导入必要的库
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
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.preprocessing import StandardScaler

import warnings
warnings.filterwarnings('ignore')

# 设置中文字体（使用 macOS 自带字体）
plt.rcParams['font.sans-serif'] = ['PingFang HK', 'Songti SC', 'Arial Unicode MS', 'Heiti TC']  # 优先使用宋体-简
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题

plt.style.use('ggplot')  # 设置图表风格
pd.set_option('display.max_columns', 50)  # 显示最多50列

# 初始化 Faker 库
fake = Faker('zh_CN')  # 使用中文数据
Faker.seed(42)  # 设置随机种子确保结果可重现

print("库导入成功！")

# 通过requests使用代理下载CSV
url = 'https://gist.githubusercontent.com/ffutop/870eb32f8ef6e1cfd1fc9e40d2fa63e5/raw/1cac678610a703bde5b366f2c312a211034af9b1/challenge-1-beijing.csv'

df = pd.read_csv(url)
df.head()

print(f"数据维度: {df.shape}")
print("\n前5行数据:")
display(df.head())
print("\n数据摘要:")
display(df.info())
print("\n描述统计:")
display(df.describe(include='all'))

target_column = "每平米价格"

features = df[df.select_dtypes(include=np.number).columns.drop([target_column])]
target = df[target_column]

features.columns

pd.concat([features, target], axis=1).head()

from sklearn.model_selection import train_test_split

# 定义特征矩阵和目标变量
X = features  # 特征矩阵，包含所有数值型预测变量
y = target    # 目标变量，即"每平米价格"

# 将数据集拆分为训练集和测试集
# test_size=0.3: 30%的数据用于测试，70%用于训练
# random_state=42: 设置随机种子以确保结果可复现
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

# 打印各个数据集的形状，以验证拆分是否正确
# X_train: 训练集特征矩阵的形状 (行数,列数)
# y_train: 训练集目标变量的形状 (行数,)
# X_test: 测试集特征矩阵的形状 (行数,列数)
# y_test: 测试集目标变量的形状 (行数,)
print(X_train.shape, y_train.shape, X_test.shape, y_test.shape)

from sklearn.linear_model import LinearRegression

model = LinearRegression(fit_intercept=True)
model.fit(X_train, y_train)

print("\n参数个数: ", len(model.coef_))

# 检查每个预测模型的系数：表示每个特征变化一个单位时，目标变量（每平米价格）的变化量
for name, coef in zip(X_train.columns, model.coef_):
    print(f"{name}: {coef}")

# 检查截距项
print("\nIntercept: ", model.intercept_)

# 获取测试集的真实目标值（转为numpy数组，便于后续计算）
y_true = y_test.values

# 使用训练好的线性回归模型对测试集特征进行预测，得到预测值
y_preds = model.predict(X_test)

from sklearn import metrics

def evaluate_regression(y_true, y_preds):
    """计算并打印多个回归评估指标"""
    mse = metrics.mean_squared_error(y_true, y_preds)
    rmse = np.sqrt(mse)
    mae = metrics.mean_absolute_error(y_true, y_preds)
    mape = metrics.mean_absolute_percentage_error(y_true, y_preds)
    r2 = metrics.r2_score(y_true, y_preds)
    
    print("回归模型评估指标:")
    print(f"MSE:  {mse:,.4f}")
    print(f"RMSE: {rmse:,.4f}")
    print(f"MAE:  {mae:.4f}")
    print(f"MAPE:  {mape*100:.2f}%")
    print(f"R²:   {r2*100:.2f}%")
    
    return {"mse": mse, "rmse": rmse, "mae": mae, "r2": r2}

# 使用示例
metrics_dict = evaluate_regression(y_true, y_preds)

# MAE：衡量预测误差和绝对值大小的指标
# MSE：相比于MAE，MSE对大误差的惩罚更重，MSE会比MAE更大
# MAPE：模型的平均百分比误差，评估预测结果的相对准确性

metrics.get_scorer_names()

import statsmodels.api as sm

# 添加常数项，以拟合截距
X_train_const = sm.add_constant(X_train)

# 创建模型并拟合
model2 = sm.OLS(y_train, X_train_const)
result2 = model2.fit()
# Coef: 显示每个自变量的估计系数，即每个自变量对因变量的影响大小
print("Model parameters:", result2.params)

# 输出模型摘要
print(result2.summary())

from pycaret.regression import *
exp_setup = setup(data = df,  target = '每平米价格')
all_metrics = get_metrics()

plt.rcParams['font.sans-serif'] = ['PingFang HK', 'Songti SC', 'Arial Unicode MS', 'Heiti TC'] 
plt.rcParams['font.family'] = ['sans-serif']
plt.rcParams['axes.unicode_minus'] = False  

# 创建模型
lr_model = create_model('lr')  # 线性回归模型

# 残差图 - 检查误差分布
plot_model(lr_model, plot='residuals')

# 预测误差图 - 查看预测与实际值的对比
plot_model(lr_model, plot='error')

# 特征重要性图 - 了解哪些特征影响最大，使用PyCaret的plot_model时指定字体
plot_model(lr_model, plot='feature')

# 学习曲线 - 评估模型的泛化能力
plot_model(lr_model, plot='learning')

# 决策树可视化(如果使用决定树模型)
dt_model = create_model('dt')
plot_model(dt_model, plot='tree')

compare_models(include = ['lr', 'xgboost', 'rf', 'et', 'lightgbm'])

import shap

print("Model coefficients:\n")
for i in range(X.shape[1]-1):
    print(X.columns[i], "=", model.coef_[i].round(5))

# 提取系数
coef = model.coef_
feature_names = X_test.columns

# 创建系数DataFrame并排序
coef_df = pd.DataFrame({'特征': feature_names, '系数': coef})
coef_df = coef_df.reindex(coef_df['系数'].abs().sort_values(ascending=False).index)

# 绘制水平条形图
plt.figure(figsize=(10, 8))
plt.barh(coef_df['特征'], coef_df['系数'])
plt.xlabel('系数值')
plt.ylabel('特征')
plt.title('线性模型系数大小')
plt.axvline(x=0, color='gray', linestyle='--')
plt.tight_layout()
plt.show()

# 2. SHAP值可视化
# SHAP值提供了更精细的特征贡献分析，特别是对多样本的整体解释。

import matplotlib.pyplot as plt
import shap

plt.rcParams['font.sans-serif'] = ['PingFang HK', 'Songti SC', 'Arial Unicode MS', 'Heiti TC'] 

plt.rcParams.update({'font.size': 10}) # Set desired font size

# 创建解释器
explainer = shap.LinearExplainer(model, X_test)
shap_values = explainer.shap_values(X_test)

# 摘要图 - 展示所有特征的整体影响
shap.summary_plot(shap_values, X_test, plot_type="bar")  # 条形图版本
shap.summary_plot(shap_values, X_test)  # 点图版本

# 特征依赖图 - 分析单个特征效应
shap.dependence_plot("面积", shap_values, X_test)  # 替换"地铁"为您的特征名

# features = pd.DataFrame(features, columns=["楼层"])  # 根据实际情况修改列名

shap.partial_dependence_plot(
    "楼层",
    model.predict,
    features,
    ice=False,
    model_expected_value=True,
    feature_expected_value=True,
)

from faker import Faker
import random

# 使用 Faker 库创建更真实的北京市住房数据集
np.random.seed(42)
n_samples = 1000

# 北京各区域及其特征
beijing_districts = {
    '东城区': {'center_distance': 2, 'metro_density': 0.9, 'school_quality': 0.9},
    '西城区': {'center_distance': 3, 'metro_density': 0.8, 'school_quality': 0.95},
    '朝阳区': {'center_distance': 8, 'metro_density': 0.7, 'school_quality': 0.6},
    '海淀区': {'center_distance': 12, 'metro_density': 0.6, 'school_quality': 0.9},
    '丰台区': {'center_distance': 15, 'metro_density': 0.5, 'school_quality': 0.4},
    '石景山区': {'center_distance': 18, 'metro_density': 0.4, 'school_quality': 0.3},
    '通州区': {'center_distance': 25, 'metro_density': 0.3, 'school_quality': 0.2},
    '大兴区': {'center_distance': 30, 'metro_density': 0.2, 'school_quality': 0.1}
}

# 生成特征数据
data = {
    '区域': [fake.random_element(elements=list(beijing_districts.keys())) for _ in range(n_samples)],
    '面积': [fake.random_int(min=30, max=200) for _ in range(n_samples)],  # 房屋面积（平方米）
    '房间数': [fake.random_int(min=1, max=5) for _ in range(n_samples)],  # 房间数
    '楼层': [fake.random_int(min=1, max=30) for _ in range(n_samples)],  # 楼层
    '建造年份': [fake.random_int(min=1990, max=2023) for _ in range(n_samples)],  # 建造年份
    '装修程度': [np.random.choice([1, 2, 3, 4, 5], p=[0.1, 0.2, 0.3, 0.3, 0.1]) for _ in range(n_samples)],  # 装修程度
    '学区房': [fake.boolean(chance_of_getting_true=30) for _ in range(n_samples)],  # 是否为学区房
    '地铁站数量': [fake.random_int(min=0, max=5) for _ in range(n_samples)],  # 附近地铁站数量
    '商业设施': [fake.random_int(min=0, max=10) for _ in range(n_samples)],  # 附近商业设施数量
    '绿化率': [fake.random_int(min=10, max=80) for _ in range(n_samples)]  # 小区绿化率
}

# 根据区域计算距离市中心的距离
def calculate_center_distance(district):
    return beijing_districts[district]['center_distance'] + fake.random_int(min=0, max=5)

def calculate_metro_distance(district):
    base_distance = 5 - beijing_districts[district]['metro_density'] * 4
    return max(0.1, base_distance + fake.random_int(min=0, max=3))

# 添加计算的距离特征
data['距离市中心'] = [calculate_center_distance(district) for district in data['区域']]
data['地铁距离'] = [calculate_metro_distance(district) for district in data['区域']]

# 创建DataFrame
df = pd.DataFrame(data)

# 生成更真实的房价（基于北京房价特征）
def calculate_house_price(row):
    """基于特征计算房价"""
    base_price = 30000  # 基础价格（元/平方米）
    
    # 区域影响（不同区域基础价格不同）
    district_multiplier = {
        '东城区': 1.8, '西城区': 1.9, '朝阳区': 1.6, '海淀区': 1.7,
        '丰台区': 1.2, '石景山区': 1.1, '通州区': 0.8, '大兴区': 0.7
    }
    
    # 面积影响（面积越大，单价可能略低）
    area_factor = 1.0 if row['面积'] <= 90 else 0.95 if row['面积'] <= 120 else 0.9
    
    # 楼层影响（中间楼层更贵）
    floor_factor = 1.0
    if 5 <= row['楼层'] <= 15:
        floor_factor = 1.05
    elif row['楼层'] > 25:
        floor_factor = 0.95
    
    # 年份影响（越新越贵）
    year_factor = 1.0 + (2023 - row['建造年份']) * 0.01
    
    # 距离影响
    center_factor = max(0.5, 1.0 - row['距离市中心'] * 0.02)
    metro_factor = max(0.7, 1.0 - row['地铁距离'] * 0.1)
    
    # 学区房影响
    school_factor = 1.3 if row['学区房'] else 1.0
    
    # 装修程度影响
    decoration_factor = 0.8 + row['装修程度'] * 0.05
    
    # 地铁站数量影响
    metro_station_factor = 1.0 + row['地铁站数量'] * 0.02
    
    # 商业设施影响
    commercial_factor = 1.0 + row['商业设施'] * 0.01
    
    # 绿化率影响
    green_factor = 0.9 + row['绿化率'] * 0.002
    
    # 计算最终价格
    price_per_sqm = (base_price * district_multiplier[row['区域']] * 
                    area_factor * floor_factor * year_factor * 
                    center_factor * metro_factor * school_factor * 
                    decoration_factor * metro_station_factor * 
                    commercial_factor * green_factor)
    
    # 添加一些随机波动
    price_per_sqm *= np.random.uniform(0.9, 1.1)
    
    return price_per_sqm * row['面积']

# 计算房价
df['房价'] = df.apply(calculate_house_price, axis=1)

# 添加一些额外的真实特征
df['小区名称'] = [fake.company() + '小区' for _ in range(n_samples)]
df['房屋朝向'] = [fake.random_element(elements=['南', '北', '东', '西', '东南', '西南', '东北', '西北']) for _ in range(n_samples)]
df['建筑类型'] = [np.random.choice(['板楼', '塔楼', '板塔结合'], p=[0.4, 0.4, 0.2]) for _ in range(n_samples)]
df['产权年限'] = [fake.random_element(elements=[70, 50, 40]) for _ in range(n_samples)]

print(f"数据集创建完成！共 {len(df)} 条记录")
print(f"房价范围: {df['房价'].min():.0f} - {df['房价'].max():.0f} 元")
print(f"平均房价: {df['房价'].mean():.0f} 元")
print(f"房价中位数: {df['房价'].median():.0f} 元")
print(f"数据集形状: {df.shape}")
print("\n数据集前5行:")
print(df.head())

# 查看数据集基本信息
print("=== 数据集概览 ===")
print(f"数据集形状: {df.shape}")
print(f"特征列数: {df.shape[1]}")
print(f"样本数量: {df.shape[0]}")

print("\n=== 前5行数据 ===")
print(df.head())

print("\n=== 数据类型 ===")
print(df.dtypes)

print("\n=== 基本统计信息 ===")
print(df.describe())

print("\n=== 区域分布 ===")
print(df['区域'].value_counts())

print("\n=== 房价统计 ===")
print(f"最高房价: {df['房价'].max():,.0f} 元")
print(f"最低房价: {df['房价'].min():,.0f} 元")
print(f"平均房价: {df['房价'].mean():,.0f} 元")
print(f"房价中位数: {df['房价'].median():,.0f} 元")
print(f"房价标准差: {df['房价'].std():,.0f} 元")

# 数据探索性分析
print("数据集基本信息:")
print(df.describe())

# 可视化分析
plt.figure(figsize=(20, 16))

# 房价分布
plt.subplot(3, 3, 1)
plt.hist(df['房价'], bins=50, alpha=0.7, color='skyblue')
plt.title('房价分布')
plt.xlabel('房价 (元)')
plt.ylabel('频次')

# 面积与房价关系
plt.subplot(3, 3, 2)
plt.scatter(df['面积'], df['房价'], alpha=0.6, color='green')
plt.title('面积与房价关系')
plt.xlabel('面积 (平方米)')
plt.ylabel('房价 (元)')

# 房间数与房价关系
plt.subplot(3, 3, 3)
df.boxplot(column='房价', by='房间数', ax=plt.gca())
plt.title('房间数对房价的影响')
plt.xlabel('房间数')
plt.ylabel('房价 (元)')

# 学区房对房价的影响
plt.subplot(3, 3, 4)
df.boxplot(column='房价', by='学区房', ax=plt.gca())
plt.title('学区房对房价的影响')
plt.xlabel('学区房 (0:否, 1:是)')
plt.ylabel('房价 (元)')

# 距离市中心与房价关系
plt.subplot(3, 3, 5)
plt.scatter(df['距离市中心'], df['房价'], alpha=0.6, color='red')
plt.title('距离市中心与房价关系')
plt.xlabel('距离市中心 (公里)')
plt.ylabel('房价 (元)')

# 地铁距离与房价关系
plt.subplot(3, 3, 6)
plt.scatter(df['地铁距离'], df['房价'], alpha=0.6, color='orange')
plt.title('地铁距离与房价关系')
plt.xlabel('地铁距离 (公里)')
plt.ylabel('房价 (元)')

# 装修程度对房价的影响
plt.subplot(3, 3, 7)
df.boxplot(column='房价', by='装修程度', ax=plt.gca())
plt.title('装修程度对房价的影响')
plt.xlabel('装修程度')
plt.ylabel('房价 (元)')

# 建造年份与房价关系
plt.subplot(3, 3, 8)
plt.scatter(df['建造年份'], df['房价'], alpha=0.6, color='purple')
plt.title('建造年份与房价关系')
plt.xlabel('建造年份')
plt.ylabel('房价 (元)')

# 特征相关性热力图
plt.subplot(3, 3, 9)
# 选择数值列进行相关性分析
numeric_columns = df.select_dtypes(include=[np.number]).columns
correlation_matrix = df[numeric_columns].corr()
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', center=0, fmt='.2f')
plt.title('特征相关性热力图')

plt.tight_layout()
plt.show()

# 准备特征和目标变量
X = df.select_dtypes(include=[np.number]).drop('房价', axis=1)
y = df['房价']

# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 特征标准化
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# 训练线性回归模型
model = LinearRegression()
model.fit(X_train_scaled, y_train)

# 预测
y_pred = model.predict(X_test_scaled)

# 模型评估
mse = mean_squared_error(y_test, y_pred)
rmse = np.sqrt(mse)
r2 = r2_score(y_test, y_pred)

print("模型评估结果:")
print(f"均方误差 (MSE): {mse:,.2f}")
print(f"均方根误差 (RMSE): {rmse:,.2f}")
print(f"决定系数 (R²): {r2:.3f}")

# 特征重要性分析
feature_importance = pd.DataFrame({
    '特征': X.columns,
    '系数': model.coef_
}).sort_values('系数', key=abs, ascending=False)

print("\n特征重要性 (按系数绝对值排序):")
display(feature_importance)

# 可视化预测结果
plt.figure(figsize=(15, 10))

# 预测vs实际
plt.subplot(2, 3, 1)
plt.scatter(y_test, y_pred, alpha=0.6, color='blue')
plt.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 'r--', lw=2)
plt.xlabel('实际房价')
plt.ylabel('预测房价')
plt.title('预测vs实际房价')
plt.text(0.05, 0.95, f'R² = {r2:.3f}', transform=plt.gca().transAxes, 
         bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))

# 残差分析
plt.subplot(2, 3, 2)
residuals = y_test - y_pred
plt.scatter(y_pred, residuals, alpha=0.6, color='green')
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel('预测房价')
plt.ylabel('残差')
plt.title('残差分析')

# 残差分布
plt.subplot(2, 3, 3)
plt.hist(residuals, bins=30, alpha=0.7, color='orange')
plt.xlabel('残差')
plt.ylabel('频次')
plt.title('残差分布')

# 特征重要性可视化
plt.subplot(2, 3, 4)
feature_importance.plot(x='特征', y='系数', kind='barh', ax=plt.gca())
plt.title('特征重要性')
plt.xlabel('系数值')

# 预测误差分析
plt.subplot(2, 3, 5)
error_percentage = np.abs((y_test - y_pred) / y_test) * 100
plt.hist(error_percentage, bins=30, alpha=0.7, color='purple')
plt.xlabel('预测误差百分比 (%)')
plt.ylabel('频次')
plt.title('预测误差分布')

# 房价区间预测准确性
plt.subplot(2, 3, 6)
price_ranges = ['< 50万', '50-100万', '100-200万', '200-300万', '> 300万']
range_accuracy = []
for i, (min_price, max_price) in enumerate([(0, 500000), (500000, 1000000), 
                                          (1000000, 2000000), (2000000, 3000000), (3000000, float('inf'))]):
    mask = (y_test >= min_price) & (y_test < max_price)
    if mask.sum() > 0:
        accuracy = 1 - np.mean(np.abs((y_test[mask] - y_pred[mask]) / y_test[mask]))
        range_accuracy.append(max(0, accuracy))
    else:
        range_accuracy.append(0)

plt.bar(price_ranges, range_accuracy, color='skyblue')
plt.title('不同房价区间的预测准确性')
plt.ylabel('准确性')
plt.xticks(rotation=45)

plt.tight_layout()
plt.show()