常用气象统计方法（一）

不想去气象局

写在前面：今天看到气象家园上有小伙伴需要一些基础的统计分析的内容，所以借此整理一下这学期气象条件方法的内容，我把这篇文章写到的方法都列在前面方便大家查询这篇文章是否包含你想要的内容，后续如果有添加的话就开新的文章。这篇文章包含的方法有：简单相关系数、自相关系数、落后交叉相关系数、偏相关系数、一元线性回归、多元线性回归、滑动平均、累积距平、EOF、合成分析。

下面按照实习课程的顺序来进行整理，篇幅原因仅保留关键代码，如果需要完整代码的可以在后台告诉我需要哪个脚本的例子。

简单相关系数

import numpy as np

from scipy.stats import pearsonr

import pandas as pd

filename='C:\\Users\\DELL\\Desktop\\NUIST\\学科\\气象统计方法\\实习\\program2\\data2.xlsx'

df=np.array(pd.read_excel(filename))

r,p=pearsonr(df[:,1], df[:,2])

这个数据是存放在excel表里的，因为大家用的nc文件比较多，所以这里的excel文件的读取也可以留个备忘，r是相关系数，p是置信度，如果是二维的图，95%打点的时候用p就很方便

自相关系数

import statsmodels.tsa.api as smt

filename='C:\\Users\\DELL\\Desktop\\NUIST\\学科\\气象统计方法\\实习\\program2\\data2.xlsx'

df=pd.read_excel(filename)

yr_cor=smt.stattools.acf(df.values[:,1],nlags=10)

同样也是excel里的数据，nlags表示滞后时间，都是有现成的库的，顺便聊几句，直接调用跟我气象统计方法的老师的要求不太一样，他要我们自己写代码，所以我最后实习报告的分也就不高了。。

落后交叉相关系数

def cross_r(x1,x2,n):

   n1=np.array(x1.shape)

   n1=int(n1)

   r,p=pearsonr(x1[:n1-n], x2[n:])

return r,p

filename='C:\\Users\\DELL\\Desktop\\NUIST\\学科\\气象统计方法\\实习\\program3\\data3.xlsx'

df=pd.read_excel(filename)

da=np.array(df)

for i in range(10):

   r1,p1=cross_r(da[:,1],da[:,2],i+1)

r2,p2=cross_r(da[:,1],da[:,3],i+1)

其实我比较喜欢找现成的库，不喜欢自己写代码（不论多简单），除非实在找不到资源才会自己写。但是这个最后还是自己写了。Cross函数里的n表示滞后时间，函数使用的前提是x1,x2都是一维的，而且数据一一对应。

偏相关系数

r_ab,p=pearsonr(da[:,1],da[:,2])

r_ac,p=pearsonr(da[:,1],da[:,3])

r_bc,p=pearsonr(da[:,2],da[:,3])

r_ab_c=(r_ab-r_ac*r_bc)/(((1-r_ac**2)**0.5)*((1-r_bc**2)**0.5))

这个偏相关系数没啥好说的，就是用公式。

一元线性回归

import numpy as np

from sklearn.linear_model import LinearRegression

import pandas as pd

from statsmodels.formula.api import ols

filename='C:\\Users\\DELL\\Desktop\\NUIST\\学科\\气象统计方法\\实习\\program4\\data4.xlsx'

df=pd.read_excel(filename)

year=df['year']

tem=df['temperature']

index=df['index']

reg = LinearRegression()

tem=np.array(tem)

index=np.array(index)

reg.fit(index.reshape(-1,1),tem.reshape(-1,1))

data = pd.DataFrame({'x':index, 'y':tem})

model = ols('y~x', data).fit()

print(model.summary())

reshape(-1,1)是把数据转置成(n,1)这样的格式。利用ols可以直接输出一元线性回归的各种信息，如下图。但是一般来说自己用的时候不会用到ols，截距和系数可以直接用reg.coef_和reg.intercept_来输出，显著性检验的话一般用相关系数的p，因为一元线性回归的显著性和相关系数的显著性检验本质上是一样的。

多元线性回归

reg = LinearRegression()

reg.fit(data_x,data_y.reshape(-1,1))

和上述一元线性回归基本没区别，但是这里要注意data_x的存放顺序。

线性倾向率

这个就直接跳过了，就是一元线性回归呗。

滑动平均

def moving_ave(x,n):

   year=x.shape

   year=np.array(year)

   year=int(year)

   n1=int((n-1)/2)

   new=[]

   for i in range(n1,year-n1):

       new.append(np.nanmean(x[i-n1:i+n1+1]))

   return np.array(new)

这是自定义的函数，x是变量的时间序列，n是滑动的尺度。

累积距平

data_mean=np.mean(data)

data_a=data-data_mean

data_lei=np.zeros(85)

data_lei[0]=data_a[0]

for i in range(1,85):

   data_lei=data_lei[i-1]+data_a

这个应该也不需要过多的解释，文章太单调了，下面放几张图吧。

最后的EOF分解在上一篇文章里提到过了，可以前往https://mp.weixin.qq.com/s/o0semESB1mxX4wQw31ECbw

合成分析

合成分析的话其实还是检验先看看教材，说白了就是检验两种情况有没有显著差异。均值是t检验，方差是F检验，公式放在下面了，理解的话应该就会写代码了（其实是我又懒了）。不知道为啥气象家园这下面的图显示不出来，但是下载下来可以显示的。

最后在整理的时候发现了当时弄了好久的用Python打开DAT文件，打算之后做一下相关内容，当然flag立在这里，什么时候能写出来也还不清楚哈哈

                               
登录/注册后可看大图

最后又到了给自己打广告的时候了，欢迎大家关注我的公众号啦，等以后关注的人多了会推一些独门干货嘻嘻

这篇文章同样给了回帖奖励！！

                               
登录/注册后可看大图

scmeiqy

普及基础方法，赞一个！

wxylndd

好厉害，楼主6666

Doubleegg

太强了大佬

灭火器

很有用的总结！

爱到花开

统计基础较差，看了程序之后学习比较方便

独孤酒见

牛人啊。挽救了类似我等这种人了。感谢牛人的付出和奉献。回头加公众号。

Mr.Sun

感谢楼主，学习了！！！

yangwawa

感谢楼主，学习了！！！

庆姜3

楼主善良，基础知识最重要