暑期数学建模导论----数据处理与拟合模型

今天是自学的第一天,让我们跟随着夏令营的节奏一起进步吧!

                                                                            声明:学习内容来源于开源学习平台Datawhale

        现在进入到正文部分啦,我现在就读的专业是数据科学大数据专业,所以对这一块的内容也算是有一定的基础,所以应该会学习的比较轻松,这篇笔记(心得)是在学习的过程中,边学边记录的,如果我的理解有误,还望各位路过的佬们批评指正!!!

1.数据科学的研究对象

大数据科学的基础框架:

数据的获取和存储:爬虫、软件定义存储、硬件存储......

数据的处理:分布式计算、并行计算、数据流知识、Hadoop、Spark......

数据分析:统计学、数据挖掘、机器学习、计算机视觉、自然语言处理......

数据的管理:现代数据库系统及其架构等内容。 

数据应用:数据可视化、数据相关软件开发、报表分析以及如何将数据挖掘得到的结果还原为解决实际问题的方案。

2.数据预处理

2.1数据预处理的引入

        数据预处理的环节不止数学建模中的异常值处理,缺失值检测和处理,还有更多对数据进行预处理的环节,就像机器学习将所有数据转换成同一纬度的StandardScaler等等;况且在现实生活中我们收集到的数据也不能直接拿来使用,比如今年我打的一个比赛---能源经济学术创意大赛,在通过国家统计局官网等平台收集数据的时候收集到的数据往往格式都是比较混乱的(这里说的混乱是没有达到程序代码可以读入的情形),这时我们就需要对数据进行预处理,能让程序更好的读入数据,其次国家给的数据一般都是经过严格的检测才会公布的,如果像一些小平台流出的数据,这时还有异常值和nan值得情况出现,所以在采集完数据后最重要的步骤便是数据预处理。

        数据和特征决定了模型效果的上限,如果一份数据非常不好,不管你的模型再好也是没用的,请牢记这段话!

        以excel为例子,大家应该都使用过excel,最上面的一行我们称之为表头,表头的每一格我们称之为字段,每个字段描述了这一列数据表示的意义;表格的体量的意思是这个表格有多少行多少列,比如这个表格的体量是3行10列;如果列数超过了行数的1/2则可以说是有点稀疏了,如果列数是行数的3倍那它就是严重稀疏的数据,此时3行10列表格的体量就是严重稀疏的;这个数据的每一列称为一个属性,有多少个属性就可以称之有多少维。

注意:稀疏还有一种定义:表格里面很多都是空白的或者全是0。

        属性有离散属性和连续属性之分。比如a属性只有{1,2,3}三个有限的取值,我们称a属性为离散属性;相反地如果b属性取值是一个连续的浮点数,这种属性称之为连续属性。                                 需要说明的是,连续属性和离散属性的处理方法截然不同。

Id即是每一行数据的标识信息,可以根据这个编号对数据进行检索、排序等操作,称之为主码

2.2 Pandas库基础

(1)Python 创建一个数据框DataFrame
import pandas as pd    #导入pandas库
import numpy as np    #导入numpy库
data={'animal':['cat','cat','snake','dog','dog','cat','snake','cat','dog','dog'],
     'age':[2.5,3,0.5,np.nan,5,2,4.5,np.nan,7,3],
     'visits':[1,3,2,3,2,3,1,1,2,1],
     'priority':['yes','yes','no','yes','no','no','no','yes','no','no']}
labels=['a','b','c','d','e','f','g','h','i','j']
df=pd.DataFrame(data)

        没有pandas和numpy的可以通过以下命令来创建:

        win+R 然后输入cmd

pip install numpy
pip install pandas

        如果是jupyter笔记本的话请复制以下代码并运行:

!pin install numpy
!pin install pandas
(2)显示该DataFrame及其数据相关的基本信息
df.describe()

        此时出来的是这份数据的基本信息(数值型),我们创建的数值型的数据列表一共只有age和visits两列,所以结果为:

agevisits
count8.00000010.000000

mean 

3.4375001.900000
std2.0077970.875595
min0.5000001.000000
25%2.3750001.000000
50%3.0000002.000000
75%4.6250002.750000
max7.0000003.000000
(3)返回DataFrame的前5列数据:
df.head(5)
animalagervisitspriority
0cat2.51yes
1cat3.03yes
2snake0.52no
3dogNaN3yes
4dog5.02no
(4)从DataFrame df 选择标签列为animal和age的列
df[['animal','age']]
(5)在3,4,8行中,选择列为animal,age的数据
df.loc[[3,4,8],['animal','age']]
(6)选择列为visits中等于3的行
df.loc[df['visits']==3,:]
(7)选择age为缺失值的行
df.loc[df['age'].isna(),:]
(8)选择animal是cat且age小于3的行
df.loc[(df['animal']=='cat')&(df['age']<3),:]
(9)选择age在2到4之间的数据(包含边界值)
df.loc[(df['animal']=='cat')&(df['age']<3),:]
(10)将f行中的age改为1.5
df.index=labels    #先将labels赋值到索引行对应的位置
df.loc[['f'],['age']]=1.5
(11)对visits列的数据求和
df['visits'].sum()
(12)计算每种animal age的平均值
df.groupby(['animal'])['age'].mean()

2.3 Pandas进阶

import pandas as pd
import numpy as np
data={'animal':['cat','cat','snake','dog','dog','cat','snake','cat','dog','dog'],
     'age':[2.5,3,0.5,np.nan,5,2,4.5,np.nan,7,3],
     'visits':[1,3,2,3,2,3,1,1,2,1],
     'priority':['yes','yes','no','yes','no','no','no','yes','no','no']}
labels=['a','b','c','d','e','f','g','h','i','j']
df=pd.DataFrame(data)
df.describe()
df.head(5)
df[['animal','age']]
df.loc[[3,4,8],['animal','age']]
df.loc[df['visits']==3,:]
df.loc[df['age'].isna(),:]
df.loc[(df['animal']=='cat')&(df['age']<3),:]
df.loc[(df['age']>=2)&(df["age"]<=4),:]
df.index=labels
df.loc[['f'],['age']]=1.5
df['visits'].sum()
df.groupby(['animal'])['age'].mean()
df=pd.DataFrame({
    'From_To':['LoNDon_paris','MAdrid_miLAN','londoN_StockhOlm','Budapest_PaRis','Brussels_londOn'],
    'FlightNumber':[10045,np.nan,10065,np.nan,10085],
    'RecentDelays':[[23,47],[],[24,43,87],[13],[67,32]],
    'Airline':['KLM(!)','<Air France> (12)','British Airways.','12.Air France','"Swiss Air"']
})
df['FlightNumber']=df['FlightNumber'].interpolate().astype(int)
#在该列中添加10055和10075填充该缺失值
#astype是转化为int类型
#由于From_To代表从地点A到地点B,因此可以将这列拆分成两列,并赋予为列From与To
temp = df['From_To'].str.split('_',expand=True)
temp.columns = ['From','To']
#将列From和To转化成只有首字母大写的形式
temp['From'] = temp['From'].str.capitalize()
temp['To'] = temp['To'].str.capitalize()
#将列From_To从df中去除,并把列From和To添加到df中
df.drop('From_To',axis=1,inplace=True)
df[['From','To']] = temp
#清楚列中的特殊字符,只留下航空公司的名字
df['Airline'] = df['Airline'].str.extract(pat=r'([a-zA-Z\s]+)',expand=False).str.strip()
#expand表示是否把series类型转化为DataFrame类型
#extract:pat=r'([a-zA-Z\s]+)'表示只保留a到z,A到Z的字符串形式的字符,+代表不止一个
#在RecentDelays列中,值已作为列表输入到DataFrame中。我们希望每个第一个值在它自己的列中,每个的第二个值在它自己的列中(新的),以此类推。
#如果没有第N个值,则该值应为NaN.
#将Series列表展开名为delays的DataFrame,重命名列delay_1,delay_2等等,并将不需要的RecentDelays列替换df为delays
delays = df['RecentDelays'].apply(pd.Series)
delays.columns = ['delay_%s' % i for i in range(1,len(delays.columns)+1)]
df = df.drop('RecentDelays',axis=1).join(delays,how='left')
#将delay_i列的控制nan都填为自身的平均值
for i in range(1,4):
    df[f'delay_{i}'] = df[f'delay_{i}'].fillna(np.mean(df[f'delay_{i}']))
df = df._append(df.loc[df['FlightNumber']==10085,:],ignore_index=True)
#在df中增加一行,值与FlightNumber=10085的行保持一致
#对df进行去重
df = df.drop_duplicates()

2.4 数据的规约

(1)min-max规约

x_{new}=\frac{x-min(x)}{max(x)-min(x)} 

(2)Z-score规约

 x_{new}=\frac{x-x\bar{}}{std(x)}

3.常见的统计分析模型

3.1回归分析与分类分析

        回归分析和分类分析都是一种基于统计模型的统计分析方法。它们都研究因变量(被解释变量)与自变量(解释变量)之间存在的潜在关系,并通过统计模型的形式将这些潜在关系进行显示的表达。

        不同的是,回归分析中因变量是连续变量,如工资、销售额;而分类分析中因变量是属性变量,如判断邮件“是or否”垃圾邮件。     下面试通过几个简单的例子来说明回归分析与分类分析在数学建模中的应用。

(1)回归分析

        通过gpa1的数据来分析大学成绩GPA是否与大学入学考试的成绩和高考成绩有关,并且这两个哪个相关程度更高。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import seaborn as sns
from IPython.display import display
import statsmodels.api as sm
#加载数据
gpa1=pd.read_stata('./data/gpa1.dta')
gpa1
#在数据集中提取自变量(ACT为入学考试成绩,hsGPA为高考成绩)
X1=gpa1.ACT
X2=gpa1[['ACT','hsGPA']]
X1
X2
#提取因变量
y=gpa1.colGPA
y
#为自变量增添截距项
X1=sm.add_constant(X1)
X2=sm.add_constant(X2)
display(X2)
#在 Python 中,display 是一个用于在 Jupyter Notebook 和 IPython 环境中显示对象的函数。它与 print 的作用类似,但功能更为强大。
#display 可以用来显示更丰富的内容,比如 HTML、Markdown、数据表格、图像等。
#拟合两个模型
gpa_lm1 = sm.OLS(y,X1).fit()
gpa_lm2 = sm.OLS(y,X2).fit()
p1=pd.DataFrame(gpa_lm1.pvalues,columns=['pvalue'])
c1=pd.DataFrame(gpa_lm1.params,columns=['params'])
p2=pd.DataFrame(gpa_lm2.pvalues,columns=['pvalue'])
c2=pd.DataFrame(gpa_lm2.params,columns=['params'])
display(c1.join(p1,how='right'))
display(c2.join(p2,how='right'))

此时得到结果为 :

1.gpa_lm1模型:

paramspvalue
const2.4029798.798591e-16
ACT0.0270641.389927e-02

2.gpa_lm2模型:

paramspvalue
const1.2863280.000238
ACT0.0094260.383297
hsGPA0.4534560.000005

        从模型I可以看到,ACT变量的pvalue为0.01<0.05,说明大学入学考试成绩ACT能显著影响大学GPA; 从模型II可以看到,高考成绩haGPA的pvalue为0.000005远小于0.05,说明高考成绩haGPA显著影响大学成绩GPA,但大学入学考试成绩ACT的pvalue为0.38大于0.05,不能说明大学入学考试成绩ACT显著影响大学成绩GPA。

        为什么会出现这样矛盾的结论呢?原因有可能是:在没有加入高考成绩haGPA 时,大学入学考试成绩ACT确实能显著影响大学成绩GPA,但是加入高考成绩haGPA后,有可能高考成 绩haGPA高的同学大学入学考试也高分,这两者的相关性较高,显得“大学入学考试成绩ACT并不是那么重要了。

        从模型II的系数params数值可以看到,当高考成绩haGPA显著影响大学成绩GPA,params为0.45代表每1单位的高考成绩的增加将导致大学GPA平均增加0.45个单位

(2)分类模型

        分类分析:ST是我国股市特有的一项保护投资者利于的决策,当上市公司因财务状况不佳导致投 资者难以预测其前景时,交易所会标记该公司股票为ST,并进行一系列限制措施。

        我们想研究被ST的公司其背后的因素,并尝试通过利用公司的财务指标提前预测某上市公司在未来是否会被ST,是一件很有意义的举措。而在这项任务中,因变量就是公司是否会被ST。该例中自变量是一些财务指标,如ARA、ASSET等,它们具体的意义在此不做赘述。

#读取数据
ST = pd.read_csv('./data/ST.csv')
ST.head()
#这段代码使用了 statsmodels 库来执行一个逻辑回归(Logistic Regression)分析,并打印出回归结果的摘要。
#ST:因变量,即你要预测的目标变量。
#ARA, ASSET, ATO, ROA, GROWTH, LEV, SHARE:自变量(或解释变量),即用于预测因变量的特征。
st_logit=sm.formula.logit('ST~ARA+ASSET+ATO+ROA+GROWTH+LEV+SHARE',data=ST).fit()
print(st_logit.summary())

结果如下图所示:

其中pvalue为第四列标记为(P>|z| )的数值。

        与回归分析类似,一般认为,p值小于0.05的自变量是显著的,也就是说,有统计学的证据表明该变量会影响因变量为1的概率(即顾客购买杂志的概率)。

3.2 假设检验

        学过数理统计的人应该都知道H1,H2;假设检验是统计分析中重要的一环,它贯穿着统计分析的所有环节。在数学建模中,我们既能通过假设检验对数据进行探索性的信息挖掘,以给予我们选择模型充分且客观的依据;又能在建模完成后,通过特定的假设检验验证模型的有效性。因此,储备基本的假设检验知识,学会根据数据特点、任务需求选择相对应的假设检验十分重要。

        假设检验的本质是:根据样本信息与已知信息,对一个描述总体性质的命题进行“是或否”的检验与回答。即:假设检验验证的不是样本本身的性质,而是样本所在总体的性质。

        假设检验大体上可分为两种:参数假设检验非参数假设检验。 

        顾名思义,若假设是关于总体的一个参数或者参数 的集合,则该假设检验就是参数假设检验;若假设不能用一个参数集合来表示,则该假设检验就是非参数假设检验,一个典型就是正态性检验。

(1)正态性检验

        若数据是正态分布的,我们应该使用参数检验,此时对数据进行正态性检验就非常有必要了。在这里,有三种方法判断数据的正态性:可视化判断-正态分布概率图;Shapiro-Wilk检验;D'Agostino's K-squared检验。

#生成1000个服从正态分布的数据
data_norm = stats.norm.rvs(loc=10,scale=10,size=1000)
#rvs(loc,scale,size):生成服从指定分布的随机数,loc:期望;scale:标准差;size:数据个数
#生成1000个服从卡方分布的数据
data_chi = stats.chi2.rvs(2,3,size=1000)
#画出两个概率图
# 创建一个图形对象,设置大小为12x6英寸
fig = plt.figure(figsize=(12, 6))
# 添加第一个子图,位置为1行2列中的第1个
ax1 = fig.add_subplot(1, 2, 1)
# 使用probplot绘制服从正态分布数据的概率图,并将图绘制在ax1上
plot1 = stats.probplot(data_norm, plot=ax1)
# 添加第二个子图,位置为1行2列中的第2个
ax2 = fig.add_subplot(1, 2, 2)
# 使用probplot绘制服从卡方分布数据的概率图,并将图绘制在ax2上
plot2 = stats.probplot(data_chi, plot=ax2)
# 显示绘制的图形
plt.show()

结果如下图所示:

        可以看到,正态性数据在正态分布概率图中十分接近直线y=x,而卡方分布数据则几乎完全偏离了该直线。可见,概率图可视化确实可以起到初步判断数据正态性的作用。实际应用中,由于数据的复杂性, 仅使用一种方法判断正态性有可能产生一定的误差,因此我们通常同时使用多种方法进行判断。

        如果不同方法得出的结论不同,此时就需要仔细观察数据的特征,寻找结果不一致的原因。如:若Shapiro Wilk test(夏皮罗-威尔克检验)显著(非正态),D'Agostino's K-squared test不显著(正态),则有可能是因为样本量较大,或者样本中存在重复值现象,如果事实确实如此,那么我们就应该采纳D'Agostino's K-squared test的结论而非Shapiro-Wilk test的结论。在这里,我们假设数据服从正态分布,,p<0.05则拒绝原假设(即p<0.05则说明该数据不服从正态分布,如果p>=0.05则没有证据说明数据不服从正态分布)

#接下来,在python中定义一个函数,
#将概率图,Shapiro Wilk test,D'Agostino's K-squared test结合在一起
data_small = stats.norm.rvs(0,1,size=30)    #小样本正态性数据集
data_large = stats.norm.rvs(0,1,size=6000)     #大样本正态性数据集
from statsmodels.stats.diagnostic import lilliefors
from typing import List
'''
输入参数----------
data : numpy数组或者pandas.Series
show_flag : 是否显示概率图
Returns------
两种检验的p值;概率图
'''
#函数 check_normality 的返回类型标注为 List[float]
#即函数返回一个包含两种正态性检验的 p 值的列表。
def check_normality(data: np.ndarray, show_flag: bool=True) -> List[float]:
    if show_flag:
         _ = stats.probplot(data, plot=plt)
         plt.show()
    pVals = pd.Series(dtype='float64')
    # D'Agostino's K-squared test
    _, pVals['Omnibus'] = stats.normaltest(data) 
    # Shapiro-Wilk test
    _, pVals['Shapiro-Wilk'] = stats.shapiro(data)
    print(f'数据量为{len(data)}的数据集正态性假设检验的结果 : ----------------')
    print(pVals)
 check_normality(data_small,show_flag=True)
 check_normality(data_large,show_flag=False)

 

         当样本量大于5000会出现警告:

        经验证,两正态分布小样本和大样本的数据均服从正态分布,p值均大于0.05!

(2)单组样本均值假定的检验

        必胜中学里,陈老师班结束了一次英语考试。由于班级人数较多,短时间内很难完成批改与统计,陈老师又很想知道此次班级平均分与级长定的班级均分137的目标是否有显著区别,于是他随机抽取了已经改好的10名同学的英语成绩: 136,136,134,136,131,133,142,145,137,140

        问:陈老师可以认为此次班级平均分与级长定的班级均分137的目标没有显著区别吗?

        很明显,这就是一个典型的单组样本均值假定的检验,比较的是这个样本(10个同学的英语成绩)所代 表的总体均值(班级英语成绩均值)是否与参考值137相等。        

        那么对于这类问题,我们有两种检验可以使用:单样本t检验wilcoxon检验

data=np.array([136,136,134,136,131,133,142,145,137,140])
#定义一个单组样本均值检验函数,使它可以同时输出t检验与mannwhitneyu检验的p值
'''
输入参数
    ---------
    data : numpy数组或者pandas.Series
    checkvalue : 想要比较的均值
    significance : 显著性水平
    alternative : 检验类型,这取决于我们备择假设的符号:
    two-sided为双侧检验、greater为右侧检验、
    less为左侧检验
输出
------
在两种检验下的p值
在显著性水平下是否拒绝原假设
'''
def check_mean(data,checkvalue,significance=0.05,alternative='two-sided'):
    pVal=pd.Series(dtype='float64')
    # 正态性数据检验-t检验
    _, pVal['t-test'] = stats.ttest_1samp(data, checkvalue,alternative=alternative)
    print('t-test------------------------')
    if pVal['t-test'] < significance:
         print(('目标值{0:4.2f}在显著性水平{1:}下不等于样本均值(p={2:5.3f}).'.format(checkvalue,significance,pVal['t-test'])))
    else:
         print(('目标值{0:4.2f}在显著性水平{1:}下无法拒绝等于样本均值的假设.(p={2:5.3f})'.format(checkvalue,significance,pVal['t-test'])))
    # 非正态性数据检验-wilcoxon检验
    _, pVal['wilcoxon'] = stats.wilcoxon(data-checkvalue,alternative=alternative)
    print('wilcoxon------------------------')    
    if pVal['wilcoxon'] < significance:
         print(('目标值{0:4.2f}在显著性水平{1:}下不等于样本均值(p={2:5.3f}).'.format(checkvalue,significance,pVal['wilcoxon'])))
    else:
         print(('目标值{0:4.2f}在显著性水平{1:}下无法拒绝等于样本均值的假设.(p={2:5.3f})'.format(checkvalue,significance,pVal['wilcoxon'])))
    return pVal
 check_mean(data,137,0.05)

注明:

_, pVal['t-test'] = stats.ttest_1samp(data, checkvalue, alternative=alternative)

        这串代码单独拿出来学习一下 ,scipy.stats.ttest_1samp 函数返回两个值,第一个值是 t 统计量(t-statistic),第二个值是 p 值(p-value),假设你在计算 t 检验时只关心 p 值:

t_statistic, p_value = stats.ttest_1samp(data, checkvalue, alternative=alternative)

         但如果你只关心 p 值,可以这样写:

_, p_value = stats.ttest_1samp(data, checkvalue, alternative=alternative)

        这样,_ 表示你不打算使用 t_statistic,只使用 p_value。 _ 是一个常见的 Python 编程惯例,用于表示不打算使用的变量。它有助于提高代码的清晰度和可维护性,同时避免不必要的变量创建。在你的例子中,_ 用于忽略 t 统计量,只关注并存储 p 值。

结果为: 

(3) 两组样本的均值相等性检验

        在进行两组之间的均值比较之前,我们需要进行一个重要的判断:这两组样本之间是否独立呢?这个问题的答案将决定着我们使用何种类型的检验。

        陈老师在年级有一个竞争对手:王老师。王老师在他们班级在同一时间也举行了一次英语考试,且两个班用的是同一份卷子,因此两个班的英语成绩就具有比较的意义了。和陈老师一样,王老师也来不及批改和统计他们班级的英语成绩,不过他手头上也有12份已经改好的试卷,成绩分别为: 134,136,135,145,147,140,142,137,139,140,141,135

        问:我们可以认为两个班级的均分是相等的吗?

        这是一个非常典型的双独立样本的均值检验。

        值得注意的是,这里的独立指的是抽样意义上的独立,而 不是统计意义的独立,什么意思呢?即我们只需要保证这两个样本在选取的时候是“现实上”的互不影响 就可以了。至于两者在数值上是否独立(通过独立性检验判断的独立性)我们并不关心。对于抽样意义 上的独立的解释,各种教材众说纷纭。        

        在这里我们选取一种容易理解的说法:两个样本中,一个样本中 的受试不能影响另一个样本中的受试。

         在本例中,两个班级的授课老师不同,因此两个班学生的英语学习不会受到同一个老师的影响;两个班 级考试同时进行,意味着不存在时间差让先考完的同学给后考完的同学通风报信(进而影响受试)。因此,我们可以认为这是两个独立的样本。 事实上,两个样本是否在抽样意义上独立,是没有固定答案的,因为在很多情况下,我们既不能保证两个样本间完全独立,也很难判断出两者是否存在相关性。我们只能在抽样的时候使用更加科学的抽样方法,尽可能地避免样本间的相互影响。

        若两个样本的总体都服从正态分布,那么我们可以使用双样本t检验。如果不服从正态分布,则可以使用Mannwhitneyu秩和检验,Mannwhitneyu秩和检验是一种非参数检验。故在下方的代码中采用t检验与Mannwhitneyu秩和检验两种检验的方式,在采用t检验的时候需要先使stats.levene 函数进行 Levene 检验,以确定两组样本的方差是否相等。以便确认t检验过程中equal_var=False这个参数是传入False还是True。

""" 
输入参数----------
group1/2 : 用于比较的两组数据
significance : 显著性水平
alternative : 检验类型,这取决于我们备择假设的符号:two-sided为双侧检验、greater为右侧
检验、less为左侧检验
输出------
在两种检验下的p值
在显著性水平下是否拒绝原假设
"""
# 定义一个单组样本均值检验函数,使它可以同时输出t检验与mannwhitneyu检验的p值
def unpaired_data(group1:np.ndarray,group2:np.ndarray,significance=0.05,alternative='two-sided'):
    pVal=pd.Series(dtype='float64')
    # 先进行两组数据的方差齐性检验
    #使用 stats.levene 函数进行 Levene 检验,以确定两组样本的方差是否相等。
    _,pVal['levene']=stats.levene(group1,group2)
    # t检验-若数据服从正态分布
    if pVal['levene']<significance:
         print('在显著性水平{0:}下,两组样本的方差不相等(p={1:.4f}),因此需要使用方差不等的t检验'.format(significance,pVal['levene']))
         print('------------------------------------')
         _, pVal['t-test'] = stats.ttest_ind(group1, group2,equal_var=False,alternative=alternative) # 因为方差不相等,因此是False
         print('t检验p值:{}'.format(pVal['t']))
    else:
         print('在显著性水平{0:}下,不能拒绝两组样本方差相等的假设(p={1:.4f}),因此需要使用方差相等的t检验'.format(significance,pVal['levene']))
         print('------------------------------------')
         _, pVal['t-test'] = stats.ttest_ind(group1,group2,equal_var=True,alternative=alternative) # 因为方差相等,因此是True     
         print('t检验p值:{:.3f}'.format(pVal['t-test']))   
    # mannwhitneyu检验-数据不服从正态检验
    #使用 stats.mannwhitneyu 函数进行 Mann-Whitney U 检验,用于比较两组样本的分布。
    _, pVal['mannwhitneyu'] = stats.mannwhitneyu(group1,group2,alternative=alternative)
    print('Mann-Whitney检验p值:{:.3f}'.format(pVal['mannwhitneyu']))
    # --- >>> STOP stats <<< --
    # 两组样本均值的散点图可视化
    print('------------------------------------')
    print('两组样本均值的散点图可视化')
    plt.plot(group1, 'bx', label='group1')
    plt.plot(group2, 'ro', label='group2')
    plt.legend(loc=0)
    plt.show()
    return pVal
# 陈老师班
group1=data
# 王老师班
group2=np.array([134,136,135,145,147,140,142,137,139,140,141,135])
unpaired_data(group1,group2)

结果如下图所示: 

        两个检验的结果都表明两个班级的平均分没有显著差异。(即没有证据不能拒绝原假设)

(4)成对检验

        在进行两组间均值比较的时候,有一种特殊情况——两个样本“故意”不独立。这种情况多出现在两个样本分别为同一个受试个体的情况。对于同一个受试个体下的不同时间的受试结果,那么说这两个样本是“成对”的,是彼此紧密相连的。对这样两个样本进行均值比较检验,就是成对检验。

        为了让大家更好地理解“成对”的概念,我们再把王老师的故事讲长一点。得知了王老师班的考试均值与自己班的没有显著差异,陈老师非常生气:“我堂堂优秀教师,居然不能和隔壁老王拉开差距,岂有此理!”,于是陈老师开始了为期一周的魔鬼训练,并在一个星期后又进行了一次全班的测验。陈老师依旧抽取了上次十位同学的成绩,依次为:                                 139,141,137,136,135,132,141,148,145,139

        问:这次班级均分与上次是否存在显著差异呢?

        显然,两个样本分别为相同的同学前后两次的考试成绩,是非常典型的成对数据,因此我们可以使用成对检验。成对检验也分为两种:
        

 1.若总体服从正态分布,则使用成对t检验;

 2.若总体不服从正态分布,则使用成对wilcoxon秩和检验。

def paired_data(group1:np.ndarray,group2:np.ndarray,significance,alternative='two-sided'):
    pVal=pd.Series(dtype='float64')
    # 配对t检验-样本服从正态分布
    _, pVal['t-test'] = stats.ttest_1samp(post - pre,0,alternative=alternative)
    print('t-test------------------------')
    if pVal['t-test'] < significance:
         print(('在显著性水平{0:}下,两组配对样本的均值不相等(p={1:5.3f}).'.format(significance,pVal['t-test'])))
    else:
         print(('在显著性水平{0:}下无法拒绝等于样本均值的假设.(p={1:5.3f})'.format(significance,pVal['t-test'])))    
    # wilcoxon秩和检验
    _, pVal['wilcoxon'] = stats.wilcoxon(group1,group2,mode='approx',alternative=alternative)
    print('wilcoxon------------------------')
    if pVal['wilcoxon'] < significance:
         print(('在显著性水平{0:}下,两组配对样本的均值不相等(p={1:5.3f}).'.format(significance,pVal['wilcoxon'])))
    else:
         print(('在显著性水平{0:}下无法拒绝等于样本均值的假设.(p={1:5.3f})'.format(significance,pVal['wilcoxon'])))    
    return pVal
# 第一次测验
pre=data
# 第二次测验
post=np.array([139,141,137,136,135,132,141,148,145,139])
paired_data(pre,post,0.05)

结果如下图所示:

         结论:两种检验都显示两次测验的均分在显著性水平0.05下有显著差异。

        个人理解:

在双边检验中,p 值是为了检测任一方向上的显著性。如果我们有理由认为第二次测验的均值显著高于第一次测验,可以转换成单边检验,双边检验的 p 值是 0.039370/0.048997,单边检验的 p 值会更小,进一步支持结论成立,由于两种检验都显示显著差异,再结合实际数据的观察和上下文,可以推测第二次测验的均值显著高于第一次测验的均值。

(5)方差分析-多组样本间的均值相等性检验

        多组样本间的均值相等性检验可以转化成观察-各样本的样本均值的-差异程度,如果各样本的均值差异程度很大的话,那么它们的总体均值也有很大可能存在差异。(样本间均值的差异程度。)

        各样本的样本内差异程度也很重要,在相同的样本间的差异程度下,样本内的差异程度越大,各总体间均值存在的差异的可能性就越小。

        举个栗子:样本内的差异程度很大了,就代表偶然性很大,就很难去判断两个总体间不相等的均值是否真的不相等。

        比如,我考试均分是91,我的同桌考试均分是89,假如取一个非常极端的情况就是,我们考试分数的标准差都为0,那我一直都是91,我同桌一直都是89,那现在很容易判断出我们两个的均分存在明显的差异。

        但如果此时标准差为6(方差为36的话)我们的成绩都很不稳定,波动都很大,在36的方差下,2分的均值差没有什么说服力了。

        因此,想要综合这两个评判指标,解决的办法就是两者相除(大/小),即样本间均值的差异程度除以样本内的差异程度。这也是方差分析最根本的思想。

         

对象(非正态性的数据)大样本小样本
方差分析稳健存在误差
kruskalwallis检验(非参数检验)X

        下面对altman_910.txt数据集进行方差分析,该数据记录了3组心脏搭桥病人给予不同水平的一氧化氮通气下,他们的红细胞内叶酸水平。

        注:三组样本都是正态性样本。

data=np.genfromtxt('./data/altman_910.txt',delimiter=',')
data
group1 = data[data[:,1]==1,0]    #先data的第二列的数据为1的所有数据提取出来然后取第一列
group2 = data[data[:,1]==2,0]
group3 = data[data[:,1]==3,0]
group1
from typing import Tuple
def anova_oneway()->Tuple[float,float]:
    pVal=pd.Series(dtype='float64')
    # 先做方差齐性检验
    _,pVal['levene'] = stats.levene(group1, group2, group3)
    if pVal['levene']<0.05: #这里假设显著性水平为0.05
         print('警告: 方差齐性检验的p值小于0.05: p={},方差分析结果在小样本下可能不准确'.format(pVal['levene']))
         print('-------------------------------')
     # 单因素方差分析-假设样本服从正态分布
    _, pVal['anova_oneway_normal'] = stats.f_oneway(group1, group2, group3) # 在这里输入待分析的数据
    print('若样本服从正态分布,单因素方差分析的p值为{}'.format(pVal['anova_oneway_normal']))
    if pVal['anova_oneway_normal'] < 0.05:
         print('检验在0.05的显著性水平下显著,多组样本中至少存在一组样本均值与其它样本的均值不相等。')
    print('---------------------------------')
     # 单因素方差分析-假设样本不服从正态分布
    _, pVal['anova_oneway_notnormal'] = stats.mstats.kruskalwallis(group1,group2, group3) # 在这里输入待分析的数据
    print('若样本不服从正态分布,单因素方差分析的p值为{}'.format(pVal['anova_oneway_notnormal']))
    if pVal['anova_oneway_notnormal'] < 0.05:
         print('检验在0.05的显著性水平下显著,多组样本中至少存在一组样本均值与其它样本的均值不相等。')    
    return pVal
anova_oneway()

3.3 随机过程与随机模拟

        随机过程就是在随机变量的基础上加入了时间维度(值得注意的是,时间维度不是随机变量,只是普通变量),随机过程与概率论有很强的关联性,又有不同。随机过程是一门专业性很强的学科,统计学专业有开,很不幸,我这个专业没开。。。

        但是话又说回来我可以自学!!!

案例1:为了改善道路的路面情况(道路经常维修),需要统计一天中有多少车辆经过,因为每天的车辆数是随机的,一般来说有两种技术解决这个问题:

~在道路附近安装一个计数器或安排一个技术人员,在一段长时间的天数(如365天)每天24h统计 通过道路的车辆数。

~使用仿真技术大致模拟下道路口的场景,得出一个近似可用的仿真统计指标

        由于方案1的过程太过繁琐,虽然能得到准确的结果,但是需要耗费大量的时间和人力。因此,我们选择第二个:使用仿真技术大致模拟下道路口的场景,得出一个近似可用的仿真统计指标

        通过一周的简单调查,得到每天的每 个小时平均车辆数:[30, 20, 10, 6, 8, 20, 40, 100, 250, 200, 100, 65, 100, 120, 100, 120, 200, 220, 240, 180, 150, 100, 50, 40],通过利用平均车辆数进行仿真。

#模拟仿真研究该道路口一天平均有多少车经过
import simpy
class Road_Crossing:
    def __init__(self,env):
        self.road_crossing_container =simpy.Container(env,capacity=1e9,init=0)
def come_across(env, road_crossing, lmd):
    while True:
        body_time = np.random.exponential(1.0/(lmd/60))  # 经过指数分布的时间后,泊松过程记录数+1
        yield env.timeout(body_time)  # 经过body_time个时间
        yield road_crossing.road_crossing_container.put(1)
hours = 24  # 一天24h
minutes = 60  # 一个小时60min
days = 3   # 模拟3天
lmd_ls = [30, 20, 10, 6, 8, 20, 40, 100, 250, 200, 100, 65, 100, 120, 100, 120, 200, 220, 240, 180, 150, 100, 50, 40]   # 每隔小时平均通过车辆数
car_sum = []  # 存储每一天的通过路口的车辆数之和
print('仿真开始:')
for day in range(days):
    day_car_sum = 0   # 记录每天的通过车辆数之和
    for hour, lmd in enumerate(lmd_ls):
         env = simpy.Environment()
         road_crossing = Road_Crossing(env)
         come_across_process = env.process(come_across(env, road_crossing, lmd))
         env.run(until = 60)  # 每次仿真60min
         if hour % 4 == 0:
             print("第"+str(day+1)+"天,第"+str(hour+1)+"时的车辆数:", road_crossing.road_crossing_container.level)
         day_car_sum += road_crossing.road_crossing_container.level
    car_sum.append(day_car_sum)
print("每天通过交通路口的的车辆数之和为:", car_sum)

案例2: 假设每个小时进入商店的平均人数为:[10, 5, 3, 6, 8, 10, 20, 40, 100, 80, 40, 50, 100, 120, 30, 30, 60, 80, 100, 150, 70, 20, 20, 10],每位顾客的平均花费为:10元(大约一份早餐吧),请问每天商店的营业额是多少?(仿真复合泊松过程)

# 模拟仿真研究该商店一天的营业额
import simpy
class Store_Money:
    def __init__(self, env):
         self.store_money_container = simpy.Container(env, capacity = 1e8, init = 0)
def buy(env, store_money, lmd, avg_money):
    while True:
        body_time = np.random.exponential(1.0/(lmd/60))  # 经过指数分布的时间后,泊松过程记录数+1
        yield env.timeout(body_time) 
        money = np.random.poisson(lam=avg_money)
        yield store_money.store_money_container.put(money)
hours = 24  # 一天24h
minutes = 60  # 一个小时60min
days = 3   # 模拟3天
avg_money = 10
lmd_ls = [10, 5, 3, 6, 8, 10, 20, 40, 100, 80, 40, 50, 100, 120, 30, 30, 60, 80, 100, 150, 70, 20, 20, 10]   # 每个小时平均进入商店的人数
money_sum = []  # 存储每一天的商店营业额总和
print('仿真开始:')
for day in range(days):
    day_money_sum = 0   # 记录每天的营业额之和
    for hour, lmd in enumerate(lmd_ls):
         env = simpy.Environment()
         store_money = Store_Money(env)
         store_money_process = env.process(buy(env, store_money, lmd, avg_money))
         env.run(until = 60)  # 每次仿真60min
         if hour % 4 == 0:
             print("第"+str(day+1)+"天,第"+str(hour+1)+"时的营业额:",store_money.store_money_container.level)
         day_money_sum += store_money.store_money_container.level
    money_sum.append(day_money_sum)
print("每天商店的的营业额之和为:", money_sum)

案例3: 

        艾滋病发展过程分为四个阶段(状态),急性感染期(状态 1)、无症状期(状态 2), 艾滋 病前期(状态 3), 典型艾滋病期(状态 4)。

        艾滋病发展过程基本上是一个不可逆的过程,即:状态1 -> 状态2 -> 状态3 -> 状态4。现在收集某地600例艾滋病防控数据,得到以下表格:

转移后
转移前状态1234总计
110625380
21409357290
318040220
41010
总计600

        现在,我们希望计算若一个人此时是无症状期(状态2)在10次转移之后,这个人的各状态的概率是多少?

        数学过程:

1.状态转移矩阵

        给定的状态转移矩阵 P 描述了从一个状态转移到另一个状态的概率。假设有 4 个状态(状态 1 到 状态 4),则状态转移矩阵 P 是一个 4×4的矩阵:

 2.初始状态分布

        初始状态分布 p0 表示在开始时每个状态的概率。假设初始时刻病人在无症状期(状态 2),即:

 3.状态分布计算

        在经过 N 次转移之后的状态分布p_{N}可以通过以下公式计算:

        其中P^{N}表示状态转移矩阵P的N次幂

        按照数学公式推导实现编程过程:

import numpy as np
# 研究无症状期病人在10期转移后的状态分布
#numpy.matmul 函数返回两个数组的矩阵乘积。
def get_years_dist(p0, P, N):
    P1 = P
    for i in range(N):
        P1 = np.matmul(P1, P)
    return np.matmul(p0, P1)
# 初始状态分布:在无症状期
p0 = np.array([0, 1, 0, 0])
# 状态转移矩阵
P = np.array([
    [10.0/80, 62.0/80, 5.0/80, 3.0/80],
    [0, 140.0/290, 93.0/290, 57.0/290],
    [0, 0, 180.0/220, 40.0/220],
    [0, 0, 0, 1]
])
# 转移次数
N = 10
# 计算经过 N 次转移后的状态分布
result = get_years_dist(p0, P, N)
# 输出结果
print(str(N) + "期转移后,状态分布为:", np.round(result, 4))
#np.round 将结果四舍五入到小数点后四位。

        得到最终结果:

4.数据可视化 

        首先通过一个简单的例子来说明什么才算是一个好的数据可视化图表:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.sans-serif']=['SimHei','Songti SC','STFangsong']
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
from plotnine import *

# 创建 DataFrame
df = pd.DataFrame({
    'variable': ['gender', 'gender', 'age', 'age', 'age', 'income', 'income', 'income', 'income'],
    'category': ['Female', 'Male', '1-24', '25-54', '55+', 'Lo', 'Lo-Med', 'Med', 'High'],
    'value': [60, 40, 50, 30, 20, 10, 25, 25, 40],
})

# 设置分类变量
df['variable'] = pd.Categorical(df['variable'], categories=['gender', 'age', 'income'])
df['category'] = pd.Categorical(df['category'], categories=df['category'])

# 创建堆叠柱状图
plot = (
    ggplot(df, aes(x='variable', y='value', fill='category')) + geom_col()
)
#ggsave(filename='stacked_bar_chart.png', plot=plot, dpi=300)   这个方法不得行我也不知道为啥
fig = plot.draw()
fig.savefig('示例图.png',dpi=600)
plot.show()

#柱状图
plot1=(
 ggplot(df, aes(x='variable', y='value', fill='category'))+
 geom_col(stat='identity', position='dodge')
)
fig1 = plot1.draw()
fig1.savefig('./示例图柱状图.png',dpi=800)
plot.show()

        大家现在可以先忽略可视化的代码部分,着重观察两个图表在表达信息上的差异。以上两个图其实都是表示:在gender、age和income上,不同类别的占比,但是第一个图明显比第二个图表达的信息更多,因为第一个图能表达每个类别的大致占比情况,而第二个图只能比较不同类别下的数量情况。

4.1 Python三大数据可视化工具库的简介

(1) Matplotlib

        Matplotlib是一个面向对象的绘图工具库,pyplot是Matplotlib最常用的一个绘图接口,调用方法如下:

import matplotlib.pyplot as plt
#创建一个图形对象,并设置图形对象的大小(可以想象在白纸中添加一个图,并设置图的大小)
plt.figure(figsize=(12,8))
#在纸上的坐标系中绘制散点
plt.scatter(x=x,y=y)
#设置x轴的标签label
plt.xlabel('x')
#设置y轴的标签label
plt.ylabel('y')
#设置图表的标题
plt.title('y=sin(x)')
#展示图标
plt.show()

举个例子:

#创建数据
x = np.linspace(-2*np.pi,2*np.pi,100)
y = np.sin(x)
plt.figure(figsize=(8,6))
plt.scatter(x=x,y=y)
plt.xlabel('x')
plt.ylabel('y')
plt.title('y=sin(x)')
plt.show()

        这是最基础的,最简单的散点图,但是以上的方法没有体现Matplotlib面向对象的特性,下面举一个例子来体会这个库面向对象绘图的特性:

# 准备数据
x = np.linspace(-2*np.pi, 2*np.pi, 100)
y1 = np.sin(x)
y2 = np.cos(x)
# 绘制第一个图:
fig1 = plt.figure(figsize=(6,4), num='first')
fig1.suptitle('y = sin(x)')
plt.scatter(x=x,y=y1)
plt.xlabel('x')
plt.ylabel('y')
#绘制第二个图:
fig2 = plt.figure(figsize=(6,4), num='second')
fig2.suptitle('y = cos(x)')
plt.scatter(x=x, y=y2)
plt.xlabel('x')
plt.ylabel('y')
plt.show()

        根据这个例子可以理解到Matplotlib的每一句代码都是往纸上添加一个绘图特征的这句话。

        正因为如此,Matplotlib优缺点也非常显而易见,优点是非常简单,缺点是绘制复杂的图表时候要一步步绘制,代码量十分巨大。

         Seaborn是 在Matplotlib的基础上的再次封装,是对Matplotlib绘制统计图表的简化。下面,我们一起看看 Seaborn的基本绘图逻辑。

(2)Seaborn

        Seaborn主要用于统计分析绘图,它是基于Matplotlib进行了更高级的API封装。

        Seaborn比Matplotlib更加易用,尤其在统计图标的绘制上,因为它避免了Matplotlib中多种参数的设置。(可以将Seaborn看作是Matplotlib的补充)

        下面,我们使用一个简单的例子,来看看使用Seaborn绘图与使用Matplotlib绘图之间的代码有什么不同:

# 准备数据
x = np.linspace(-10, 10, 100)
y = 2 * x + 1 + np.random.randn(100)
df = pd.DataFrame({'x':x, 'y':y})
# 使用Seaborn绘制带有拟合直线效果的散点图
sns.lmplot(x="x", y="y", data=df)
# 显示图表
plt.show()

(3) Plotnine

        ggplot2奠定了R语言中的数据可视化在R语言数据科学的统治地位,R语言的数据可视化是大一统的(提到R语言数据可视化首先想到的就是ggplot2)

        数据可视化一直是Python的短板,即使拥有Matplotlib和Seaborn等数据可视化包,也无法与R语言的ggplot2相媲美。原因在于当绘制复杂的图表时,Matplotlib和Seaborn由于每一句代码都是往纸上添加一个绘图特征的特性而需要大量的代码语句。而Plotnine可以说是ggplot2在Python上的移植版,使用的基本语法和R语言的ggplot2语法一模一样,使得Python的数据可视化能力大幅度提升。

        问题又来了,为什么ggplot2和Plotnine会更适合数据可视化呢?

        原因类似于PhotoShop绘图和PPT绘图的区别,与PPT一笔一画的绘图方式不同的是,PhotoShop绘图采用了‘图层’的概念,每一句代码都相当于往图像中添加一个图层,一个图层就是一类绘图动作,这无疑给数据可视化工作大大减负,同时更符合绘图者的认知。

        下面,我们通过一个案例来看看Plotnine的图层概念以及Plotnine的基本绘图逻辑:

from plotnine import *    #讲Plotnine所有模块引入
from plotnine.data import *    #引入Plotnine自带的数据集
mpg.head()

        mpg数据集记录了美国1999年和2008年部分汽车的制造厂商,型号,类别,驱动程序和耗油量。其中,manufacturer 生产厂家 ;model 型号类型 ;year 生产年份 ;cty 和 hwy分别记录城市和高速公路驾驶耗油量 ;cyl 气缸数; displ 表示发动机排量 ;drv 表示驱动系统:前轮驱动、(f),后轮驱动®和四轮驱动(4) ;class 表示车辆类型,如双座汽车,suv,小型汽车 ;fl (fuel type),燃料类型。

#绘制汽车在不同驱动系统下,发动机排量与耗油量的关系
p1=(
    ggplot(mpg,aes(x='displ',y='hwy',color='drv'))    
    #设置数据映射图层,数据集使用mpg,x数据使用mpg['displ'],y数据使用mpg['hwy'],颜色映射使用mpg['drv']
    + geom_point()    #绘制散点图图层
    + geom_smooth(method='lm')      #绘制平滑先图层
    + labs(x='displacement',y='horsepower')    #绘制x,y标签图层
)
#print(p1)    #展示p1图像
p1.show()     #这也是展示图像的方法

        从上面的案例,我们可以看到Plotnine的绘图逻辑是:一句话一个图层。因此,在Plotnine中少量的代码就能画图非常漂亮的图表,而且可以画出很多很复杂的图表,就类似于PhotoShp能轻松画出十分复杂的图但是PPT需要大量时间也不一定能达到同样的效果。

        当然这个库也是有缺陷的,Plotnine虽然具有ggplot2的图层特性,但是由于开发时间较晚,目前这个工具包还有一些缺陷,其中最大的缺陷就是:没有实现除了直角坐标以外的坐标体系,如:极坐标。因此,Plotnine无法绘制类似于饼图、环图等图表。

        为了解决这个问题,在绘制直角坐标系的图表时,我们可以使用Plotnine进行绘制,当涉及极坐标图表时,我们使用Matplotlib和Seaborn进行绘制。有趣的是,Matplotlib具有ggplot风格,可以通过设置ggplot风格绘制具有ggplot 风格的图表。

plt.style.use("ggplot")   
#风格使用ggplot

        但是值得注意的是,这里所说的绘制ggplot风格,是看起来像ggplot表格,实际上Matplotlib还是不具备图层风格,只是图片外观看起来像。

4.2 基本图标 Quick Start

(1) 类别型图表

        通常称X类别下Y数值之间的比较的图表为类别型图表,主要表现形式为X轴为类别型数据,Y轴为数值型数据。

        类别型图表常常有:柱状图、条形图(横向柱状图)、堆叠柱状图、极坐标的柱状图、词云图、雷达图、桑基图等等。

Matplotlib绘制不同房价对比的柱状图

"""
x 表示x坐标,数据类型为int或float类型,也可以为str
height 表示柱状图的高度,也就是y坐标值,数据类型为int或float类型
width 表示柱状图的宽度,取值在0~1之间,默认为0.8
bottom 柱状图的起始位置,也就是y轴的起始坐标
align 柱状图的中心位置,"center","lege"边缘
color 柱状图颜色
edgecolor 边框颜色
linewidth 边框宽度
tick_label 下标标签
log 柱状图y周使用科学计算方法,bool类型
orientation 柱状图是竖直还是水平,竖直:"vertical",水平条:"horizontal"
 """
data = pd.DataFrame({
    'city': ['深圳', '上海', '北京', '广州', '成都'],
    'house_price(w)': [3.5, 4.0, 4.2, 2.1, 1.5]
})
#创建一个绘图对象 fig,并设置绘图区域的大小为 10x6 英寸。
fig = plt.figure(figsize=(10, 6))
#在绘图对象上添加一个子图 ax1,指定子图的位置和大小为 [left, bottom, width, height],即左下角坐标 (0.15, 0.15),宽度和高度分别为 0.7。
ax1 = fig.add_axes([0.15, 0.15, 0.7, 0.7])
#使用 ax1.bar() 函数创建柱状图。data['city'] 是 x 轴数据(城市),data['house_price(w)'] 是 y 轴数据(房价),柱子的宽度设置为 0.6,柱子对齐方式为居中,标签为‘城市’。
bars = ax1.bar(data['city'], data['house_price(w)'], width=0.6, align='center', label='城市')
ax1.set_title('不同城市的房价对比图')
ax1.set_xlabel('城市')
ax1.set_ylabel('房价/w')
#启用网格线,参数 visible=True 表示显示网格线,which='both' 表示同时显示主网格线和次网格线。
ax1.grid(visible=True, which='both')
#为每个柱子添加高度值标签。通过循环遍历每个柱子,获取其高度,并使用 ax1.text() 函数在柱子顶部添加文本标签。标签位置设置为柱子中心偏左一点,且高于柱子顶部一点。
for bar in bars:
    height = bar.get_height()
    ax1.text(bar.get_x() + bar.get_width() / 2 - 0.1, height + 0.01, f'{height}', ha='center', va='bottom', fontsize=13)
#添加图例,使用 ax1.legend() 函数显示图例。
ax1.legend()
plt.show()

Plotnine绘制不同房价对比的柱状图

from plotnine.themes import element_text
# 数据准备
data = pd.DataFrame({
    'city': ['深圳', '上海', '北京', '广州', '成都'],
    'house_price(w)': [3.5, 4.0, 4.2, 2.1, 1.5]
})

# 创建单系列柱状图
p_single_bar = (
    ggplot(data, aes(x='city', y='house_price(w)', fill='city', label='house_price(w)')) +
    geom_bar(stat='identity') +  # 使用geom_bar绘制柱状图,stat='identity'表示使用原始数据
    labs(x="城市", y="房价(w)", title="不同城市的房价对比图") +  # 设置x轴标签、y轴标签和图表标题
    geom_text(nudge_y=0.08) +  # 在柱子上方添加文本标签,nudge_y用来调整文本位置
    theme(text=element_text(family="SimHei"))  # 设置主题,指定文本字体为SimHei(黑体),以支持中文显示
)

# 打印图表
# print(p_single_bar)
p_single_bar.show()

 Matplotlib绘制多系列柱状图:不同城市在不同年份的房价对比

 ## Matplotlib绘制多系列柱状图:不同城市在不同年份的房价对比
data = pd.DataFrame({
 '城市':['深圳', '上海', '北京', '广州', '成都', '深圳', '上海', '北京', '广州', '成都'],
 '年份':[2021,2021,2021,2021,2021,2022,2022,2022,2022,2022],
 '房价(w)':[3.5, 4.0, 4.2, 2.1, 1.5, 4.0, 4.2, 4.3, 1.6, 1.9]
 })
fig = plt.figure(figsize=(10,6))
ax1 = fig.add_axes([0.15,0.15,0.7,0.7])  
# [left, bottom, width, height], 它表示添加到画布中的矩形区域的左下角坐标(x, y),以及宽度和高度
plt.bar(
np.arange(len(np.unique(data['城市'])))-0.15, 
data.loc[data['年份']==2021,'房价(w)'], 
width=0.3, 
align='center', 
orientation='vertical', 
label='年份:2021'
)
plt.bar(
np.arange(len(np.unique(data['城市'])))+0.15, 
data.loc[data['年份']==2022,'房价(w)'], 
width=0.3, 
align='center', 
orientation='vertical', 
label='年份:2022'
)
plt.title("不同城市的房价对比图")   # 在axes1设置标题
plt.xlabel("城市")    
# 在axes1中设置x标签
plt.ylabel("房价/w")    
# 在axes1中设置y标签
plt.xticks(np.arange(len(np.unique(data['城市']))), np.array(['深圳', '上海', '北京', '广州', '成都']))
plt.grid(visible=True, which='both')  # 在axes1中设置设置网格线
data_2021 = data.loc[data['年份']==2021,:]
for i in range(len(data_2021)):
     plt.text(i-0.15-0.05, data_2021.iloc[i,2]+0.05, 
data_2021.iloc[i,2],fontsize=13)   
# 添加数据注释
data_2022 = data.loc[data['年份']==2022,:]
for i in range(len(data_2022)):
     plt.text(i+0.15-0.05, data_2022.iloc[i,2]+0.05, data_2022.iloc[i,2],fontsize=13) 
     #iloc[i, 2] 表示使用整数位置访问 data_2022 DataFrame 的第 i 行和第 2 列的元素。
plt.legend()
plt.show()

        其中对于这串代码的解释:

for i in range(len(data_2022)):
     plt.text(i+0.15-0.05, data_2022.iloc[i,2]+0.05, data_2022.iloc[i,2],fontsize=13) 
     #iloc[i, 2] 表示使用整数位置访问 data_2022 DataFrame 的第 i 行和第 2 列的元素。(0,1,2)

        在 Matplotlib 中,plt.text 用于在指定的坐标位置添加文本注释。函数原型如下:

matplotlib.pyplot.text(x, y, s, **kwargs)
  • x:文本的 x 坐标。
  • y:文本的 y 坐标。
  • s:要显示的文本字符串。
  • **kwargs:其他可选参数,例如字体大小、颜色等。
  • 通过循环遍历 data_2022 DataFrame 的每一行。
  • i 是当前柱子的索引(位置)
  • 0.15-0.05 是对 x 坐标的调整,以便将文本准确地放在柱子的顶部中央。
  • data_2022.iloc[i, 2] 是当前柱子的高度(第三列的值)。
  • + 0.05 是对 y 坐标的调整,使文本略高于柱子的顶部,避免文本与柱子重叠。
  • 文本内容:data_2022.iloc[i, 2]:显示当前柱子的高度值。

  • 其他参数:fontsize=13;设置文本字体大小为 13。

Plotnine绘制多系列柱状图:不同城市在不同年份的房价对比

## Plotnine绘制多系列柱状图:不同城市在不同年份的房价对比
data = pd.DataFrame({'城市':['深圳', '上海', '北京', '广州', '成都', '深圳', '上海', '北京', '广州', '成都'],
 '年份':[2021,2021,2021,2021,2021,2022,2022,2022,2022,2022],
 '房价(w)':[3.5, 4.0, 4.2, 2.1, 1.5, 4.0, 4.2, 4.3, 1.6, 1.9]
 })
data['年份'] = pd.Categorical(data['年份'], ordered=True, categories=data['年份'].unique())
p_mult_bar = (
 ggplot(data, aes(x='城市', y='房价(w)', fill='年份'))+
 geom_bar(stat='identity',width=0.6, position='dodge')+
 scale_fill_manual(values = ["#f6e8c3", "#5ab4ac"])+
 labs(x="城市", y="房价(w)", title="不同城市的房价对比图")+
 geom_text(aes(label='房价(w)'), position = position_dodge2(width = 0.6,preserve = 'single'))+
 theme(text = element_text(family = "SimHei"))
 )
p_mult_bar.show()

 

Matplotlib绘制堆叠柱状图:不同城市在不同年份的房价对比

# Matplotlib绘制堆叠柱状图:不同城市在不同年份的房价对比
data = pd.DataFrame({
 '城市':['深圳', '上海', '北京', '广州', '成都', '深圳', '上海', '北京', '广州', '成都'],
 '年份':[2021,2021,2021,2021,2021,2022,2022,2022,2022,2022],
 '房价(w)':[3.5, 4.0, 4.2, 2.1, 1.5, 4.0, 4.2, 4.3, 1.6, 1.9]
 })
#将城市和年份设置为索引,并使用房价(w)的值进行unstack 操作,使得年份成为列标签,房价(w) 为相应的值。
tmp=data.set_index(['城市','年份'])['房价(w)'].unstack()
data=tmp.rename_axis(columns=None).reset_index()
data.columns = ['城市','2021房价','2022房价']
# 创建数据框,包含城市、年份和房价数据
data = pd.DataFrame({
    '城市': ['深圳', '上海', '北京', '广州', '成都', '深圳', '上海', '北京', '广州', '成都'],
    '年份': [2021, 2021, 2021, 2021, 2021, 2022, 2022, 2022, 2022, 2022],
    '房价(w)': [3.5, 4.0, 4.2, 2.1, 1.5, 4.0, 4.2, 4.3, 1.6, 1.9]
})
# 将数据框的'城市'和'年份'列设置为索引,然后对'房价(w)'列进行数据透视
# `unstack()` 将'年份'列转化为列标签
tmp = data.set_index(['城市', '年份'])['房价(w)'].unstack()
# 重置索引,并将列名称重命名为'城市'、'2021房价'和'2022房价'
data = tmp.rename_axis(columns=None).reset_index()
data.columns = ['城市', '2021房价', '2022房价']
# 创建一个图形对象,设置图形的尺寸
plt.figure(figsize=(10, 6))
# 绘制柱状图,第一个系列(2021年的房价)
# 'data['2021房价']' 提供了2021年每个城市的房价数据
plt.bar(
    data['城市'],  
    data['2021房价'], 
    width=0.6,  # 设置柱的宽度
    align='center',  # 设置对齐方式
    orientation='vertical',  # 设置柱状图方向为垂直
    label='年份:2021'  # 设置图例标签
)
# 绘制柱状图,第二个系列(2022年的房价)
# `bottom` 参数使得2022年的柱子堆叠在2021年的柱子上面
plt.bar(
    data['城市'], 
    data['2022房价'], 
    width=0.6, 
    align='center', 
    orientation='vertical', 
    bottom=data['2021房价'],  # 堆叠在2021年房价柱子上
    label='年份:2022'  # 设置图例标签
)
# 设置图形的标题
plt.title("不同城市2021-2022年房价对比图")  # 这里将标题中的“2121”更正为“2021”
# 设置x轴和y轴的标签
plt.xlabel("城市")    
plt.ylabel("房价/w")    
# 添加图例以标识不同年份的柱子
plt.legend()
# 显示图形
plt.show()

# 创建数据框,包含城市、年份和房价数据
data = pd.DataFrame({
    '城市': ['深圳', '上海', '北京', '广州', '成都', '深圳', '上海', '北京', '广州', '成都'],
    '年份': [2021, 2021, 2021, 2021, 2021, 2022, 2022, 2022, 2022, 2022],
    '房价(w)': [3.5, 4.0, 4.2, 2.1, 1.5, 4.0, 4.2, 4.3, 1.6, 1.9]
})
# 将数据框的'城市'和'年份'列设置为索引,然后对'房价(w)'列进行数据透视
# `unstack()` 将'年份'列转化为列标签
tmp = data.set_index(['城市', '年份'])['房价(w)'].unstack()
# 重置索引,并将列名称重命名为'城市'、'2021房价'和'2022房价'
data = tmp.rename_axis(columns=None).reset_index()
data.columns = ['城市', '2021房价', '2022房价']
# 创建一个图形对象,设置图形的尺寸
plt.figure(figsize=(10, 6))
# 绘制柱状图,第一个系列(2021年的房价)
# 'data['2021房价']/(data['2021房价']+data['2022房价'])' 计算了每个城市2021年房价占总房价的比例
plt.bar(
    data['城市'],  
    data['2021房价'] / (data['2021房价'] + data['2022房价']), 
    width=0.4,  # 设置柱的宽度
    align='center',  # 设置对齐方式
    orientation='vertical',  # 设置柱状图方向为垂直
    label='年份:2021'  # 设置图例标签
)
# 绘制柱状图,第二个系列(2022年的房价)
# `bottom` 参数用来堆叠第二个系列在第一个系列上面
plt.bar(
    data['城市'], 
    data['2022房价'] / (data['2021房价'] + data['2022房价']), 
    width=0.4, 
    align='center', 
    orientation='vertical', 
    bottom=data['2021房价'] / (data['2021房价'] + data['2022房价']),
    label='年份:2022'  # 设置图例标签
)
# 设置图形的标题
plt.title("不同城市2021-2022年房价对比图")
# 设置x轴和y轴的标签
plt.xlabel("城市")
plt.ylabel("房价占比")
# 添加图例以标识不同的柱子
plt.legend()
# 显示图形
plt.show()

 使用Mtaplotlib绘制雷达图:英雄联盟几位英雄能力对比

#使用Mtaplotlib绘制雷达图:英雄联盟几位英雄能力对比
data = pd.DataFrame({
    '属性':['血量','攻击力','攻速','物抗','魔抗'],
    '艾希':[3,7,8,2,2],
    '诺手':[8,6,3,6,6]
})
plt.figure(figsize=(8,8))
theta = np.linspace(0,2*np.pi,len(data),endpoint=False)    #每个坐标点的位置
theta =np.append(theta,theta[0])    #让数据封闭
aixi = np.append(data['艾希'].values,data['艾希'][0])    #让数据封闭
nuoshou = np.append(data['诺手'].values,data['诺手'][0])    #让数据封闭
shuxing = np.append(data['属性'].values,data['属性'][0])    #让数据封闭
plt.polar(theta,aixi,'ro-',lw=2,label='艾希')    #画出雷达图的点和线
plt.fill(theta,aixi,facecolor='r',alpha=0.5)    #填充
plt.polar(theta,nuoshou,'bo-',lw=2,label='诺手')     #画出雷达图的点和线
plt.fill(theta,nuoshou,facecolor='b',alpha=0.5)     #填充
plt.thetagrids(theta/(2*np.pi)*360,shuxing)    #为每个轴添加标签
plt.ylim(0,10)
plt.legend()
plt.show()

(2)关系型图表

        关系型图标一般表现为X数值与Y数值之间的关系。例如:是否是线性关系,是否有正向相关关系等等。 

        一般来说关系可以分为:数值型关系、层次型关系和网络型关系。

使用Matplotlib和四个图说明相关关系

#使用Matplotlib和四个图说明相关关系
x=np.random.randn(100)*10
y1=np.random.randn(100)*10
y2=2*x+1+np.random.randn(100)
y3=-2*x+1+np.random.randn(100)
y4=x**2+1+np.random.randn(100)
plt.figure(figsize=(12,12))

plt.subplot(2,2,1)    #创建两行两列的子图,并绘制第一个子图
plt.scatter(x,y1,c='dodgerblue',marker='.',s=50)
plt.xlabel('x')
plt.ylabel('y1')
plt.title('y1与x不存在关联关系')

plt.subplot(2,2,2)    #创建两行两列的子图,并绘制第二个子图
plt.scatter(x,y2,c='tomato',marker='o',s=10)
plt.xlabel('x')
plt.ylabel('y2')
plt.title('y2与x存在关联关系')

plt.subplot(2,2,3)    #创建两行两列的子图,并绘制第三个子图
plt.scatter(x,y3,c='magenta',marker='o',s=10)
plt.xlabel('x')
plt.ylabel('y3')
plt.title('y3与x存在关联关系')

plt.subplot(2,2,4)    #创建两行两列的子图,并绘制第四个子图
plt.scatter(x,y4,c='deeppink',marker='s',s=10)
plt.xlabel('x')
plt.ylabel('y4')
plt.title('y4与x存在关联关系')

plt.show()

使用Plotnine和四个图说明相关关系 

#使用Plotnine和四个图说明相关关系
x=np.random.randn(100)*10
y1=np.random.randn(100)*10
y2=2*x+1+np.random.randn(100)
y3=-2*x+1+np.random.randn(100)
y4=x**2+1+np.random.randn(100)

df = pd.DataFrame({
    'x':np.concatenate([x,x,x,x]),
    'y':np.concatenate([y1,y2,y3,y4]),
    'class':['y1']*100+['y2']*100+['y3']*100+['y4']*100
})
df

 

p1=(
    ggplot(df)+
    geom_point(aes(x='x',y='y',fill='class',shape='class'),color='black',size=2)+
    scale_fill_manual(values=('#00AFBB','#FC4E07','#00589F','#F68F00'))+
    theme(text = element_text(family='SimHei'))
)
p1.show()

使用Matplotlib绘制具备趋势线的散点图

# 生成数据
x = np.linspace(-10, 10, 100)
y = np.square(x) + np.random.randn(100) * 100
# 线性回归
x_reshape = x.reshape(-1, 1)
linear_regressor = LinearRegression()
y_linear_pred = linear_regressor.fit(x_reshape, y).predict(x_reshape)
# 二次多项式回归
poly_features = PolynomialFeatures(degree=2)
x_poly2 = poly_features.fit_transform(x_reshape)
poly_regressor = LinearRegression()
y_poly_pred = poly_regressor.fit(x_poly2, y).predict(x_poly2)
# 指数回归
y_exp_pred = LinearRegression().fit(np.exp(x_reshape), y).predict(np.exp(x_reshape))
# LOESS回归
y_loess_pred = sm.nonparametric.lowess(y, x, frac=2/3)[:, 1]
# 绘图
plt.figure(figsize=(8, 8))
# 带线性趋势线的散点图
plt.subplot(2, 2, 1)
plt.scatter(x, y, c='tomato', marker='o', s=10)
plt.plot(x, y_linear_pred, c='dodgerblue')
plt.xlabel('x')
plt.ylabel('y')
plt.title('带线性趋势线的散点图')
# 带二次多项式趋势线的散点图
plt.subplot(2, 2, 2)
plt.scatter(x, y, c='tomato', marker='o', s=10)
plt.plot(x, y_poly_pred, c='dodgerblue')
plt.xlabel('x')
plt.ylabel('y')
plt.title('带二次多项式趋势线的散点图')
# 带指数趋势线的散点图
plt.subplot(2, 2, 3)
plt.scatter(x, y, c='tomato', marker='o', s=10)
plt.plot(x, y_exp_pred, c='dodgerblue')
plt.xlabel('x')
plt.ylabel('y')
plt.title('带指数趋势线的散点图')
# 带LOESS趋势线的散点图
plt.subplot(2, 2, 4)
plt.scatter(x, y, c='tomato', marker='o', s=10)
plt.plot(x, y_loess_pred, c='dodgerblue')
plt.xlabel('x')
plt.ylabel('y')
plt.title('带LOESS趋势线的散点图')
plt.tight_layout()
plt.show()

使用Matplotlib绘制聚类散点效果图

#使用Matplotlib绘制聚类散点图
from sklearn.datasets import load_iris     #加载鸢尾花数据集
iris = load_iris()
X=iris.data
label = iris.target
feature = iris.feature_names
df = pd.DataFrame(X,columns=feature)
df['label']=label
df

label_unique = np.unique(df['label']).tolist()
#这行代码获取标签列中所有唯一的标签,并将它们转换为一个列表。 np.unique(df['label']) 返回一个包含唯一标签的NumPy数组, tolist() 方法将其转换为Python列表。

plt.figure(figsize=(10,6))
for i in label_unique:
    df_label =df.loc[df['label']==i,:]
    plt.scatter(x=df_label['sepal length (cm)'],y=df_label['sepal width (cm)'],s=20,label=i)
plt.xlabel('sepal length(cm)')
plt.ylabel('sepal width(cm)')
plt.title('sepal width(cm)~sepal length(cm)')
plt.legend()
plt.show()

使用Plotnine绘制相关系数矩阵图

        首先需要了解的是,corr()函数用于计算DataFrame中不同列之间的相关性系数。其基本语法如下:

DataFrame.corr(method='pearson', min_periods=1, numeric_only=True)

1.method:用于指定计算相关性系数的方法,默认为’pearson’(皮尔逊相关系数)。其他可选方法包括’kendall’(肯德尔等级相关系数)和’spearman’(斯皮尔曼等级相关系数)。
2.min_periods:用于指定在计算相关性之前,每对列需要至少有多少个非空值。默认为1。
3.numeric_only:指定是否仅对数值型列进行计算。默认为True。

import numpy as np
import pandas as pd
from plotnine import *
from plotnine.data import mtcars

# 去除第一列
mtcars_numeric = mtcars.iloc[:, 1:]
# 计算相关系数矩阵并重置索引
corr_mat = np.round(mtcars_numeric.corr(), 1).reset_index()
# 将宽数据转换为长数据
df = pd.melt(corr_mat, id_vars='index', var_name='variable', value_name='corr_xy')
# 计算相关系数的绝对值
df['abs_corr'] = np.abs(df['corr_xy'])

p1=(
    ggplot(df,aes(x='index',y='variable',fill='corr_xy',size='abs_corr'))+
    geom_point(shape='o',color='black')+
    scale_size_area(max_size=11,guide=None)+ # 这里设置 guide=None 来去掉 size 图例
    scale_fill_cmap(name='RdYIBu_r')+
    coord_equal()+
    labs(x='Variable',y='Variable')+
    theme(dpi=100,figure_size=(4.5,4.55))
)

p2=(
    ggplot(df,aes(x='index',y='variable',fill='corr_xy',size='abs_corr'))+
    geom_point(shape='s',color='black')+
    scale_size_area(max_size=10,guide=None)+ # 这里设置 guide=None 来去掉 size 图例
    scale_fill_cmap(name='RdYIBu_r')+
    coord_equal()+
    labs(x='Variable',y='Variable')+
    theme(dpi=100,figure_size=(4.5,4.55))
)

p3=(
    ggplot(df,aes(x='index',y='variable',fill='corr_xy',label='corr_xy'))+
    geom_tile(color='black')+
    #scale_size_area(max_size=11,guide=None)+ # 这里设置 guide=None 来去掉 size 图例
    geom_text(size=8,color='white')+
    scale_fill_cmap(name='RdYIBu_r')+
    coord_equal()+
    labs(x='Variable',y='Variable')+
    theme(dpi=100,figure_size=(4.5,4.55))
)
print([p1,p2,p3])

p1 

p2 

 

p3

 

使用Matplotlib/Seaborn绘制相关系数矩阵图

#使用Matplotlib/Seaborn绘制相关系数矩阵图
uniform_data = np.random.rand(10,12)
sns.heatmap(uniform_data)

(3) 分布型图表

        所谓数据的分布,就是看数据的密集稀疏程度;描述数据密集或者稀疏的情况实际上可以用频率或者概率。

  • 9
    点赞
  • 13
    收藏
    觉得还不错? 一键收藏
  • 2
    评论
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值