python均值检验函数
ttest检验原理公式
pytho的显著线检验一般都使用SciPy程序包，当检验两组样本均值时，有如下几种方法：
from scipy import stats
import numpy as np
tvalue,pvalue = stats.ttest_ind_from_stats(mean1, std1, numb1, \
mean2, std2, numb2, equal_var=True, alternative='two-sided')
# 输入两组数据的均值、方差和样本量以进行显著性检验
# equal_var=True表示两组数据具有相近的总体方差（默认）
# 若False，表示两组数据不具有相近的总体方差
# 一样的输入，若equal_var=True，则得到的pvalue偏小，更容易通过显著性检验
tvalue,pvalue = stats.ttest_ind(data1,data2, axis=0, equal_var=True,\
nan_policy='propagate', permutations=None, random_state=None, \
alternative='two-sided', trim=0)
# 也是用于检验两组样本均值是否显著差异，但该函数需要将两个样本全部输入，
# 若样本量较大，则比较占用内存，因此个人还是更喜欢上面那个函数
# 关于各参数的具体介绍可以看
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ttest_ind.html
# 此外 SciPy 也有提供检验两组样本总体方差是否相等的函数levene与bartlett
# 输出统计值和pvalue，当pvalue大时，说明两组样本的总体方差差异大
# 相同的两组数据，bartlett的pvalue大于levene的
print(stats.bartlett(data1, data2))
print(stats.levene(data1, data2))
siglvl = 0.05 # 95% 显著性水平
var.values=np.ma.array(var.values,mask=(pvalue>siglvl))
mask = np.array([pu>siglvl,pv>siglvl]).all(axis=0)
uwnd.values=np.ma.array(uwnd.values,mask=mask)
vwnd.values=np.ma.array(vwnd.values,mask=mask)
参考资料：
- https://github.com/scipy/scipy/blob/135b734994a512334055a97a0f2fae8f0990887b/scipy/stats/_stats_py.py
- https://blog.csdn.net/yuqiuwang929/article/details/103821349?spm=1001.2101.3001.6650.1&utm_medium=distribute.pc_relevant.none-task-blog-2%7Edefault%7ECTRLIST%7Edefault-1.no_search_link&depth_1-utm_source=distribute.pc_relevant.none-task-blog-2%7Edefault%7ECTRLIST%7Edefault-1.no_search_link
动态计算样本均值和方差
由于现在需要对两组逐6小时数据做均值检验，样本量高达10000多。没有办法一次性将所有数据都读出，只能逐年读出。而两组样本的均值检验需要两组样本的均值、方差、样本量，因此需要动态计算样本均值和方差，公式推导如下：
假设前n年的样本量为$M_n$，均值 $\bar E_n$, 方差 $V_n$, 第n的样本量为$m_n$,样本元素用为$a_i$表示
样本量的变化：$M_n=M_{n-1}+m_n$
均值$\bar E_n$递推公式:
\begin{split}
\bar E_n &= \frac{(M_n-m_n)\bar E_{n-1} + \sum_{i=1}^{m_n} a_i}{M_n} \\
&=\bar E_{n-1} + \frac{\sum_{i=1}^{m_n} a_i - m_n\bar E_{n-1}}{M_n}
\end{split}
方差$V_n$递推公式:
\begin{split}
V_n &= \frac{M_{n-1}{V_{n-1}+(\bar E_{n}-\bar E_{n-1})^2}+ \sum_{i=1}^{m_{n}} (a_{i}-\bar E_n)^2}{M_n}
\end{split}
\begin{split}
M_nV_n - M_{n-1}V_{n-1} &= \sum_{i=1}^{M_{n}} (A_{i}-\bar E_n)^2-\sum_{i=1}^{M_{n-1}} (A_{i}-\bar E_{n-1})^2\\
&=\sum_{i=1}^{M_{n-1}} (A_{i}-\bar E_n)^2+\sum_{i=1}^{m_{n}} (a_{i}-\bar E_n)^2-\sum_{i=1}^{M_{n-1}} (A_{i}-\bar E_{n-1})^2\\
&=\sum_{i=1}^{m_{n}} (a_{i}-\bar E_n)^2+M_{n-1}(\bar E_n^2-\bar E_{n-1}^2)-2(\bar E_{n}-\bar E_{n-1})M_{n-1}\bar E_{n-1}\\
&=\sum_{i=1}^{m_{n}} (a_{i}-\bar E_n)^2+M_{n-1}(\bar E_{n}-\bar E_{n-1})^2
\end{split}
python脚本：
def mean_vari_dynamic(mean0,vari0,numb0,samp):
print("numb0 %d, mean0[1,12,12] %.2f, vari0[1,12,12] %.2f"
%(numb0,mean0[1,12,12],vari0[1,12,12]))
numb1 = numb0 + len(samp)
mean1 = mean0 + (np.sum(samp,axis=0)-len(samp)*mean0)/numb1
vari1 = (numb0*(vari0+np.power(mean1-mean0,2))+
np.sum(np.power(np.subtract(samp,mean1),2),axis=0))/numb1
return mean1,vari1,numb1
mean = np.zeros( [len(levc),len(ilat),len(ilon)],dtype=float )
vari = np.zeros( [len(levc),len(ilat),len(ilon)],dtype=float )
allt = 0
for year in range(1979,2021,1):
ds = xr.open_dataset("%s%s/ERA5_NH_%s_%d.nc"%(datapath,varname,varname,year))
term = ds[varname].sel(time=ds.time.dt.month.isin(months[nm]),
level=levc,longitude=ilon,latitude=ilat)
mean, vari, allt = mean_vari_dynamic(mean,vari,allt,term)
支付宝鼓励 
鸡腿鸡腿