数据拟合代码

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy import interpolate
from sklearn.metrics import r2_score, mean_squared_error
import warnings
warnings.filterwarnings('ignore')
# ==================== 1. 基础设置 ====================
# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体为黑体
plt.rcParams['axes.unicode_minus'] = False   # 解决负号'-'显示为方块的问题

# ==================== 2. 准备数据 ====================
data_points = np.array([
    [0.1, 0.03],
    [0.15, 0.033],
    [0.2, 0.038],
    [0.25, 0.038],
    [0.3, 0.045],
    [0.4, 0.052],
    [0.5, 0.055],
    [0.6, 0.065],
    [1.0, 0.08]
])
x_data = data_points[:, 0]
y_data = data_points[:, 1]

# ==================== 3. 定义拟合函数库 ====================
def linear(x, a, b):
    return a * x + b

def quadratic(x, a, b, c):
    return a * x**2 + b * x + c

def cubic(x, a, b, c, d):
    return a * x**3 + b * x**2 + c * x + d

def poly_4th(x, a, b, c, d, e):
    return a * x**4 + b * x**3 + c * x**2 + d * x + e

def poly_5th(x, a, b, c, d, e, f):
    return a * x**5 + b * x**4 + c * x**3 + d * x**2 + e * x + f

def power_law(x, a, b):
    return a * (x ** b)

def exponential(x, a, b, c):
    return a * np.exp(b * x) + c

def logarithmic(x, a, b):
    return a * np.log(x) + b

def rational(x, a, b, c):
    return (a * x) / (b + x) + c

def sigmoid(x, a, b, c, d):
    return a / (1. + np.exp(-b * (x - c))) + d

def fourier_3(x, a0, a1, b1, a2, b2, a3, b3, w):
    return (a0 + a1*np.cos(x*w) + b1*np.sin(x*w) +
                     a2*np.cos(2*x*w) + b2*np.sin(2*x*w) +
                     a3*np.cos(3*x*w) + b3*np.sin(3*x*w))

# 函数库字典
function_lib = {
    '线性': (linear, 2),
    '二次': (quadratic, 3),
    '三次': (cubic, 4),
    '四次多项式': (poly_4th, 5),
    '五次多项式': (poly_5th, 6),
    '幂律': (power_law, 2),
    '指数': (exponential, 3),
    '对数': (logarithmic, 2),
    '有理函数': (rational, 3),
    'S型生长': (sigmoid, 4),
    '傅里叶3阶': (fourier_3, 8),
}

# ==================== 4. 自动拟合与详细参数输出 ====================
print("="*80)
print("🤖 开始自动拟合测试 - 详细参数输出")
print("="*80)

# 存储所有结果
all_results = {}
best_name = ""
best_r2 = -np.inf

for name, (func, param_count) in function_lib.items():
    try:
        # 初始参数猜测
        p0 = [1] * param_count
        if name == '幂律':
            p0 = [0.1, 0.5]
        elif name == '指数':
            p0 = [0.03, 2.0, 0.01]
        elif name == 'S型生长':
            p0 = [0.08, 10, 0.5, 0.03]
        elif name == '傅里叶3阶':
            p0 = [0.05, 0.01, 0.01, 0.005, 0.005, 0.002, 0.002, 3.14]
        
        # 执行拟合
        params, _ = curve_fit(func, x_data, y_data, p0=p0, maxfev=5000)
        
        # 计算预测值和评估指标
        y_pred = func(x_data, *params)
        r2 = r2_score(y_data, y_pred)
        mse = mean_squared_error(y_data, y_pred)
        
        # 生成可读的函数表达式
        if name == '线性':
            expr = f"y = {params[0]:.6f} * x + {params[1]:.6f}"
        elif name == '二次':
            expr = f"y = {params[0]:.6f}*x² + {params[1]:.6f}*x + {params[2]:.6f}"
        elif name == '三次':
            expr = f"y = {params[0]:.6f}*x³ + {params[1]:.6f}*x² + {params[2]:.6f}*x + {params[3]:.6f}"
        elif name == '四次多项式':
            expr = f"y = {params[0]:.6f}*x⁴ + {params[1]:.6f}*x³ + {params[2]:.6f}*x² + {params[3]:.6f}*x + {params[4]:.6f}"
        elif name == '五次多项式':
            expr = f"y = {params[0]:.6f}*x⁵ + {params[1]:.6f}*x⁴ + {params[2]:.6f}*x³ + {params[3]:.6f}*x² + {params[4]:.6f}*x + {params[5]:.6f}"
        elif name == '幂律':
            expr = f"y = {params[0]:.6f} * x^{params[1]:.6f}"
        elif name == '指数':
            expr = f"y = {params[0]:.6f} * exp({params[1]:.6f}*x) + {params[2]:.6f}"
        elif name == '对数':
            expr = f"y = {params[0]:.6f} * ln(x) + {params[1]:.6f}"
        elif name == '有理函数':
            expr = f"y = ({params[0]:.6f}*x)/({params[1]:.6f} + x) + {params[2]:.6f}"
        elif name == 'S型生长':
            expr = f"y = {params[0]:.6f} / (1 + exp(-{params[1]:.6f}*(x-{params[2]:.6f}))) + {params[3]:.6f}"
        elif name == '傅里叶3阶':
            expr = f"y = {params[0]:.6f} + Σ[a_i*cos(i*w*x) + b_i*sin(i*w*x)], w={params[7]:.6f}"
        else:
            expr = f"y = f(x) with {param_count} parameters"
        
        # 存储结果
        all_results[name] = {
            'func': func,
            'params': params,
            'r2': r2,
            'mse': mse,
            'expr': expr,
            'y_pred': y_pred
        }
        
        # 详细输出
        print(f"\n📊 {name}:")
        print(f"  R² = {r2:.6f}, MSE = {mse:.8f}")
        print(f"  参数值: {np.round(params, 6)}")
        print(f"  函数表达式: {expr}")
        
        # 更新最佳结果
        if r2 > best_r2:
            best_r2 = r2
            best_name = name
            
    except Exception as e:
        print(f"\n❌ {name}: 拟合失败 - {str(e)[:50]}")
        all_results[name] = None

print("\n" + "="*80)
print(f"🏆 综合最优模型: {best_name} (R² = {best_r2:.6f})")
print("="*80)

# ==================== 5. LOWESS 非参数拟合 ====================
if LOWESS_AVAILABLE:
    print("\n" + "-"*70)
    print("尝试非参数拟合：局部加权回归 (LOWESS)...")
    try:
        # 使用LOWESS进行平滑，frac控制平滑带宽（值越大越平滑）
        frac_value = 0.7  # 可调整，范围通常在0.2到0.8之间
        lowess_smoothed = lowess(y_data, x_data, frac=frac_value, it=3, return_sorted=True)
        lowess_x, lowess_y = lowess_smoothed[:, 0], lowess_smoothed[:, 1]
        
        # === 新增：生成高密度平滑插值 ===
        # 使用三次样条插值生成更平滑的曲线
        lowess_smooth_interp = interpolate.interp1d(lowess_x, lowess_y, 
                                                    kind='cubic', 
                                                    bounds_error=False, 
                                                    fill_value='extrapolate')
        
        # 生成密集的插值点（从0.1到1.0，共500个点）
        x_smooth_dense = np.linspace(min(x_data), max(x_data), 600)
        y_smooth_dense = lowess_smooth_interp(x_smooth_dense)
        
        # === 新增：保存插值数据为TXT文件 ===
        lowess_smooth_data = np.column_stack((x_smooth_dense, y_smooth_dense))
        np.savetxt('lowess_smooth_interpolation.txt', lowess_smooth_data, 
                  fmt='%.6f', delimiter='\t', 
                  header='x\t y (LOWESS smoothed)', comments='')
        print(f"✅  LOWESS平滑插值数据已保存至: lowess_smooth_interpolation.txt")
        print(f"    数据点: {len(x_smooth_dense)}个, 范围: [{min(x_data):.2f}, {max(x_data):.2f}]")
        
        # 为计算指标，需要在原始x_data位置估计y值
        lowess_y_pred = lowess_smooth_interp(x_data)
        
        # 计算评估指标
        lowess_mse = mean_squared_error(y_data, lowess_y_pred)
        lowess_r2 = r2_score(y_data, lowess_y_pred)
        
        # 更新结果存储结构，包含平滑数据
        lowess_result = (lowess_x, lowess_y, lowess_y_pred, lowess_mse, lowess_r2, 
                        frac_value, x_smooth_dense, y_smooth_dense)
        
        print(f"✅  LOWESS (frac={frac_value}) | MSE = {lowess_mse:.6f} | R² (参考) = {lowess_r2:.6f}")
        print("   (注：LOWESS为平滑曲线，无显式函数参数)")
        
    except Exception as e:
        print(f"❌  LOWESS 拟合失败 | {e}")
        lowess_result = None
else:
    print("\n(跳过 LOWESS 拟合，因为 'statsmodels' 库不可用)")
    lowess_result = None




# ==================== 6. 新的可视化方案 ====================
print("\n🎨 正在生成分析图表...")

# 准备数据：整合参数拟合与LOWESS结果
all_viz_results = []  # 统一存储用于可视化的结果
# 1. 添加所有成功的参数拟合结果
for name, res in all_results.items():
    if res is not None:
        all_viz_results.append({
            'name': name,
            'type': '参数拟合',
            'r2': res['r2'],
            'mse': res['mse'],
            'expr': res['expr'],
            'data': res  # 指向原始结果字典
        })
# 2. 如果LOWESS拟合成功，也添加进来
if LOWESS_AVAILABLE and lowess_result is not None:
    # 解包LOWESS结果（现在包含8个元素）
    lowess_x, lowess_y, lowess_y_pred, lowess_mse, lowess_r2, frac_value, x_smooth, y_smooth = lowess_result
    all_viz_results.append({
        'name': f'LOWESS (frac={frac_value})',
        'type': '非参数拟合',
        'r2': lowess_r2,
        'mse': lowess_mse,
        'expr': '平滑曲线 (无显式公式)',
        'data': {
            'lowess_x': lowess_x, 
            'lowess_y': lowess_y,
            'lowess_y_pred': lowess_y_pred,
            'x_smooth': x_smooth,  # 新增：平滑插值x
            'y_smooth': y_smooth   # 新增：平滑插值y
        }
    })
# 按 R² 从高到低排序
all_viz_results.sort(key=lambda x: x['r2'], reverse=True)
func_names_viz = [item['name'] for item in all_viz_results]
r2_values_viz = [item['r2'] for item in all_viz_results]

# 6.1 第一部分：绘制包含LOWESS的所有函数的 R² 值对比图
print("\n1. 正在生成函数拟合优度对比图 (含LOWESS)...")
plt.figure(figsize=(12, 8))

colors = []
for item in all_viz_results:
    r2 = item['r2']
    # LOWESS用特殊颜色标识，参数拟合根据R²值着色
    if 'LOWESS' in item['name']:
        colors.append('purple')  # 紫色代表非参数拟合
    elif r2 > 0.95:
        colors.append('limegreen')
    elif r2 > 0.9:
        colors.append('orange')
    else:
        colors.append('lightcoral')

bars = plt.barh(func_names_viz, r2_values_viz, color=colors, edgecolor='black')
plt.xlabel('拟合优度 R²', fontsize=14)
plt.title('所有模型拟合效果对比 (R²值)', fontsize=16, fontweight='bold')
plt.xlim([0, 1.05])
plt.axvline(x=0.9, color='gray', linestyle='--', alpha=0.7, label='R²=0.9')
plt.axvline(x=0.95, color='gray', linestyle=':', alpha=0.7, label='R²=0.95')
plt.legend()
# 在条形末端添加R²数值标签
for bar, v in zip(bars, r2_values_viz):
    width = bar.get_width()
    plt.text(width + 0.01, bar.get_y() + bar.get_height()/2., 
             f'{v:.4f}', ha='left', va='center', fontsize=11)
plt.grid(axis='x', alpha=0.3)
plt.tight_layout()
plt.show()

# 6.2 第二部分：为每个函数单独绘制拟合图和残差图 (包括LOWESS)
print("\n2. 正在为每个模型生成详细分析图...")
x_plot_dense = np.linspace(min(x_data)*0.9, max(x_data)*1.0, 500)

for item in all_viz_results:
    model_name = item['name']
    model_type = item['type']
    r2 = item['r2']
    mse = item['mse']
    
    print(f"  正在处理: {model_name}")
    
    # 创建该模型的专属画布
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
    
    # 根据模型类型设置标题和颜色
    if 'LOWESS' in model_name:
        fig.suptitle(f'模型分析: {model_name} [非参数拟合]', fontsize=16, fontweight='bold')
        line_color = 'purple'
        point_color = 'darkviolet'
        # LOWESS模型的数据
        lowess_data = item['data']
        lowess_x = lowess_data['lowess_x']
        lowess_y = lowess_data['lowess_y']
        lowess_y_pred = lowess_data['lowess_y_pred']
    else:
        fig.suptitle(f'模型分析: {model_name} [参数拟合] (R² = {r2:.6f})', fontsize=16, fontweight='bold')
        line_color = 'blue'
        point_color = 'darkred'
        # 参数拟合模型的数据
        param_data = item['data']  # 这是原始结果字典，不需要再取['data']
        func = param_data['func']
        params = param_data['params']
        y_pred = param_data['y_pred']
    
    # ---- 子图1: 拟合曲线对比 ----
    ax1.scatter(x_data, y_data, color='darkred', s=80, zorder=5, 
                label='原始数据点', edgecolors='white', linewidth=1)
    
    # 根据模型类型绘制不同的拟合曲线
    if 'LOWESS' in model_name:
        # 获取平滑插值数据
        x_smooth = lowess_data['x_smooth']
        y_smooth = lowess_data['y_smooth']
        
        # 绘制平滑曲线（使用三次样条插值结果）
        ax1.plot(x_smooth, y_smooth, color=line_color, linewidth=2.5, 
                label=f'{model_name}平滑曲线', zorder=4)
        
        # 可选：同时绘制原始LOWESS点（半透明，仅作参考）
        # ax1.scatter(lowess_x, lowess_y, color=line_color, s=30, alpha=0.5, zorder=3)
        
        # 预测值已经在lowess_y_pred中
        residuals = y_data - lowess_y_pred
    else:
        # 参数拟合：直接计算函数值
        y_plot = func(x_plot_dense, *params)
        ax1.plot(x_plot_dense, y_plot, color=line_color, linewidth=2.5, 
                label=f'{model_name}拟合曲线', zorder=4)
        # 预测值已经在y_pred中
        residuals = y_data - y_pred
    
    ax1.set_xlabel('x', fontsize=12)
    ax1.set_ylabel('y', fontsize=12)
    ax1.set_title('原始数据与拟合/平滑曲线', fontsize=14)
    ax1.legend(loc='best')
    ax1.grid(True, alpha=0.3)
    
    # 在图中添加模型信息
    if 'LOWESS' in model_name:
        info_text = f"R2 (参考) = {r2:.6f}\nMSE = {mse:.2e}\n类型: 非参数平滑"
    else:
        expr_display = item['expr'] if len(item['expr']) < 80 else item['expr'][:77] + "..."
        info_text = f"R2 = {r2:.6f}\nMSE = {mse:.2e}\n表达式: {expr_display}"
    
    ax1.text(0.02, 0.98, info_text, 
             transform=ax1.transAxes, fontsize=10, 
             verticalalignment='top', 
             bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))
    
    # ---- 子图2: 残差分析图 ----
    ax2.scatter(x_data, residuals, color='green', s=80, zorder=5, alpha=0.8, label='残差点')
    ax2.axhline(y=0, color='r', linestyle='--', linewidth=2, alpha=0.7)
    ax2.fill_between(x_plot_dense, -0.005, 0.005, color='gray', alpha=0.1, label='±0.005区间')
    ax2.set_xlabel('x', fontsize=12)
    ax2.set_ylabel('残差 (y - ŷ)', fontsize=12)
    ax2.set_title('拟合残差分析', fontsize=14)
    ax2.legend(loc='upper right')
    ax2.grid(True, alpha=0.3)
    
    # 计算并显示残差统计信息
    mean_res = np.mean(residuals)
    std_res = np.std(residuals)
    max_abs_res = np.max(np.abs(residuals))
    ax2.text(0.02, 0.98, f"残差均值: {mean_res:.6f}\n残差标准差: {std_res:.6f}\n最大绝对残差: {max_abs_res:.6f}", 
             transform=ax2.transAxes, fontsize=11, 
             verticalalignment='top', 
             bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.8))
    
    plt.tight_layout()
    plt.show()