Python科学计算(第2版)
上QQ阅读APP看书,第一时间看更新

2.4.2 求和、平均值、方差

本节介绍的函数如表2-4所示。

表2-4 本节要介绍的函数

sum()计算数组元素之和,也可以对列表、元组等与数组类似的序列进行求和。当数组是多维时,它计算数组中所有元素的和。这里我们使用random.randint()模块中的函数创建一个随机整数数组。

    np.random.seed(42)
    a = np.random.randint(0,10,size=(4,5))
            a         np.sum(a)
    -----------------  ---------
    [[6, 3, 7, 4, 6],  96       
     [9, 2, 6, 7, 4],           
     [3, 7, 7, 2, 5],           
     [4, 1, 7, 5, 1]]           

如果指定axis参数,则求和运算沿着指定的轴进行。在上面的例子中,数组a的第0轴的长度为4,第1轴的长度为5。如果axis参数为1,则对每行上的5个数求和,所得的结果是长度为4的一维数组。如果参数axis为0,则对每列上的4个数求和,结果是长度为5的一维数组。即结果数组的形状是原始数组的形状除去其第axis个元素:

    np.sum(a, axis=1)   np.sum(a, axis=0)  
    -----------------  --------------------
    [26, 28, 24, 18]   [22, 13, 27, 18, 16]

当axis参数是一个轴的序列时,对指定的所有轴进行求和运算。例如下面的程序对一个形状为(2,3,4)的三维数组的第0和第2轴求和,得到的结果为一个形状为(3,)的数组。由于数组的所有元素都为1,因此求和的结果都是8:

    np.sum(np.ones((2, 3, 4)), axis=(0, 2))
    array([ 8.,  8.,  8.])

有时我们希望能够保持原数组的维数,这时可以设置keepdims参数为True:

    np.sum(a, 1, keepdims=True)  np.sum(a, 0, keepdims=True)
    ---------------------------  ---------------------------
    [[26],                       [[22, 13, 27, 18, 16]]     
     [28],                                                  
     [24],                                                  
     [18]]                                                  

sum()默认使用和数组的元素类型相同的累加变量进行计算,如果元素类型为整数,则使用系统的默认整数类型作为累加变量,例如在32位系统中使用32位整数作为累加变量。因此对整数数组进行累加时可能会出现溢出问题,即数组元素的总和超过了累加变量的取值范围。下面的程序计算数组a中每个元素占其所在行总和的百分比。在调用sum()函数时:

●设置dtype参数为float,这样得到的结果是浮点数组,能避免整数的整除运算。

●设置keepdims参数为True,这样sum()得到的结果的形状为(4,1),能够和原始数组进行广播运算。

    pa = a / np.sum(a, 1, dtype=float, keepdims=True) * 100
                        pa                     pa.sum(1, keepdims=True)
    ------------------------------------------  ------------------------
    [[ 23.08,  11.54,  26.92,  15.38,  23.08],[[ 100.],               
     [ 32.14,   7.14,  21.43,  25.  ,  14.29],[ 100.],               
     [ 12.5 ,  29.17,  29.17,   8.33,  20.83],[ 100.],               
     [ 22.22,   5.56,  38.89,  27.78,   5.56]][ 100.]]               

对很大的单精度浮点数类型的数组进行计算时,也可能出现精度不够的现象,这时也可以通过dtype参数指定累加变量的类型。在下面的例子中,我们对一个元素都为1.1的单精度数组进行求和,比较单精度累加变量和双精度累加变量的计算结果:

    np.set_printoptions(precision=8)
    b = np.full(1000000, 1.1, dtype=np.float32) # 创建一个很大的单精度浮点数数组
    b # 1.1无法使用浮点数精确表示,存在一些误差
    array([ 1.10000002,  1.10000002,  1.10000002, ...,  1.10000002,
            1.10000002,  1.10000002], dtype=float32)

使用单精度累加变量进行累加计算,误差将越来越大,而使用双精度浮点数则能够得到较精确的结果:

    np.sum(b)  np.sum(b, dtype=np.double)
    ---------  --------------------------
    1099999.3  1100000.0238418579        

上面的例子将产生一个新的数组来保存求和的结果,如果希望将结果直接保存到另一个数组中,可以和ufunc函数一样使用out参数指定输出数组,它的形状必须和结果数组的形状相同。

mean()求数组的平均值,它的参数与sum()相同。和sum()不同的是:对于整数数组它使用双精度浮点数进行计算,而对于其他类型的数组,则使用和数组元素类型相同的累加变量进行计算:

    np.mean(a, axis=1) # 整数数组使用双精度浮点数进行计算
    array([ 5.2,  5.6,  4.8,  3.6])
 
    np.mean(b)  np.mean(b, dtype=np.double)
    ----------  ---------------------------
    1.0999993   1.1000000238418579         

此外,average()也可以对数组进行平均计算。它没有out和dtype参数,但有一个指定每个元素权值的weights参数,可以用于计算加权平均数。例如有三个班级,number数组中保存每个班级的人数,score数组中保存每个班级的平均分,下面计算所有班级的加权平均分,得到整个年级的平均分:

    score = np.array([83, 72, 79])
    number = np.array([20, 15, 30])
    print np.average(score, weights=number)
    78.6153846154

相当于进行如下计算:

    print np.sum(score * number) / np.sum(number, dtype=float)
    78.6153846154

std()和var()分别计算数组的标准差和方差,有axis、out、dtype以及keepdims等参数。方差有两种定义:偏样本方差(biased sample variance)和无偏样本方差(unbiased sample variance)。偏样本方差的计算公式为:

而无偏样本方差的公式为:

当ddof参数为0时,计算偏样本方差;当ddof为1时,计算无偏样本方差,默认值为0。下面我们用程序演示这两种方差的差别。

首先产生一个标准差为2.0、方差为4.0的正态分布的随机数组。我们可以认为总体样本的方差为4.0。假设从总体样本中随机抽取10个样本,我们分别计算这10个样本的两种方差,这里我们用一个二维数组重复上述操作100000次,然后计算所有这些方差的期望值:

    a = nr.normal(0, 2.0, (100000, 10)) 
    v1 = np.var(a, axis=1, ddof=0) #可以省略ddof=0
    v2 = np.var(a, axis=1, ddof=1)
       np.mean(v1)         np.mean(v2)    
    ------------------  ------------------
    3.6008566906846693  4.0009518785385216

可以看到无偏样本方差的期望值接近于总体方差4.0,而偏样本方差比4.0小一些。

偏样本方差是正态分布随机变量的最大似然估计。如果有一个样本包含n个随机数,并且知道它们符合正态分布,通过该样本可以估算出正态分布的概率密度函数的参数。所估算的那组正态分布参数最符合给定的样本,就称为最大似然估计。

正态分布的概率密度函数的定义如下,其中μ表示期望,σ2表示方差:

所谓最大似然估计,就是找到一组参数使得下面的乘积最大,其中x1为样本中的值:

f(x1 )f(x2 )…f(xn )

专业术语总是很难理解,下面我们还是用程序来验证:

    def normal_pdf(mean, var, x):
        return 1 / np.sqrt(2 * np.pi * var) * np.exp(-(x - mean) ** 2 / (2 * var))
    
    nr.seed(42)
    data = nr.normal(0, 2.0, size=10)                         ❶
    mean, var = np.mean(data), np.var(data)                   ❷
    var_range = np.linspace(max(var - 4, 0.1), var + 4, 100)  ❸
    
    p = normal_pdf(mean, var_range[:, None], data)            ❹
    p = np.product(p, axis=1)                                 ❺

normal_pdf()为计算正态分布的概率密度的函数。❶产生10个正态分布的随机数。❷计算其最大似然估计的参数。❸以最大似然估计的方差为中心,产生一组方差值。❹用正态分布的概率密度函数计算每个样本、每个方差所对应的概率密度。由于使用了广播运算,得到的结果p是一个二维数组,它的第0轴对应var_range中的各个方差,第1轴对应data中的每个元素。❺沿着p的第1轴求所有概率密度的乘积。product()和sum()的用法类似,用于计算数组所有元素的乘积。

下面绘制var_range中各个方差对应的似然估计值,并用一条竖线表示偏样本方差。由图2-8可以看到偏样本方差位于似然估计曲线的最大值处。

图2-8 偏样本方差位于似然估计曲线的最大值处

    import pylab as pl
    pl.plot(var_range, p)
    pl.axvline(var, 0, 1, c="r")
    pl.show()