概述
MatLab在科学计算中，具有很强大的功能，常用的计算方法是计算一组数据的平均数，标准差，方差。
1.计算一个数组的平均值
使用 mean() 函数
示例1：输入：a = [ 1 2 3;4 5 6] 输入：mean(a) 输出：[ 2.5,3.5,4.5]解析：默认是先求得每一列的的和，然后对每一个列求平均值
示例2：
输入：mean(mean(a))
解析：此时所输出的就是整个数组的平均值，需要加以注意2.计算一个数组的方差
示例：
输入：var(a) % a 采用上面的示例 输出：[4.5,4.5,4.5] 解析：依然是求得每一列的方差 解决方法：将矩阵转换为列向量 输入：a = [1;2;3;4;5;6] 输入：var(a) 输出：[3.500]3.计算一个数组的标准差
和计算方差一样，将数组转换为列向量，使用std(a)，进行求解
输入：a = [1;2;3;4;5;6] 输入：std(a) 输出：[1.8708]
将Mat类型转为IplImage类型并计算它的均值和标准差，其中frame为Mat类型。
IplImage ipl;
ipl = IplImage(frame);
CvScalar mean, std_dev;
cvAvgSdv(&ipl, &mean, &std_dev);
cout << "该图像的平均值为:" << mean.val[0] << "\t标准差为:" << std_dev.val[0] << endl;
本文主要参考 ISO 13528 的附件 C 的 3 条款。算法原理
应用此算法计算得到数据平均值和标准差的稳健值。稳健性是估计算法的特点，而不是其产生的估计值的特点，因此严格来说，称由此算法计算的平均值和标准差是稳健的是不确切的。然而，为避免使用繁琐的术语，“稳健均值”和“稳健标准差”应理解为利用稳健算法计算的总体均值和总体标准差的均值估计。
设数据为 x i , i ∈ ( 1 , 2 , ⋯ , n ) x_i, i\in(1,2,\cdots,n) xi,i∈(1,2,⋯,n)，记稳健平均值为 x ∗ , s ∗ x^*, s^* x∗,s∗，首先是计算初始值：
x ∗ = med x i s ∗ = 1.483 × med ∣ x i − x ∗ ∣ \begin{array}{c} x^{*}=\operatorname{med} x_{i} \\ s^{*}=1.483 \times \operatorname{med}\left|x_{i}-x^{*}\right| \end{array} x∗=medxis∗=1.483×med∣xi−x∗∣
对每个 x i x_i xi，有：
x i ∗ = { x ∗ − δ , 若 x i < x ∗ − δ x ∗ + δ , 若 x i > x ∗ + δ x i , 其他 x_{i}^{*}=\left\{\begin{array}{cc} x^{*}-\delta, & \text { 若 } x_{i}<x^{*}-\delta \\ x^{*}+\delta, & \text { 若 } x_{i}>x^{*}+\delta \\ x_{i}, & \text { 其他 } \end{array}\right. xi∗=⎩⎨⎧x∗−δ,x∗+δ,xi, 若 xi<x∗−δ 若 xi>x∗+δ 其他
其中： δ = 1.5 s ∗ \delta=1.5 s^{*} δ=1.5s∗，再次计算：
x ∗ = ∑ x i ∗ / p s ∗ = 1.134 ∑ ( x i ∗ − x ∗ ) 2 / ( p − 1 ) \begin{array}{c} x^{*}=\sum x_{\mathrm{i}}^{*} / p \\ s^{*}=1.134 \sqrt{\sum\left(x_{\mathrm{i}}^{*}-x^{*}\right)^{2} /(p-1)} \end{array} x∗=∑xi∗/ps∗=1.134∑(xi∗−x∗)2/(p−1)重复：
x i ∗ = { x ∗ − δ , 若 x i < x ∗ − δ x ∗ + δ , 若 x i > x ∗ + δ x i , 其他 x_{i}^{*}=\left\{\begin{array}{cc} x^{*}-\delta, & \text { 若 } x_{i}<x^{*}-\delta \\ x^{*}+\delta, & \text { 若 } x_{i}>x^{*}+\delta \\ x_{i}, & \text { 其他 } \end{array}\right. xi∗=⎩⎨⎧x∗−δ,x∗+δ,xi, 若 xi<x∗−δ 若 xi>x∗+δ 其他 直到 s ∗ s^* s∗ 的 第三位有效数字和 x ∗ x^* x∗ 的对应数字在连续两次迭代中不变。
代码
# -*- coding: utf-8 -*- import pandas as pd import numpy as np def location_corresponding(x_str, s_str): ''' s 的第三位有效数字，以及对应 x 的相应的数字 ''' digit_location = s_str.find('.') if digit_location >= 3: # 如果小数点在字符串的第 4 位或大于第四位，也即数字至少在 100 以上 # 此时第三位有效数字，肯定在小数点前。 significant_num = s_str[2] # 有效位数于小数点的相对位置 # 0 代表在小数点前， 3-1 是指从小数点前数起的位数 # sig_loc 的第一位表示有效数字在小数点前，还是后 # 第二位代表有效数字在小数点的“距离” str_len = len(s_str[:digit_location]) sig_loc = (0, str_len-3-1) elif s_str[0] == '0': # 若整数部分是0，则有效数字的位置在小数点后 # 刨除整数和小数点部分 s_without_int = s_str[digit_location:] # 若小数点后有 0，则不将 0 计入有效数字位 if '.0' in s_without_int: digit_location += 1 while '.00' in s_without_int: # 若小数点后有多个0，也不计入 s_without_int = s_without_int.replace('.00', '.0') digit_location += 1 try: # 若位数不够，如 0.0，则有效数字为 0、 significant_num = s_str[digit_location+3] except: significant_num = '0' # 1 表示有效数字在小数点后 sig_loc = (1, 3) elif digit_location == -1 and len(s_str) >= 3: # 若只有整数，没有小数部分，且整数部分大于 100，即有超过三位数 # 则直接取第三位 str_len = len(s_str) sig_loc = (0, str_len-3) significant_num = s_str[2] elif digit_location == -1 and len(s_str) < 3: # 若小于 100，则有效数字为0 significant_num = '0' else: # 若整数部分小于 2 位，且整数部分大于 0 sig_loc = (1, 3-digit_location) significant_num = s_str[sig_loc[1]+digit_location] # x 对应的数字 x_digit_location = x_str.find('.') if sig_loc[0] == 0: if x_digit_location == -1 and len(x_str) >= 3: x_significant_num = x_str[2] elif x_digit_location == -1 and len(x_str) < 3: x_significant_num = '0' else: x_significant_num = x_str[x_digit_location-sig_loc[1]] print(sig_loc[1]) print('前') else: try: x_significant_num = x_str[x_digit_location+sig_loc[1]] except: x_significant_num = '0' return x_significant_num, significant_num def coverage_critiria(x_list, s_list): ''' 收敛准则 其中 s_list 是一个长度为 3 的 list, 包含当前迭代的 s* 和之前两个迭代的 s* 其中 x_list 也一样 难点在于：如何找出 s 的第三位有效数字，对应 x 的数位呢？ 这里的解决办法是：找出 s 三位有效数字，在小数点的位置，从而应用于 x 中 ''' s_numbers = [] x_numbers = [] for i in range(3): # 连续两位不变，故需要进行 3 此迭代。 s = s_list[i] # s 的字符串 s_str = str(s) x = x_list[i] # x 的字符串 x_str = str(x) # 找出 s* 的第三位有效数字，和 x* 对应的数字。 x_sig_num, s_sig_num = location_corresponding(x_str, s_str) s_numbers.append(s_sig_num) x_numbers.append(x_sig_num) # 稳健标准差的第三位有效数字，连续两次不变，且稳健平均值的对应数字亦连续两次不变，则判断为收敛 # 函数返回 True。否则返回 False s_equal_flag = (s_numbers[0] == s_numbers[1]) and (s_numbers[0] == s_numbers[2]) x_equal_flag = (x_numbers[0] == x_numbers[1]) \ and (x_numbers[0] == x_numbers[2]) if s_equal_flag and x_equal_flag: return True else: return False def perform_algorithm_A(x): ''' x 是一个向量 ''' if type(x) is list: x = np.array(x) # 稳健平均和稳健标准差初始值 x_star = np.median(x) x_diff = np.absolute(x-x_star) s_star = 1.483*np.median(x_diff) # 上界和下界的初始值 delta = 1.5*s_star higher_bound = x_star + delta lower_bound = x_star - delta x_list = [] s_list = [] while True: x_tmp = np.copy(x) x_tmp[x_tmp>higher_bound] = higher_bound x_tmp[x_tmp<lower_bound] = lower_bound x_star = np.mean(x_tmp) s_star = 1.134*np.std(x_tmp) x_list.append(x_star) s_list.append(s_star) if len(x_list) == 4: x_list.pop(0) s_list.pop(0) # 达到收敛条件 if coverage_critiria(x_list, s_list): return x_tmp, x_list[2], s_list[2], higher_bound, lower_bound # 计算上下界 delta = 1.5*s_star higher_bound = x_star + delta lower_bound = x_star - delta if __name__ == '__main__': x = [927,952,977,995,915,962,966,950,969,949,961,940,1002,956,960,943] x, x_star, s_star, higher_bound, lower_bound = perform_algorithm_A(x) print('算法 A 收敛后，数据变为： ', x) print('稳健平均值： ', x_star) print('稳健标准差： ', s_star)运行结果如下：
用numpy函数可实现快速计算
n [23]: import numpy as np In [24]: data = [1, 2, 3] In [25]: np.mean(data) Out[25]: 2.0 In [26]: np.sum(data) Out[26]: 6 In [27]: np.max(data) Out[27]: 3 In [28]: np.min(data) Out[28]: 1 In [29]: np.std(data) Out[29]: 0.81649658092772603In [32]: data = np.array([1, 2, 3]) In [33]: data.min() Out[33]: 1 In [34]: data.max() Out[34]: 3 In [35]: data.std() Out[35]: 0.81649658092772603 In [36]: data.mean() Out[36]: 2.0
转载于:https://www.cnblogs.com/Zhu-Xueming/p/8650497.html
