关于#python#的问题：这里是用什么方法选取py_init的

Question

想问一下下面程序里面的
py_init = np.array([Y[np.abs(X - x) < (xmax - xmin) * 0.01].mean() for x in px_init])
这句话具体含义，为什么要np.abs(X - x) < (xmax - xmin) * 0.01？
并且为什么要乘以0.01？
这里是用什么方法选取py_init的？

%matplotlib inline
import numpy as np
import pylab as pl
 
def make_test_data(seg_count, point_count):
    x = np.random.uniform(2, 10, seg_count)
    x = np.cumsum(x)
    x *= 10 / x.max()
    y = np.cumsum(np.random.uniform(-1, 1, seg_count))
    X = np.random.uniform(0, 10, point_count)
    Y = np.interp(X, x, y) + np.random.normal(0, 0.05, point_count)
    return X, Y
   
from scipy import optimize
 
def segments_fit(X, Y, count):
    xmin = X.min()
    xmax = X.max()
 
    seg = np.full(count - 1, (xmax - xmin) / count)
 
    px_init = np.r_[np.r_[xmin, seg].cumsum(), xmax]
    py_init = np.array([Y[np.abs(X - x) < (xmax - xmin) * 0.01].mean() for x in px_init])
 
    def func(p):
        seg = p[:count - 1]
        py = p[count - 1:]
        px = np.r_[np.r_[xmin, seg].cumsum(), xmax]
        return px, py
 
    def err(p):
        px, py = func(p)
        Y2 = np.interp(X, px, py)
        return np.mean((Y - Y2)**2)
 
    r = optimize.minimize(err, x0=np.r_[seg, py_init], method='Nelder-Mead')
    return func(r.x)
 
X, Y = make_test_data(10, 2000)
px, py = segments_fit(X, Y, 8)
 
pl.plot(X, Y, ".")
pl.plot(px, py, "-or");
 
 

Answer 1

引自免费微信小程序：皆我百晓生

这段Python代码是用来实现分段线性插值的。它首先创建了一些测试数据（X, Y），然后通过优化算法找到最佳的分段点（px）和每个段上的线性函数的y截距（py）。最终，程序会绘制原始数据点（.）和拟合的分段线（-or）。

现在，我们逐句分析这个代码，解释其含义和目的：

1. py_init = np.array([...]): 这行代码是初始化py_init数组，用于存储每一段线性部分的初始y截距估计。对于每一个px_init中的x值，它计算与之最接近的X范围内的Y值的平均值。

   具体来说：

   np.abs(X - x) < (xmax - xmin) * 0.01
   

   是一个布尔条件数组，表示X值与当前x值之间的绝对差小于(xmax - xmin)（整个x范围）的0.01倍。这意味着它会选择那些离x值非常近的X值对应的Y值进行平均。

   乘以0.01是为了确保选择的Y值是在给定x值附近的数据点，这样可以得到更准确的初始y截距估计。0.01是一个经验值，可以根据实际数据调整以获得更好的拟合效果。

2. func(p): 它是一个内部函数，接收参数p，并根据这些参数计算出分段的x点（px）和y截距（py）。

3. err(p): 这个函数定义了优化的目标函数，即误差平方和。输入参数p包含了所有分段长度（seg）和y截距（py），函数返回的是拟合后的线性函数与原数据点的残差平方和。

4. r = optimize.minimize(...): 使用Scipy库的optimize.minimize函数对目标函数err进行最小化。这里的method='Nelder-Mead'指定了优化算法为Nelder-Mead法，这是一种全局优化方法，适合没有导数的情况。

5. return func(r.x): 最后，返回优化得到的最佳参数对应的px和py，用于绘制拟合曲线。

以下是完整的代码，已经注释了解释：

import numpy as np
import matplotlib.pyplot as plt

def make_test_data(seg_count, point_count):
    # 创建测试数据
    ...

def segments_fit(X, Y, count):
    # 初始化变量
    xmin = X.min()
    xmax = X.max()

    seg = (xmax - xmin) / count  # 计算每个段的宽度

    # 初始化px_init，包含起始点、段宽以及终止点
    px_init = np.r_[np.r_[xmin, seg].cumsum(), xmax]

    # 初始化py_init，对于每个px_init中的x值，取X范围内最接近的Y值的平均
    py_init = np.array([Y[np.abs(X - x) < (xmax - xmin) * 0.01].mean() for x in px_init])

    def func(p):
        # 根据参数计算px和py
        seg = p[:count - 1]
        py = p[count - 1:]
        px = np.r_[np.r_[xmin, seg].cumsum(), xmax]
        return px, py

    def err(p):
        # 目标函数：误差平方和
        px, py = func(p)
        Y2 = np.interp(X, px, py)
        return np.mean((Y - Y2)**2)

    # 优化
    r = optimize.minimize(err, x0=np.r_[seg, py_init], method='Nelder-Mead')

    # 返回最优解对应的px和py
    return func(r.x)

# 测试数据
X, Y = make_test_data(10, 2000)

# 拟合
px, py = segments_fit(X, Y, 8)

# 绘制
plt.plot(X, Y, ".", label="Data Points")
plt.plot(px, py, "-or", label="Fitted Segments")
plt.legend()
plt.show()

这就是这段代码的工作原理和每一部分的具体作用。如果你还有其他疑问或需要进一步的解释，请随时提问。