John | 曲

Reflection in Transition

第 18 章　理解实验数据

18 UNDERSTANDING EXPERIMENTAL DATA

曲政 / 2018-05-14

18.1 The Behavior of Springs

获得弹簧常数，实验科学家怎么做？

做法一：

知道胡克定律，用一个质量，测一个位移，然后计算，确信得到了弹簧常数。

要想保证上面的结论是正确的，下面的情况必须成立。

  • 确保实验操作完美，确信已经得到了 k。
  • 明确知道这次实验在弹性极限之内。

因此： 实验科学里很少这么做事情。

做法二：

在弹簧上悬挂多个重量不断增加的物体，并测量出弹簧每次拉伸的长度，然后绘制结果。

为什么要引入目标函数？

不论我们使用何种曲线 curve（包括直线）拟合 fit 数据，都需要某种方法确定哪条曲线才是数据的最佳拟合 best fit。这意味着我们需要定义一个目标函数 objective function，对曲线拟合数据的程度 how well the curve fit the data 提供一个定量的评价 quantitative assessment。如果我们定义了目标函数，那么找到最优拟合就可以明确表述为一个最优化问题 optimization problem：找到一条曲线，使目标函数值 the value of that function 最小 minimize（或最大 maximize）。

最小二乘作为目标函数有那两个优点？

最常用的目标函数称为最小二乘。令 observed 和 predicted 为两个同样长度的向量，observed 中是实际测量出来的点，predicted 中则是拟合曲线上相应的数据点。

那么，目标函数就可以定义为：

\[\sum^{len (\rm observed)-1}_{i=0}({\rm observed}[i]-{\rm predicted}[i])^2\]

（观察点值 - 预测点值）的平方和

  • 对 observed 和 predicted 中的点的差进行平方，使得二者之间较大的差比较小的差更重要。
  • 对差进行平方还可以消除差的正负影响。

为什么叫 “回归” 分析？

回归分析 regression analysis？

奠基者之一 K。皮尔逊（K. Pearson, 1856 1936) 在研究父母身高与其子女身高的遗传问题时，观察了 1 078 对夫妇，以每对夫妇的平均身高作为 x，而取他们的一个成年儿子的身高作为 y，将结果在平面直角坐标系上绘成散点图，发现趋势近乎一条直线。计算出的回归直线方程为：

\[ y = 33.73 + 0.516x \]

这种趋势及回归方程总的表明父母平均身高 I 每增加个单位，其成年儿子的身高 y 平均增加 0.516 个单位。这个结果表明，虽然高个子父辈的确有生高个子儿子的趋势，但父辈身高增加一个单位，儿子身高仅增加半个单位左右。反之，矮个子父辈的确有生矮个子儿子的趋势，但父辈身高减少一个单位，儿子身高仅减少半个单位左右。通俗地说，一群特高个子父辈（例如排球运动员）的儿子们在同龄人中平均仅为高个子，一群高个子父辈的儿子们在同龄人中平均仅为略高个子；一群特矮个子父辈的儿子们在同龄人中平均仅为矮个子，一群矮个子父辈的儿子们在同龄人中平均仅为略矮个子，即子代的平均高度向中心回归了。正是因为子代的身高有回到同龄人平均身高的这种趋势，才使人类的身高在一定时间内相对稳定，没有出现父辈个子高其子女更高，父辈个子矮其子女更矮的两极分化现象。这个例子生动地说明了生物学中 “种” 的概念的稳定性。正是为了描述这种有趣的现象，高尔顿引进了 “回归” 这个名词来描述父辈身高 x 与子辈身高 y 的关系。尽管 “回归” 这个名称的由来具有其特定的含义，而人们在研究大量的问题中，其变量 x 与 y 之间的关系并不总是具有这种 “回归” 的含义，但仍借用这个名词把研究变量 x 与 y 间统计关系的量化方法称为 “回归” 分析，也算是对高尔顿这位伟大的统计学家的纪念。

C.R.Rao 等在 Linear Models and Generalizations: Least Squares and Alternatives 中解释道 the literature meaning of REGRESSION is “to move in the backward direction”。

看以下两个陈述： S1: model generates data S2: data generates model

Rao 认为很明显陈述 S1 才是对的，因为模型实际上本来就是存在的只不过我们不知道 (model exists in nature but is unknown to the experimenter)。

先有模型所以我们知道 X 就能得到 Y：

先有模型 =》有了 X 就有 Y（S1）

而 “回归” 的意思就是我们通过收集 X 与 Y 来确定实际上存在的关系模型：

收集 X 与 Y =》确定模型（S2）

与 S1 相比，S2 就是一个 “回到” 模型的过程，所以就叫做 “regression”。

当然我认为 “回归” 也有一种回到平均水平的意思：

作者：Keven Howe 链接：https://www.zhihu.com/question/30123729/answer/47111877 来源：知乎 著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

ployfit 使用的算法称为线性回归 linear regression 都是线性的吗？

可以是高阶多项式回归 polynomial regression。

但是它预测出来的参数都是线性的。Although polynomial regression fits a nonlinear model to the data, the model is linear in the unknown parameters that it estimates.

不懂。

比如：

    model = pylab.polyfit (xVals, yVals, 3) = (a, b, c)

是线性的？

用一阶多项式线拟合弹簧数据？

def getData (fileName):
    """
    从有两列数据的文件中提取数据
    :param fileName: str
    :return (masses, distances): 质量和位移
    """
    dataFile = open (fileName, 'r')
    distances = []
    masses = []
    dataFile.readline () #discard header 先读一行，跨过标题行
    for line in dataFile:
        d, m = line.split ()
        distances.append (float (d))
        masses.append (float (m))
    dataFile.close ()
    return (masses, distances)
    
def labelPlot ():
    """为弹簧绘图做标题、x 坐标轴、y 坐标轴标注。"""
    pylab.title ('Measured Displacement of Spring')
    pylab.xlabel ('|Force| (Newtons)')
    pylab.ylabel ('Distance (meters)')

def plotData (fileName):
    xVals, yVals = getData (fileName)
    xVals = pylab.array (xVals)
    yVals = pylab.array (yVals)
    xVals = xVals*9.81  #acc. due to gravity
    pylab.plot (xVals, yVals, 'bo',
               label = 'Measured displacements')
    labelPlot ()
    
def fitData (fileName):
    xVals, yVals = getData (fileName)
    xVals = pylab.array (xVals)
    yVals = pylab.array (yVals)
    xVals = xVals*9.81 #get force x 轴数据变成 N 了
    pylab.plot (xVals, yVals, 'bo',
               label = 'Measured points')
    labelPlot ()                 
    a,b = pylab.polyfit (xVals, yVals, 1)
    estYVals = a*pylab.array (xVals) + b
    print ('a =', a, 'b =', b)
    pylab.plot (xVals, estYVals, 'r',
               label = 'Linear fit, k = '
               + str (round (1/a, 5)))
    pylab.legend (loc = 'best')
    
fitData ('springData.txt')

用三次多项式拟合弹簧数据，预测性能怎么样？

拟合得挺漂亮，但是远距离预测很不靠谱。

def fitData3 (fileName):
    xVals, yVals = getData (fileName)
    xVals = pylab.array (xVals)
    yVals = pylab.array (yVals)
    xVals = xVals*9.81 #get force
    pylab.plot (xVals, yVals, 'bo',
               label = 'Measured points')
    labelPlot ()
    model = pylab.polyfit (xVals, yVals, 3)
    add_point = pylab.array ([15])
    xVals = numpy.hstack ((xVals, add_point))
    estYVals = pylab.polyval (model, xVals)
    pylab.plot (xVals, estYVals, 'k--',)
    pylab.legend (loc = 'best')

fitData3 ('springData.txt')

竟然出现了负的位移，重物升到天蓬上面去了。

  • 在技术文献中，我们经常会看到类似的图，其中既有原始数据，也有一条拟合数据的曲线。然而，作者往往会将拟合曲线当作对真实情况的描述，而将原始数据看作实验误差。这是非常危险的一种做法。
  • 回忆一下，理论上 x 值和 y 值应该是线性关系，而不是三次关系。
  • 这种情况就是过拟合 overfitting 的典型例子。
  • 当模型过于复杂时，经常会出现过拟合，例如，相对于数据量，参数特别多的时候。
  • 发生过拟合时，拟合模型反映出的是数据中的噪声 capture noise，而不是数据中有意义的关系 meaningful relationship。
  • 过拟合模型的预测能力 predictive power 通常很弱，就像这个例子中一样。

可以丢弃实验数据吗？什么情况下可以？

绝对不能仅仅因为要得到更好的拟合而丢弃实验数据。

在这里，我们丢弃了最右侧的点是合理的，因为根据胡克的理论，弹簧是具有弹性极限的。但我们不能凭借这种理由丢弃数据中其他的点。

    xVals = pylab.array (xVals [:-6])
    yVals = pylab.array (yVals [:-6])

去掉 6 个点后，模型确实发生了变化：k 显著减小，而且线性拟合和三次拟合几乎没有区别。

18.2 The Behavior of Projectiles

从四组实验数据里提取一条弹丸轨迹。

import pylab


def getTrajectoryData (fileName):
    """
    从弹道数据文件中提取位置与高度，并做格式转换，不做计算。
    :param fileName: string，弹道数据文件名。
    :return: 距离，[高度 1 2 3 4] 数据列表
    """
    dataFile = open (fileName, 'r')
    distances = []
    heights1, heights2, heights3, heights4 = [],[],[],[]
    dataFile.readline () # 跨过标题行
    for line in dataFile:
        d, h1, h2, h3, h4 = line.split () # tuple 对应 list 赋值
        # str 转 float 后，推入列表
        distances.append (float (d))
        heights1.append (float (h1))
        heights2.append (float (h2))
        heights3.append (float (h3))
        heights4.append (float (h4))
    dataFile.close ()
    return (distances, [heights1, heights2, heights3, heights4])


def processTrajectories (fileName):
    """
    处理弹道数据，输出图形。
    :param fileName:
    :return: None，显示图形，用一次和二次多项式拟合弹道轨迹。
    """
    distances, heights = getTrajectoryData (fileName)
    numTrials = len (heights)
    distances = pylab.array (distances)
    #生成一个数组，用于存储每个距离的高度，并计算平均值
    totHeights = pylab.array ([0]*len (distances))
    for h in heights:
        totHeights = totHeights + pylab.array (h)
    meanHeights = totHeights/len (heights)
    # 绘图
    title = 'Trajectory of Projectile (Mean of ' + str (numTrials) + ' Trials)'
    pylab.figure (title)
    pylab.title (title)
    pylab.xlabel ('Inches from Launch Point')
    pylab.ylabel ('Inches Above Launch Point')
    pylab.plot (distances, meanHeights, 'ko')
    fit = pylab.polyfit (distances, meanHeights, 1)
    altitudes = pylab.polyval (fit, distances)
    pylab.plot (distances, altitudes, 'b', label = 'Linear Fit')
    fit = pylab.polyfit (distances, meanHeights, 2)
    # 用实验数据对应的点位 distances 来绘图，导致点少的区域上，曲线不够平滑。
    altitudes = pylab.polyval (fit, distances)
    pylab.plot (distances, altitudes, 'k:', label = 'Quadratic Fit')
    pylab.legend ()


processTrajectories ('launcherData.txt')
pylab.show ()
  • 二次曲线不够平滑没关系。
  • 二次曲线比一次要好，这很直观，但是好到什么程度？

为什么不能用最小二乘值来衡量不同拟合方式的优度？

拟合优度 goodness of fit 也就是预测的准确度 accuracy of these predictions。

回忆一下，拟合是通过最小化均方误差 minimizing the mean square error 而得到的，这说明我们可以通过均方误差来评价拟合优度。

这种方法的问题在于，尽管均方误差具有下界（0），但它没有上界 no upper bound。这意味着对于同一数据的两种拟合，我们可以使用均方误差来比较它们的相对优度 relative goodness of two fits，但很难用它衡量一个拟合的绝对优度 absolute goodness of a fit。

可以使用可决系数 coefficient of determination 计算一个拟合的绝对优度，可决系数通常写作 R2。令 yi 为第 i 个观测值，pi 为相应的模型预测值，μ 为观测值的均值，则：

\[ R^2=1-\frac {\sum_i (y_i-p_i)^2}{\sum_i (y_i-\mu)^2} \]

通过比较估计误差 estimation errors（分子 the numerator）和原始数据本身的变异性 variability of the original values（分母 the denominator），R2 可以表示在一个数据集中，有多大比例的变异性 the proportion of variability (relative to the mean) 是由于统计模型 by the statisitic model 通过拟合 provided by the fit 造成的 accounted for。

评价线性回归模型时，R2 的值总是位于 0 和 1 之间。如果 R2=1，模型就是对数据的完美拟合；如果 R2=0，那么模型预测值就与均值周围数据的分布方式没有任何联系。

def rSquared (measured, predicted):
    """ 计算线性回归模型的 R2 可决系数。
    假设 measured 为表示测量值的一维数组
        predicted 为表示预测值的一维数组
    返回可决系数 """
    # 用数组，语句简洁。
    # 计算分子，预测值与测量值差方和
    estimateError = ((predicted - measured)**2).sum ()
    # 计算原始数据方差。分母值是很大，貌似是给均方误差一个上界。
    meanOfMeasured = measured.sum ()/len (measured)
    variability = ((measured - meanOfMeasured)**2).sum ()
    return 1 - estimateError/variability
RSquare of linear fit = 0.0177433205441
RSquare of quadratic fit = 0.985765369287

简单地说，这个结果告诉我们，测量数据中只有不到 2% 的变异性 the variation in the measured data 可以用线性模型来解释 can be explained by，但超过 98% 的变异性可以由二次模型来解释。

  • 注意这个说法：测量数据的百分之多少的变异性可以用模型来解释。

用弹丸模型做什么？观摩什么通用模式？

比如求得弹丸落地时的水平速度。

def getHorizontalSpeed (quadFit, minX, maxX):
    """ 计算弹丸落地时的水平速度，打印输出。
    假设 quadFit 是二次多项式的系数
        minX 和 maxX 是用英寸表示的距离，抛物线与 x 轴的两个交点的坐标。
       返回以英尺 / 秒表示的水平速度 """
    inchesPerFoot = 12
    xMid = (maxX - minX)/2
    a,b,c = quadFit [0], quadFit [1], quadFit [2]
    yPeak = a*xMid**2 + b*xMid + c
    g = 32.16*inchesPerFoot #accel. of gravity in inches/sec/sec
    t = (2*yPeak/g)**0.5 #从最高点到目标高度所需时间，单位为秒
    print ('Horizontal speed =',
          int (xMid/(t*inchesPerFoot)), 'feet/sec')
          
# 计算落地时的水平速度，打印输出。
getHorizontalSpeed (fit, distances [-1], distances [0])

以上我们采取 work through 的一系列步骤 sequence of steps 是一种通用的模式 a common pattern：

  1. 首先进行实验，获得关于实体系统行为的数据； We started by performing an experiment to get some data about the behavior of a physical system.

  2. 然后通过计算找出描述系统行为的模型，并对模型质量进行评价； We then used computation to find and evaluate the quality of a model of the behavior of the system.

  3. 最后使用理论分析，设计一个简单的计算过程，推导出感兴趣的模型结果。 Finally, we used some theory and analysis to desigh a simple computation to derive an interesting consequence of the model.

18.3 Fitting Exponentially Distributed Data

指数分布的数据能用线性回归模拟吗？

使用 polyfit 找出数据模型的方法有个适用范围，即数据之间的联系可以使用形式为 y=base^(ax+b) 这样的公式来描述。

def fitExpData (xVals, yVals):
    """ 假设 xVals 和 yVals 是两个数值型数组，满足 yVals [i]=f (xVals (i))，这里的 f 是指数函数。
       返回 a、b、base，使得 log (f (x), base)==ax+b"""
    logVals = []
    for y in yVals:
        logVals.append (math.log (y, 2.0)) #求出以 2 为底的对数值
    fit = pylab.polyfit (xVals, logVals, 1)
    return fit, 2.0

把数据转换一个形式，就可以用已有的方法了。

18.4 When Theory is Missing

理论 theory 科学、实验 experimental 科学和计算 computational 科学之间的相互作用 interplay 怎样？

理论和计算都可以建立模型，实验和计算都可以生成数据，理论和实验互相启发和验证。

在理想世界中，我们可以进行完全可控的实验（例如在弹簧上悬挂重物），研究结果，进而回顾性地 retrospectively 构建 formulate 一个与这些结果一致的 consistent with 模型。然后，再进行新的实验（例如在同一个弹簧上悬挂另一个重物），并将新实验的结果与模型预测的结果进行比较。

不幸的是，在很多情况下，我们甚至不能进行可控的实验。举例来说，假设想构建一个研究利率如何影响股价的模型，我们几乎不可能去设定利率，然后再看看对股价有何影响 very few of us are in a position to set interest rates and see what happens。但从另一方面说，相关的历史数据可一点都不少。

在理论缺失的情况下，怎样让实验数据和计算模型良性互动？

两个类似的策略，共同点是，用数据生成模型，也用数据验证模型。

一种方法是将现有数据划分为训练集 training set 和保留集 holding set，保留集将来会作为测试集 test set 使用。

另外一种检验模型的方法使用从原始数据中随机选择的多个子集来训练模型，然后检查这些模型彼此之间的相似程度，称为交叉验证 cross validation。

以上，2018-05-14 11:23:58.