🔛 引八月初的时候，博导跟我介绍了这个项目，大概内容就是对某市供水管网中的泵站与清水池进行优化调度。第一次博导跟我探讨这个项目的时候，我还是有点忐忑，觉得万一自己干不好岂不是砸了老师的招牌。但是我导说了句特别有安全感的话：没事，只管干吧，搞不出来有我给你兜着！安全感直接拉满，奥里给，干！还得是我导儿。历时4个多月，十几次的每周PPT工作汇报，终于完成了模型的算法设计、实际运维逻辑匹配、数据库输入输出对接、线上部署等工作。在这个项目过程中，也经历过很多次的绝望、失落、喜悦，以及无数次的DEBUG，感觉自己还是学到了很多东西的。这篇博客就来记录下整个项目下来我自己的一些学习与心得体会。在上次的组会过程中，我其实已经进行了一个大概的总结，所以说接下来就按照组会PPT中的内容进行一个记录吧~组会汇报⏱️ 考虑时间成本Quote / 参考面对未知的事物，试错过程需要花费大量时间成本，凡事一定要提前准备，预备好抵抗未知错误的缓冲区。对于任何一件未知的事（这个未知是指没有别人做过类似的事可以给你参考），那么就一定要提前考虑到遇到未知错误所花费的时间成本。在项目开始前，博导说预备时间是两个半月，然后我们就先按两个月的来，给最后留半个月的时间缓冲。我当时觉得时间还是很充裕的。然而实际情况做下来，我们花费了将近4个多月，几乎是原来时间的两倍了。首先是第一个问题：项目开始前一定要做充分的调研，大方向不能错，或者说尽量少错。在项目开始前，与甲方多次沟通交流，确定了的方案是“基于历史数据驱动的智能调度模型研究”。从八月到国庆节前，我们也是一直按照这个方向来做的。结果到了国庆节前发现，我们所走的路是存在致命问题的，不得已又改变路线，选择了“基于区间流量平衡与遗传算法的智能调度模型研究”，这也就意味着我们前期浪费了巨多时间做了无用功。其实选择的这个新方案，我导在第一次跟我讲这个项目的时候就提到了这种方案。因为对于一个泵站来说，其水泵开关的组合非常适合使用遗传算法（Genetic Algorithm, GA）进行求解。只是我当时根本不知道遗传算法是个什么东西，所以没有GET到我导的意思。现在看来，我导对大方向的把握真是有点东西，如果一开始就按照我导的思路做，也许我们能减少很多无用功了。第二个感受比较深的过程是在做项目前后端数据库对接的时候。简单来说，我负责模型的构建，模型输入输出的数据需要从数据库中读取或写入，另外一个团队负责WEB前端的开发，前端显示的内容就是从数据库中读取的。这里，犯了一个很严重的逻辑问题。直接与数据库对接的是前端和模型，那么前端就应该在与模型商议数据库结构之前，先确定了用户的需求。但是实际上前端团队并不是这样做的，而是边与用户沟通加入新的需求边与模型端沟通确定写入数据格式。这就造成了一个问题，前端一边与模型一起对接数据库，一边与顾客商议需求，这就给我的模型带来了很大困扰，我对我的模型输出进行了多次很大结构上的调整，花费了很多时间。所以说，选择与靠谱的团队进行合作还是很重要的😑。🤌 嘴上说，心里想, 不如手上做Quote / 参考不知道什么行不行，那就先动手试试，不动手永远不知道结果。动手了还真可能一个个慢慢就解决了。上面也提到了，在国庆前后我们准备更换新的工作路线。我导跟我提到了GA，说GA非常适合我们当时所遇到的问题，但是当时的我对于GA一窍不通。对于未知的事物，人们的本能总是恐惧，所以总是想逃避。当时我的一听到“算法”两个字，似乎就感觉一座大山压在我的身上，要让我迅速学习一个算法并应用到实际项目，这对我来说太难了。所以，那天我还是想看看以前的路到底还有没有办法能够继续走得通。一直到最后，发现老路实在是走不下去了，我才“逼上梁山”开始学习GA。但是当我真的开始学习GA，发现它其实并没有想象中的那么难。从初识到应用到自己的项目里，前后也就是三五天的时间。通过这件事就说明两点：要对自己有信心，面对未知事物，肯定会遇到困难。但是要相信自己——只用功夫深，铁棒磨成针！！！凡事要一点点来，不用心急，每天进步一点点——水滴石穿！！！🐛 先谋大略，后顾细微Quote / 参考技术应用于实际工程，需要考虑太多太多太多的未知情况，前期尽量去考虑照顾到，大方向不要出错误，后期预留一定的时间慢慢调试实际运行中的小问题。如果一项技术的科学研究已经发展到了98%，那么能将其30%成功应用到实际就很不错了， 科学与技术之间还是有很多隔阂的。这个隔阂，其实指的就是将科学应用到实际过程中需要协调的未知情况，因为科学实验大部分是在较为理想的环境中进行的。可能并没有什么好的办法能够消除这种隔阂的，只能说是在技术应用的过程中就尽量多的去考虑实际中可能遇到的情况。然后，还要在实际部署后，预留大量的时间去DEBUG。💡 君子生非异也，善假于物也Quote / 参考擅用ChatGPT！！！特别是数据库相关内容（猜测可以推广到Numpy、Pandas的数据分析），把问题描述清楚（可以采用分步提问的方式）。从我小学接触互联网以来，互联网第一次给我的震撼是CS Online（一款在线射击游戏），第二次就是今年的ChatGPT。这个项目的全过程中，涉及到数据库的操作全部是ChatGPT帮助我完成的，对于数据库这种本身其结构性就很强的东西，ChatGPT的回答相当完美，几乎没有出现过任何错误。ChatGPT另一个比较好的应用场景是基础知识的学习，最近在学习图论的过程中，我没有看过任何视频教程，全部都是跟着ChatGPT进行学习的。这让我的学习效率得到了极大得提高。前几天毫不犹豫下单了ChatGPT 4.0，这个价格能够提供给我的帮助让我觉得非常实惠。⌨️ Python程序设计指南（非正式版）Quote / 参考项目程序设计一定要最起码有一套“自己的规范”。结构化设计这个项目核心文件将近3000行代码，算是我目前写过最大的项目了。在项目开发的过程中，越发觉得程序设计规范化的重要性。当你写10行代码时，没有什么需要注意的；但是当你写了1000行代码时，回过头来可能你自己都忘了刚开始写的是什么。这个时候就凸显了注释和结构的重要性。这里所说的结构化其实并不一定要按照Python编程规范或者说某一种固定的格式来编写代码。我个人觉得是最起码你自己要有一套自己的规范（如果你是单打独斗的话）。例如，我就喜欢对每个类设置一个run()函数来作为类对象的入口函数：if __name__ == '__main__': optimizer = Optimizer( log_level=logging.INFO # logging.INFO, logging.DEBUG ) optimizer.run()虚拟环境在大项目开发过程中，虚拟环境真的太重要了，主要体现在以下三点：有可能项目需要部署到其他电脑来测试，如果没有虚拟环境，那么就需要在新电脑重新配置一个与原电脑一样的环境，非常麻烦；避免了其他项目环境对本项目的污染，也避免了本项目对系统环境的污染；如果程序需要打包为EXE或者其他可执行文件，那么虚拟环境可以避免无用的包被打包进程序里，从而减小程序的体积。如何导出当前环境中的所有包及其版本号？使用命令： pip freeze > .\requirements.txt即可。会在根目录生成一个名为requirements.txt的文件。将该文件复制到需要配置新环境的电脑中，使用命令：pip install -r requirements.txt即可一键复制安装源环境中的所有包。EXE程序的打包PyInstaller 是最常用和最方便的打包工具，尤其是由于其对不同操作系统的广泛支持和对各种Python版本和库的兼容性。使用PyInstaller打包Python程序相对简单，社区支持也很强大。PyInstaller可以打包Python程序为独立的可执行文件，适用于Windows、macOS和Linux。它支持Python 2.7和Python 3.5及以上版本，能够打包大多数Python库，并且有很多选项来自定义打包过程，下面记录几个打包过程踩的坑。打包程序(venv) PS D:\Code\Optimizer> pyinstaller -w -i favicon.ico --add-data="data;data" main.py参数解释如下：-w：让程序以无命令行的形式运行；-i favicon.ico：将根文件夹中的favicon.ico作为程序的LOGO；--add-data="data;data"将data文件夹打包为程序下的data文件夹。异常捕获Python程序在经过打包后，如果有异常就会直接闪退，非常不利于调试错误。这个时候我们可以使用：try: # 尝试捕获异常 ... except KeyError as e: # 成功捕获异常 print(e) else: # 未发现异常 ... finally: # 最终要执行的 ...来捕获并处理异常。路径处理Python程序如果需要涉及读写文件的操作，那么就必须使用下面的函数来兼容EXE运行时的路径问题：@staticmethod def resource_path(*args): """ 获取文件路径, 兼容EXE打包问题 :param args: 任意多个参数, 文件路径. 例如访问 './data/Setting_CQ.json'文件, 需要输入两个参数'data', 'Setting_CQ.json' :return: 组合后的绝对路径 """ if hasattr(sys, '_MEIPASS'): # 在PyInstaller打包后的可执行文件中 return os.path.join(sys._MEIPASS, *args) else: # 在原始Python脚本中 return os.path.join(os.path.abspath("."), *args)💨 结总体上来说，这个过程还是挺锻炼人的，我在这个过程中学习到了很多东西，也成长了很多，也算是自己的一笔财富吧。加油！（不过我还是希望后面甲方反馈的BUG能少一点~🙏）🙏🙏🙏佛祖保佑🙏🙏🙏print(" _ooOoo_ ") print(" o8888888o ") print(" 88 . 88 ") print(" (| -_- |) ") print(" O\\ = /O ") print(" ____/`---'\\____ ") print(" . ' \\| |// `. ") print(" / \\||| : |||// \\ ") print(" / _||||| -:- |||||- \\ ") print(" | | \\\\\\ - /// | | ") print(" | \\_| ''\\---/'' | | ") print(" \\ .-\\__ `-` ___/-. / ") print(" ___`. .' /--.--\\ `. . __ ") print(" ."" '< `.___\\_<|>_/___.' >'"". ") print(" | | : `- \\`.;`\\ _ /`;.`/ - ` : | | ") print(" \\ \\ `-. \\_ __\\ /__ _/ .-` / / ") print(" =======`-.____`-.___\\_____/___.-`____.-'======= ") print(" `=------=' ") print(" ") print(" ............................................. ") print(" 佛祖镇楼 BUG消失 ") print(" 佛曰: ") print(" 写字楼里写字间，写字间里程序员； ") print(" 程序人员写程序，又拿程序换酒钱。 ") print(" 酒醒只在网上坐，酒醉还来网下眠； ") print(" 酒醉酒醒日复日，网上网下年复年。 ") print(" 但愿老死电脑间，不愿鞠躬老板前； ") print(" 奔驰宝马贵者趣，公交自行程序员。 ") print(" 别人笑我忒疯癫，我笑自己命太贱； ") print(" 不见满街漂亮妹，哪个归得程序员？")
Introduction / 引言大学期间用来打发无聊时间学的Python没想到竟然在写毕业论文的时候用处这么大，整个硕士论文所做研究，从前期的数据整理、数据分析，到最后的数据可视化我基本上都使用Python来完成，这篇博客就来分享下我毕业论文课题中所做数据分析相关的Python代码。本博文所有相关的代码都上传在GitHub仓库：Data-Analysis-for-Thesis，如果帮到你了，记得给我来个Star☺️，也可以顺便去参观下我的个人博客亚灿网志。Python环境配置环境安装首先是在Python官网下载你计算机对应的Python软件，然后安装。安装过程基本都是傻瓜式，不做过多叙述，一路回车即可。Python官网IDE选择然后就是IDE（Integrated Development Environment）的配置，最常见的肯定是VSCode，特点是足够轻量化且免费。但是我选择的是PyCharm，我选择它的理由是以前学Python的时候就用的它，因此快捷键什么的都比较熟悉。PyCharm正版是收费的，教育版可以申请一年的使用权限，破解版的自行百度即可。VSCodePyCharmWarning / 注意Pycharm也是可以申请教育版嘀，但是教育版并不支持编写Jupyter Notebook的.ipynb文件。Tips / 提示关于IDE的选择大家没必要纠结，新手不太了解想简单上手的话就选择VSCode；用过PyCharm并且熟练怎么破解安装的那就选择PyCharm。包的安装数据矩阵分析及处理：Pandas、Numpy、Math、Scipy；绘图可视化：Matplotlib、Seaborn；其他包：hues可以在控制台打印出彩色的提示信息，用法也比较简单，官方手册。hues的简单用法数据整理及预处理数据介绍论文所做内容简单来说就是对不同样本点的多个指标进行相关性分析，具体的数据格式就是：样本点变量1变量2变量nD-H1122424...L-N2324232...............基本常量的定义IMG_TYPE = '.svg' # 出图格式，我选择矢量图svg DATA_PATH, EXPORT_PATH, DPI = './data/', './export/', 300 # 数据存储路径、出图写入路径、出图DPI ROW_NUM, COL_NUM = 5, 3 # 对于多子图的图片，定义默认布图网格为5×3文件读取读取csv文件需要使用pandas的pd.read_csv()方法，具体的参数有：index_col：设置行索引为哪一列，可以使用序号或者列名称；sep：csv文件中的分隔符，默认常见的用法都可以自动识别，不需要设置；header：设置表头，参数为None就是没有表头，设置为n就是把第n行读取为表头；names：设置列名称，参数为list；usecols：仅读取文件内某几列。Quote / 参考具体用法可以参考李庆辉所著《深入浅出Pandas——利用Python进行数据处理与分析》3.2章 读取CSV（PDF P89）。数据表合并首先遇到的第一个需求就是，所有样本点的列变量存储在不同的数据表中，比如，样本点的指标分为上覆水的指标与沉积物的指标两部分，分别存储在两个或者多个数据表中，那么如何将两个或者多个数据表进行合并呢？all_df_0 = pd.concat([gene_abundance, sediment_env_df, eems_df, overlying_water_env_df, sediment_size_df], axis=1) # 将五个数据表按照行索引合并Quote / 参考具体用法可以参考李庆辉所著《深入浅出Pandas——利用Python进行数据处理与分析》7.2章 数据连接pd.concat（PDF P274）。根据行索引为每个样本点设置分类行索引也就是每个样本点的标记名，分别为“D-H1”、“L-N3”之类的，其中第一个字符的值为“D”、“L”、“W”，分别代表枯水期（Dry Season）、平水期（Level Season）、丰水期（Wet Season）；第二个字符为一个分隔符“-”；第三个字符的值为“N”、“P”、“H”，分别代表三条不同的河流，南淝河、派河、杭埠河；最后的数字代表是第几个样本点。例如：“D-H1”代表枯水期杭埠河第一个样本点、“L-N3”代表平水期南淝河第三个样本点。那么问题来了，我想要为合并后的数据表新增两列“River”、“Period”，分别来反应这个样本点的属性，应该如何实现呢？样本点变量1变量2变量nRiverPeriodD-H1122424...Hangbu RiverDry SeasonL-N2324232...Nanfei RiverLevel Season..................思路其实也很简单，就是使用apply函数分别对每一行（也就是每一个样本点）进行处理，获取该行的行索引，然后对行索引的字符进行判断即可：all_df['Period'] = all_df.apply(lambda x: 'Dry' if 'D' in x.name else ('Wet' if 'W' in x.name else 'Level'), axis=1) + ' Season' all_df['River'] = all_df.apply(lambda x: 'Nanfei' if 'N' in x.name else ('Pai' if 'P' in x.name else 'Hangbu'), axis=1) + ' River'Tips / 提示这里使用了Python列表推导式相关的知识，具体的讲解可以看之前的博文：「Python」列表推导式。在对每一行的样本点添加River、Period变量后，会有一个问题，River、Period两列的数据都是Object字符串类型。这种数据类型有两个问题：如果数据矩阵有几十万行，那么这两列会占用很大的内存空间；对数据进行绘图过程中，我想把River变量按照Nanfei River、Pai River、Hangbu River的顺序排列，那么就很麻烦。因为字符串变量默认是按照首字母的顺序来进行排序的，默认排序是Hangbu River、Nanfei River、Pai River。为了解决这两个问题，我们可以将这两列的数据由原来的object类型转换为Category类型，Category的好处就是，当数据量较大时，可以显著减小数据所占用的内存；第二还可以对数据类型进行排序。具体的处理方法：# 对季节、河流两列进行排序，首先定义category类型顺序 river_order = CategoricalDtype( # 河流的顺序定义为南淝河、派河、杭埠河 ['Nanfei River', 'Pai River', 'Hangbu River'], ordered=True ) period_order = CategoricalDtype( # 时期的顺序定义为枯水期、平水期、丰水期 ['Dry Season', 'Level Season', 'Wet Season'], ordered=True ) # 将两列Object类型数据转换为category类型并排序 all_df['River'] = all_df['River'].astype(river_order) all_df['Period'] = all_df['Period'].astype(period_order) all_df = all_df.sort_values(by=['Period', 'River', 'SortNum'])Warning / 注意在将一列数据转化为Category对象后，如果数据表中没有某个Category，但是绘图的时候还是会占用一个位置，下面举例说明：这个数据表中的Period列已经不包含Level Season的数据，但是使用.value_counts()方法统计的时候还是会对其进行统计。如果这个时候进行绘图，会出现一个问题：中间数据明明是空的，但是还是会占位。解决这个问题的方法是使用：nanfei_dry_and_wet_gene['Period'].cat.remove_unused_categories()删除Period列中没有用到的Category，这个时候再次进行绘图就可以正常显示了。异常值处理缺失值的填充Pandas中缺失值的填充所用方法时pd.fillna()，具体的参数可以填写：In [16]: pd.DataFrame.fillna Out[16]: <function pandas.core.frame.DataFrame.fillna(self, value: 'object | ArrayLike | None' = None, method: 'FillnaOptions | None' = None, axis: 'Axis | None' = None, inplace: 'bool' = False, limit=None, downcast=None) -> 'DataFrame | None'>value：直接将缺失值填充为字符串或者数字；method：填充方式，method='ffill' 向前填充，method='bfill'向后填充，也就是说用前面的值来填充NA或用后面的值来填充NA。另外，在使用读取pd.read_csv()读取csv文件的时候，也可以通过参数：na_values=Nonekeep_default_na=Truena_filter=True的设置来对NA值进行过滤或者识别。删除缺失值使用pd.DataFrame.dropna()方法完成缺失值的删除：In [17]: pd.DataFrame.dropna Out[17]: <function pandas.core.frame.DataFrame.dropna(self, axis: 'Axis' = 0, how: 'str' = 'any', thresh=None, subset: 'IndexLabel' = None, inplace: 'bool' = False)>通过参数how的属性值来设置：any：当每一行有一个缺失值时就删除这一行；all：当一行所有的数据都时缺失值时再删除这一行。重复值的删除使用pd.DataFrame.drop_duplicates()方法完成缺失值的删除：In [18]: pd.DataFrame.drop_duplicates Out[18]: <function pandas.core.frame.DataFrame.drop_duplicates(self, subset: 'Hashable | Sequence[Hashable] | None' = None, keep: "Literal['first'] | Literal['last'] | Literal[False]" = 'first', inplace: 'bool' = False, ignore_index: 'bool' = False) -> 'DataFrame | None'>通过参数keep的属性值来设置：first：所有重复行删除，保留第一行；last：所有重复行删除，保留最后一行。数据处理与可视化绘图前的小准备画图格式的定义如何在Matplotlib中显示中文：plt.rcParams['font.sans-serif'] = ['SimHei'] # 用于显示中文 plt.rcParams['axes.unicode_minus'] = False如何定义Seaborn绘图的默认字体属性：sns.set( style='darkgrid', # 绘图风格 font='Times New Roman', # 默认字体 font_scale=3 # 默认字体比例 )如何实现子图编号CHAR = [chr(i) for i in range(97, 123)] # 获取26个英文字母，用于给子图编号定义一个26个英文字母的list，循环绘制子图的时候直接调用即可。重复代码的打包每次进行数据分析我都会新建一个.ipynb文件，而数据分析前都需要经过数据表合并、数据清洗等工作，那么最好的方式其实是将数据分析前的准备工作进行一个打包，然后在.ipynb文件的第一行引入包即可。例如：我新建一个ResearchMain.py文件，然后将所有数据表合并、数据清洗的代码都放在这个文件里：# 引入数据分析常用的包 ... # 读取文件 ... # 合并文件 ... # 数据表清洗 ...然后在每次新建.ipynb文件进行数据分析时，我都会在第一行使用：from ResearchMain import *来引入所有ResearchMain.py文件中定义的变量与方法。QPCR标准曲线可视化Tips / 提示什么是标准曲线？简单来说，自变量x与因变量y之间存在某种线性关系——$y=ax+b$，那么我们可以通过多次改变自变量x的值，然后观察y的值并记录，得到几组对应的$x_1$、$x_2$、$x_3$、$x_4$、$x_5$、...与$y_1$、$y_2$、$y_3$、$y_4$、$y_5$、...，那么我们就可以通过得到的这几组数据来对自变量x与因变量y进行线性拟合，从而得到一个标准曲线$y=ax+b$，有了标准曲线之后，我们就可以直接输入任意的自变量x值，计算出因变量y的值。在Numpy中，拟合函数直接有现成的，可以直接调用：In [37]: x, y = [1, 2, 3, 4, 5], [2, 4, 6, 8, 11] # 需要进行拟合的自变量与因变量 In [38]: np.polyfit(x, y, 1) # 对自变量x与因变量y进行拟合，且拟合为一次函数 Out[38]: array([ 2.2, -0.4]) # 拟合结果为y=2.2x-0.4拟合完毕如何拼接拟合方程？In [25]: fitting_arr = np.polyfit(x, y, 1) In [26]: fitting_equ = np.poly1d(fitting_arr) # 获取拟合方程 In [27]: fitting_equ([6, 7, 8, 9, 10]) # 利用拟合方程计算任意自变量对应的因变量 Out[27]: array([12.8, 15. , 17.2, 19.4, 21.6]) In [28]: fitting_equ(11) Out[28]: 23.799999999999997 # 23.799999999999997=11*2.2-0.4获取拟合的$R^2$：In [43]: np.corrcoef(y, fitting_equ(x))[0, 1] ** 2 Out[43]: 0.9918032786885246$R^2$的计算原理其实就是把自变量带入拟合方程，将计算出来的因变量与原始的因变量进行比较，计算其相关性。搞懂上面的原理之后，进行标准曲线的可视化其实就很简单了。图中可以看出，还生成了一个拼接的一元一次方程，方程的拼接可以直接用我写好的函数，函数的具体用法以及讲解已经在注释里说的很清楚了：Tips / 提示函数的主要作用就是传入np.polyfit(X, Y, DEG)返回的list对象，比如返回的是[2, 3, -4]，那么就输出2x^2+3x-4，可以自动识别n元函数。def make_fit_equ_str(paras): """ 拼接拟合方程 :param paras: 传入参数为np.polyfit(X, Y, DEG)返回的对象 :return: 返回一个拼接拟合方程的字符串(LaTex格式) """ fit_equ_str = '' for i in range(len(paras)): # 遍历常数项的过程中，主要需要进行两步“组装”：①如果该常数项不是第一个常数项，且该常数项大于0，需要转化为字符串并在前面添加一个“+”； if paras[i]: # 如果常数项不为0(为0直接跳过这一项) cur_item = str(round(paras[i], 4)) # 每个常数项都保留两位小数，并转化为字符串 # 如果不是第一个常数项，且该常数项大于0，需要转化为字符串前面添加一个“+” if paras[i] > 0 and fit_equ_str != '': # 注意两种特殊情况：[0, 0, -5, 7]、[0, -5, 7]如何排除 cur_item = '+' + cur_item # ②为每一个常数项添加x的幂指数 if i == (len(paras) - 1) and cur_item != '0': # 最后一个项只有一个常数，不需要添加x了（如果为0就不用+了） fit_equ_str += cur_item else: # 前面的每一项都需要在最后加上一个x^幂数 idempotent_num = len(paras) - 1 - i # 对应的幂指数 if idempotent_num == 1: # 拟合为一次多项式 fit_equ_str += cur_item + 'x' else: # 拟合为二次及二次以上的多项式 fit_equ_str += cur_item + 'x^' + str(idempotent_num) #方程拼接 # fit_equ_str = '$' + fit_equ_str + '$' # 将多项式转化为LaTex公式 return fit_equ_strQPCR数据处理上一步的操作已经制作出来的功能基因的标准曲线，接下来就是把仪器测得的CT值（自变量）转化为丰度值（因变量），就是带进拟合的标准曲线计算下就OK了。具体绘图代码直接看GitHub代码即可，没有什么难度。三维荧光数据可视化读取数据表使用日立F-7000荧光光谱仪对沉积物中溶解性有机质（Dissolved Organic Matter, DOM）结构特征和组成成分进行表征。仪器得到的数据是.txt格式，且有用的数据表是从Data Points这一行后面开始的。得到的txt文件所以说，我们要先读取.txt文件，循环读取每一行，直到读取到Data Points这一行，说明已经到数据表了。def get_skip_rows(path): """ 读取txt文件，并在文件中查找含有'Data Points'的行，数据矩阵就在这一行的下面 :param path: 文件路径 :return: 数据矩阵开始的行号 """ f = open(path) for index, line in enumerate(f.readlines()): if 'Data Points' in line: return index + 1这个函数的作用就是输入.txt文件的路径，然后会返回需要的数据表是在第几行开始的。然后使用pd.read_table()方法读取.txt，并通过设置skiprows的值，来跳过前面无用的数据。例如利用get_skip_rows()函数获取到.txt文件中数据表从第156行开始：df = pd.read_table(search_info['Path'], skiprows=156, index_col=0)这样的话就可以完美跳过.txt文件前面无用的数据，直接读取所需的数据表。消除瑞利散射瑞利散射消除前后对比瑞利散射的消除其实很简单，观察数据表就可以看出来，瑞利散射其实就是不该出现在某个区间内数据峰，我们只需要慢慢根据$E_x$与$E_m$的设置范围来进行消除就行了：for i in range(modify_df.shape[0]): # 遍历ex for j in range(modify_df.shape[1]): # 遍历em ex = 200 + i * 5 # ex范围为200~450，间隔为5 em = 280 + j * 5 # em范围为280~550，间隔为5 if ex + 455 > em > ex + 145: modify_df.iloc[i, j] = 0 # 不能使用0, 因为0是有意义的数据光谱图的分区直接使用：cur_ax.axvline(250, color='white', linestyle='--', linewidth=5) # 垂直线 cur_ax.axhline(330, color='white', linestyle='--', linewidth=5) # 水平线就可以在图上画线了。需要讲解的其实就这几个部分，别的直接看代码就行了。粒径数据可视化沉积物粒径百分比分布使用Malvern Mastersizer 2000型激光粒度仪进行分析。得到数据后需要手动整理为.csv格式。数据格式首先来讲解下数据格式，每一列代表一个样本，每一行代表对应粒径所占百分比。例如图中红方框所示就是代表D-N4样本点对应粒径为0.955 μm颗粒占比为0.03%。清除空行上图中可以看出，数据有很多空行，那么首先第一步就是清除掉这些空行：df.dropna(how='all', inplace=True) # 删除缺失值（行全为空） df.dropna(how='all', axis=1, inplace=True) # 删除缺失值（列全为空）绘图还是直接看代码吧。常规指标数据可视化# 子图摆放设置 ROW_NUM = 2 COL_NUM = int(len(df.columns) / ROW_NUM) # 绘图设置 sns.set( style='darkgrid', font='Times New Roman', font_scale=2 ) # 布局设置 plt.figure(dpi=DPI) fig, ax_arr = plt.subplots(ROW_NUM, COL_NUM, sharex='col', figsize=(20, 10)) fig.subplots_adjust(hspace=0.2, wspace=0.2) # 循环绘图 for index, s in enumerate(df.iteritems()): cur_df = pd.concat([s[1], affix_info], axis=1) # 计算行列数 row_num = index // COL_NUM col_num = index - row_num * COL_NUM cur_ax = ax_arr[row_num][col_num] cur_plot = sns.boxplot( # 绘图 x="Period", # X轴 y=s[0], # Y轴 hue="River", # 颜色分类 data=cur_df, # 数据表 dodge=True, palette="Set1", # 配色 ax=cur_ax # 绘图坐标轴 ) # 子图图名、x轴、y轴、图例 cur_ax.tick_params(axis='x', rotation=90) if len(df.columns) - (index + 1) >= COL_NUM: cur_plot.set_xlabel(None) cur_plot.set_title(f'({CHAR[index]}) ${s[0][2:]}$ ({unit[s[0]]})',fontproperties="Times New Roman", fontsize=20) cur_ax.get_legend().remove() # 移除子图图例 cur_plot.set_ylabel(None) # 子图的图例相同，获取最后一个子图的图例 lines, labels = fig.axes[-1].get_legend_handles_labels() fig.legend(lines, labels, ncol=3, loc='lower center', bbox_to_anchor=(0.5, -0.16)) # 保存图片 # fig.savefig(EXPORT_PATH + '沉积物环境因子时空变化' + IMG_TYPE, dpi=DPI, bbox_inches='tight' )数据分析主成分分析主成分分析（Prin­ci­pal Com­po­nent Analy­sis, PCA）的具体逻辑及Python实现方法可以看：PCA主成分析原理、理解和代码实现。聚类分析K-means聚类与层次聚类相关代码直接查看GitHub项目Data-Analysis-for-Thesis。Mantel Test见博文：Mantel Test算法原理讲解、代码实现、绘图可视化。
Tips / 提示这篇博文将是你可以在所有中文搜索引擎里找得到的最详细的关于Mantel Test的中文分析讲解。起源在一次课题组师兄汇报的时候，我第一听说了Mantel Test，当时第一眼就被这个漂亮的图形所吸引，所以就想着以后也能用到自己的文章里，便自己花时间了解了下。全网能找到这张图的资料，最早的时间是2015年，在Science中发表的一篇文章里第一次出现了这张图。SUNAGAWA S, COELHO L P, CHAFFRON S, et al. Structure and function of the global ocean microbiome[J]. Science, 2015, 348(6237): 1261359.只能说不愧能发Science，文章里面的每一张图都是经过精心处理的。介绍1967年，美国国立卫生研究院（National Institutes of Health）的生物统计学家Nathan Mantel提出了Mantel Test分析方法。在统计学中，传统相关系数只能用于计算分析一个数据矩阵中每两列变量之间的相关性，而在面对两个矩阵之间的相关性时就一筹莫展。Mantel Test的出现克服了传统相关系数的这一缺陷，该方法本质上是分析两个矩阵之间的相关性。自Mantel Test推出以来，其在微生物群落生态学、动植物种群遗传学等领域都得到快速发展与应用，在分析微生物菌群与环境因子之间的相关性方面更是得到了广泛的应用。在使用Mantel Test分析环境因子与微生物群落结构之间的相关性时，通常对微生物群落OTU数据矩阵使用Bray-Curtis相异度（Bray-Curtis dissimilarity）来计算微生物群落结构之间的差异性；而对环境变量数据矩阵采用欧氏距离（Euclidean Distances）进行差异性计算。Mantel Test的分析过程主要包括：分别使用各自的距离公式计算两个数据矩阵的距离矩阵，然后将两个距离矩阵进行压缩得到两个压缩距离列，然后计算这两列的相关性（一般都采用皮尔逊pearson相关性指数）；在完成一次计算后，对原数据矩阵中的一列或者两列进行置换，重新计算距离公式以及压缩距离公式，计算新的相关性系数（r值）；经过成千上万次的置换后，观察实际数据的r值在经过多次置换后所得的r值分布中的位置，如果跟随机置换得到的结果站队较近，则说明相关性不显著，相反则说明两个矩阵具有较强的显著相关性。讲解下面结合我查阅的一些中外文资料，对Mantel Test做出自己的理解，如果有错误的地方，恳请指正交流（文末有微信）。图形讲解右侧上三角首先来看图形右半部分，这部分大家都很常见，是一个相关性热图，它代表了一个数据矩阵中每两列之间的相关性。而计算相关性的算法一般都选择Pearson相关。皮尔逊（Pearson）相关（r），它测量两个变量（x和y）之间的线性相关性。它也称为参数相关性检验，因为它取决于数据的分布。仅当x和y来自正态分布时才可以使用它。y = f（x）的图称为线性回归曲线。举个例子：Oxygen与Temperature的相关性为深蓝色，对比最后边色条可以发现，深蓝色对应的相关性为-0.8~-1.0之间，说明其两者呈现出很强的相关性（相关性强弱看绝对值），且呈现负相关（相关性方向看正负）。左侧下三角再来看图形的左侧，左侧反映出来的也是相关性（r值），只不过除了相关性之外还反映出了一个显著性（p值），在相关性分析中，r值代表相关性，p（sig）值代表显著性。p值是判断r值相关系数是否具有统计学意义，判定标准一般为0.05，如果p值小于0.05说明两者之间具有较强的相关性，即r值有意义。然后再看r值，|r|值越大，相关性越好。正数指正相关，即一个变量随着另一个变量增大而增大；负数指负相关，即一个变量随另一个变量的增大而减小。可以看出右侧上三角中的r值与左侧下三角中的r值是一样的，具有相同的统计学概念。Pearson相关系数的一点小补充Pearson主要是把两个序列看成两个变量，计算两个变量是否线性相关。Pearson公式如：$$ p_{XY} = \frac{cov(X,Y)}{\sigma _{X} \sigma _{Y}}=\frac{E(XY)-EXEY}{\sigma _{X} \sigma _{Y}} $$分子是协方差的计算公式，协方差我们知道，是计算两个变量的相关性，大于0是正相关，小于是负相关，相等不相关。除以分母相当于归一化到[-1,1]之间。所以，Pearson相关系数的计算结果也等于将数据矩阵进行标准化后再求协方差，此时求出的协方差就等于源数据矩阵中各列的相关性。代码演示：import pandas as pd # 引包 import seaborn as sns df = sns.load_dataset('iris') # 加载机器学习经典数据集——鸢尾花数据集 df.info() # 查看数据信息 <class 'pandas.core.frame.DataFrame'> RangeIndex: 150 entries, 0 to 149 Data columns (total 5 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 sepal_length 150 non-null float64 1 sepal_width 150 non-null float64 2 petal_length 150 non-null float64 3 petal_width 150 non-null float64 4 species 150 non-null object # 这一列为花的分类，为Object文字对象，需要删除 dtypes: float64(4), object(1) memory usage: 6.0+ KB df.pop('species') # 删除分类列species文字列 df # 查看数据集 Out[14]: sepal_length sepal_width petal_length petal_width 0 5.1 3.5 1.4 0.2 1 4.9 3.0 1.4 0.2 .. ... ... ... ... [150 rows x 4 columns] X = (df - df.mean()) / df.std(ddof=0) # 将数据矩阵标准化 SIGMA = (X.T @ X) / X.shape[0] # 计算数据矩阵的协方差矩阵 SIGMA Out[20]: sepal_length sepal_width petal_length petal_width sepal_length 1.000000 -0.117570 0.871754 0.817941 sepal_width -0.117570 1.000000 -0.428440 -0.366126 petal_length 0.871754 -0.428440 1.000000 0.962865 petal_width 0.817941 -0.366126 0.962865 1.000000 df.corr() # 计算数据矩阵中每两列之间的Pearson相关性 Out[21]: sepal_length sepal_width petal_length petal_width sepal_length 1.000000 -0.117570 0.871754 0.817941 sepal_width -0.117570 1.000000 -0.428440 -0.366126 petal_length 0.871754 -0.428440 1.000000 0.962865 petal_width 0.817941 -0.366126 0.962865 1.000000由上述代码可以发现，一个数据矩阵的Pearson相关性系数矩阵就等于其协方差矩阵。算法讲解Mantel Test算法的讲解使用Python语言实现，编程语言只是工具，与算法概念无关，选择Python只是因为这门语言我比较顺手~数据准备正如上面的讲解所说，Mantel Test是对两个数据矩阵之间的相关性进行分析。因此，需要先准备好两个数据矩阵。举个栗子，我想要对一个微生物数据矩阵、一个环境因子数据矩阵进行分析：微生物OTU矩阵环境因子矩阵注意看，上方就是两个进行检验的矩阵，需要注意的是，两个数据矩阵的行索引应该是相互对应的。计算距离矩阵数据准备好之后，要进行的第一步操作便是计算两个数据矩阵各自的距离矩阵。一般情况下，对微生物或者基因数据矩阵采用bc距离公式计算，而对于环境因子数据矩阵采用欧氏距离公式计算。bc距离计算计算公式：$$ D_{Bray-Curtis} = 1-2 \frac{\sum min(S_{A,i}, S_{B,i})}{\sum S_{A,i}+ \sum S_{B,i}} $$详细计算过程：$min(S_{A,i}, S_{B,i})$ = 7+3+4+1+0 = 15$\sum S_{A,i}$ = 10+8+4+1+1 = 24$\sum S_{B,i}$ = 7+3+4+8+4+0 = 22D = 1 - 2*15/(24+22) = 0.3478D值越小表示二者组成差异小。也可以直接使用Python科学计算库scipy包中的工具进行计算：from scipy.spatial.distance import braycurtis as bc import pandas as pd OTU = pd.DataFrame([['A', 10, 8, 4, 1, 1], ['B', 7, 3, 8, 4, 0]], columns=['community', 'OTU1', 'OTU2', 'OTU3', 'OTU4', 'OTU5']) OTU.set_index('community', inplace=True) OTU Out[4]: OTU1 OTU2 OTU3 OTU4 OTU5 community A 10 8 4 1 1 B 7 3 8 4 0 bc(OTU.iloc[0,:], OTU.iloc[1,:]) Out[5]: 0.34782608695652173欧式距离欧式距离就不用介绍了，就是初中就学过的那个欧氏距离，直接自己写下就行了。非冗余距离公式经过计算后，得到的距离公式应该是这样的：这个时候会发现，其实这个距离公式是有瑕疵的，为啥这么说呢？可以看到1→2的距离其实是与2→1的距离一样的，且1→1、2→2的距离都是0。也就是说，这个距离公式中有一半的元素（上三角或者下三角）以及对角线上的元素其实都是无效的数据，因为我们根本不需要它或者用不到。这个时候我们就需要对距离公式进行压缩，使其成为一个非冗余距离矩阵，应该为：使其上三角或者下三角全为0，这个时候，这个距离矩阵才是非冗余的。自己封装了一个可以计算欧式距离或者bc距离公式，并可以得到冗余距离矩阵以及非冗余距离矩阵的函数：def calc_distance(df, method='euc', redundant=True): """ 计算欧式距离(Euclidean)和BC(Bray-Curtis)距离 :param df: 输入矩阵 :param method: euc：欧氏距离（默认值）；bc：Bray-Curtis距离；Others：Undefined :param redundant: 是否返回冗余距离矩阵？冗余矩阵就是上三角和下三角一样的一个矩阵 :return: 返回距离公式 """ method = method.lower() distance_df = np.zeros((df.shape[0], df.shape[0])) # 创建一个用于存储距离的矩阵 if redundant: # 计算冗余矩阵 for i in range(len(df)): # 遍历矩阵元素 for j in range(len(df)): # 由于距离矩阵中，对角线代表每个元素到自身的距离（必定为0），且上三角和下三角是重复的，因此精简算法，只对下三角进行填充 if method == 'euc': # 欧氏距离 distance_df[i, j] = np.sqrt(np.sum((df.iloc[i] - df.iloc[j]) ** 2)) elif method == 'bc': # 计算BC距离 # distance_df[i, j] = bc(df.iloc[i], df.iloc[j]) # 引包计算 _ = pd.DataFrame([df.iloc[i], df.iloc[j]]) distance_df[i, j] = 1 - 2 * np.sum(_.min()) / np.sum(np.sum(_)) else: # 报错：未定义的距离计算公式！ raise ValueError(f'Undefined distance formula: {method}!') else: # 计算非冗余矩阵（算法都一样） for i in range(1, len(df)): # 遍历矩阵元素 for j in range(i): # 由于距离矩阵中，对角线代表每个元素到自身的距离（必定为0），且上三角和下三角是重复的，因此精简算法，只对下三角进行填充 if method == 'euc': # 欧氏距离 distance_df[i, j] = np.sqrt(np.sum((df.iloc[i] - df.iloc[j]) ** 2)) elif method == 'bc': # 计算BC距离 distance_df[i, j] = bc(df.iloc[i], df.iloc[j]) else: # 报错：未定义的距离计算公式！ raise ValueError(f'Undefined distance formula: {method}!') return pd.DataFrame(distance_df)距离向量其实，上一步得到的非冗余距离矩阵还是可以再次压缩，使之成为向量的格式：这样的话，对两个数据矩阵计算得到的距离矩阵经过一系列的操作后，压缩成为了两个距离向量，这个时候，才完成了距离数据的准备。在上一步的距离矩阵上，可以得到距离向量的函数。def condensed_distance_matrices(df): """ 输入一个距离矩阵，对该矩阵进行压缩 :param df: 一个下三角距离矩阵 :return: 返回压缩后的距离数据（Series） """ pointer = [] # 存放距离指向 value = [] # 存放距离 for i in range(1, df.shape[0]): # 遍历矩阵元素 for j in range(i): pointer.append(str(f'{i}➡{j}')) value.append(df.iloc[i, j]) return pd.DataFrame({'Pointer': pointer, 'Value': value})相关性的计算Mantel Test所计算的相关性，其本质上是对两个数据矩阵的距离矩阵相关性进行检验。在得到两个距离向量之后，就可以使用Pearson相关系数来反应两个距离向量之间的相关性。而且我本人经过检验（R、Python），确实是这么一回事儿。显著性的计算按照上述方法，我们已经得到了相关性r值。但是，这个r值是否可靠还需要显著性p值来进行检验。而想要得到p值，就先必须理解统计学里的一个假设-检验（hypothesis testing）的概念。其基本思路就是：我想要证明数据矩阵A与数据矩阵B具有相关性；那么我先假设数据矩阵A与数据矩阵B之间没有相关性，预先设定的检验水准为0.05，当检验假设为真，但被错误地拒绝的概率，记作α，通常取α=0.05或α=0.01；对样本观察值按相应的公式计算出统计量的大小（例如T值等）；根据统计量的大小及其分布确定检验假设成立的可能性P的大小并判断结果。而这个假设检验体现在Mantel Test里面就是，按照上述得到r值的方法，对两个源数据矩阵的两行或者两列进行置换，再次计算出一个r值。就这样，进行千千万万次的置换，就可以得到千千万万个r值。这个得到的千千万万个r值，理论上是服从正态分布的。如果我们第一次计算的r值，能够落在接受原假设的区域上，那就说明原假设为真，即两个数据矩阵之间没有相关性；如果r值落在了拒绝原假设的域上，就说明可以拒绝原假设，即两个数据矩阵之间存在相关性。举个例子：选定检验水准α为0.05，如果第一次计算得到的r值为0.1516028059，且该r值在符合正态分布的经过千千万万次置换得到的r值list中可以落到拒绝原假设的域上，那就说明两个数据矩阵之间存在相关性。至于是单边检验还是双边检验，看了下源代码：@property def p(self): if self._p is None: if self.tail == "upper": # 右侧检验 self._p = sum(self.correlations >= self.r) / self.perms # self.perms就是千千万万次置换得到的r值得list的长度 elif self.tail == "lower": # 左侧检验 self._p = sum(self.correlations <= self.r) / self.perms else: # 双边检验 self._p = sum(abs(self.correlations) >= abs(self.r)) / self.perms return self._pR语言计算安装包# 安装包 ## install.packages("devtools") # 如果没有安装开发者工具箱的话，需要先安装 devtools::install_github("Hy4m/linkET", force = TRUE) 加载包# 加载包 library(linkET) library(ggplot2)加载数据集并进行Mantel检验# Mantel检验结果 data("varechem", package = "vegan") data("varespec", package = "vegan") mantel_test(varespec, varechem) # A tibble: 14 x 4 # 输出结果 spec env r p * <chr> <chr> <dbl> <dbl> 1 spec N 0.152 0.029 2 spec P 0.221 0.008 3 spec K 0.213 0.009 4 spec Ca 0.0961 0.138 5 spec Mg 0.0634 0.22 6 spec S 0.137 0.042 7 spec Al 0.278 0.004 8 spec Fe 0.167 0.045 9 spec Mn 0.221 0.014 10 spec Zn 0.0400 0.33 11 spec Mo 0.00330 0.48 12 spec Baresoil 0.145 0.045 13 spec Humdepth 0.188 0.023 14 spec pH -0.0282 0.573举例说明：微生物矩阵与环境因子矩阵中N元素相关性的显著性p值为0.029<0.05，说明其显著性较强，即相关性具有统计学意义；r值为0.152，说明两者呈现出相关性，且为正相关。绘制Mantel Test相关性图# Mantel绘图 library(dplyr) library(ggplot2) mantel <- mantel_test(varespec, varechem, spec_select = list(Spec01 = 1:7, Spec02 = 8:18, Spec03 = 19:37, Spec04 = 38:44)) %>% mutate(rd = cut(r, breaks = c(-Inf, 0.2, 0.4, Inf), labels = c("< 0.2", "0.2 - 0.4", ">= 0.4")), pd = cut(p, breaks = c(-Inf, 0.01, 0.05, Inf), labels = c("< 0.01", "0.01 - 0.05", ">= 0.05"))) qcorrplot(correlate(varechem), type = "lower", diag = FALSE) + geom_square() + geom_couple(aes(colour = pd, size = rd), data = mantel, curvature = 0.1) + scale_fill_gradientn(colours = RColorBrewer::brewer.pal(11, "RdBu")) + scale_size_manual(values = c(0.5, 1, 2)) + scale_colour_manual(values = color_pal(3)) + guides(size = guide_legend(title = "Mantel's r", override.aes = list(colour = "grey35"), order = 2), colour = guide_legend(title = "Mantel's p", override.aes = list(size = 3), order = 1), fill = guide_colorbar(title = "Pearson's r", order = 3))
Matplotlib绘图相关中文与负数中负号-的显示：plt.rcParams['font.sans-serif'] = ['SimHei'] # 中文显示 plt.rcParams['axes.unicode_minus'] = False # 负号显示当同时绘制的两组数据差距过大时，可以使用：plt.yscale('log') # 对y轴进行log缩放，减小两组图视觉上的差异，更方便观察（★★★★★）进行对数转换，减小不同组数据之间的量级差异。处理前处理后数据透视表三个分组统计函数size()、sum()、count()三个统计函数，能分别统计分组数量、不同列的分组和、不同列的分组数量。举个栗子：import pandas as pd df = pd.DataFrame({ 'key1': ['a', 'a', 'b', 'b', 'a'], 'key2': ['one', 'two', 'one', 'two', 'one'], 'x': np.random.randn(5), 'y': np.random.randn(5)}) df Out[4]: key1 key2 x y 0 a one -0.343202 0.992670 1 a two -0.996222 -0.181044 2 b one 1.380729 1.091965 3 b two 0.520865 -1.400465 4 a one -0.635766 -1.154578count()首先是count()函数，这个函数真的就是数数的，类似于pd.Series.value_counts()。df.groupby('key1').count() # key1列下面有3个a，2个b Out[5]: key2 x y key1 a 3 3 3 b 2 2 2 df.groupby(['key1', 'key2']).count() # key1列和key2列的values共有四种随机组合：a-one、a-two、：b-one、b-two，然后数一数他们的个数。 Out[6]: x y key1 key2 a one 2 2 two 1 1 b one 1 1 two 1 1 df.key1.value_counts() # 等同于pd.Series.value_counts() Out[7]: a 3 b 2 Name: key1, dtype: int64​ 所以综上所述，分组聚合时，使用count()函数好像根本没有意义~sum()不同于无聊的count()函数，sum()函数在分组聚合中感觉还是有点用的：df.groupby('key1').sum() # 将df按照key1列的value分份儿，然后求和。 Out[8]: # x列对应key1=a的value和为-1.975191，y列对应key1=b的value和为-0.308500 x y key1 a -1.975191 -0.342953 b 1.901594 -0.308500 df.groupby(['key1', 'key2']).sum() # 将df按照key1、key2列value的随机组合分份儿，然后求和。 Out[9]: # x列对应key1=a、key2=one的value和为-0.978968，y列对应key1=b、key2=two的value和为-1.400465 x y key1 key2 a one -0.978968 -0.161908 two -0.996222 -0.181044 b one 1.380729 1.091965 two 0.520865 -1.400465size()size()函数，顾名思义，就是求大小，这个大小不是指数值上的大小，而是个数上的大小。df.groupby('key1').size() # 先看key1列有几个value，然后看每个value有几行。 Out[10]: key1 a 3 b 2 dtype: int64 df.groupby(['key1', 'key2']).size() # 先看key1列与key2value的随机组合有几个，然后看每个组合有几行。 Out[11]: key1 key2 a one 2 two 1 b one 1 two 1 dtype: int64修改行索引列名称的方法# 方法一 df.rename_axis('Index', inplace=True) # 方法二 df.index.name = 'Index'提取时间格式列中的月份In [8]: demo = pd.DataFrame({ ...: 'date': ['2022-07-23', '2021-07-23'] ...: }) ...: demo['date'] = pd.to_datetime(demo['date']) ...: demo['astype'] = demo['date'].astype('datetime64[M]') ...: demo['dt.month'] = demo['date'].dt.month ...: demo Out[8]: date astype dt.month 0 2022-07-23 2022-07-01 7 1 2021-07-23 2021-07-01 7注意看两种提取方式的不同，一种可以把不同年份的月份都单独提取出来，另一种只能提取出来月份的数字。验证三列数据中对应行是否都是NaNnan_01 = pd.isna(df['Coupon_id']) nan_02 = pd.isna(df['Discount_rate']) nan_03 = pd.isna(df['Date_received']) np.all(nan_01 == nan_03) and np.all(nan_01 == nan_02) # 说明df['Coupon_id']列与df['Date_received']列的空值完全相等
预准备引包import numpy as np import pandas as pd import matplotlib.pyplot as plt绘图全局设置# 行内绘图 # % matplotlib inline plt.style.use('ggplot') # 绘图风格 plt.rcParams['font.sans-serif'] = ['SimHei'] # 显示中文字体读取数据并整体观察数据结构In [3]: col = ['user', 'date', 'product', 'amount'] ...: df_0 = pd.read_table('./data/CDNOW_master.txt', names=col, sep='\s+') # sep='\s+ 列于列之间的分隔为至少一个空格 ...: df = df_0.copy() ...: df Out[3]: user date product amount 0 1 19970101 1 11.77 1 2 19970112 1 12.00 2 2 19970112 5 77.00 3 3 19970102 2 20.76 4 3 19970330 2 20.76 ... In [4]: df.info() <class 'pandas.core.frame.DataFrame'> RangeIndex: 69659 entries, 0 to 69658 Data columns (total 4 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 user 69659 non-null int64 1 date 69659 non-null int64 2 product 69659 non-null int64 3 amount 69659 non-null float64 dtypes: float64(1), int64(3) memory usage: 2.1 MB分析：观察数据结构可得，数据表一共有69659行、4列，且无空行，其中4列包括user（用户索引列）、date（用户消费时间列）、product（购买产品数量列）、amount（购买总金额列）。数据的预处理观察date（用户消费时间列）可发现，其时间格式Pandas未能识别，需要手动将其转换成时间格式列（datetime），方便后续操作。In [5]: df['date'] = pd.to_datetime(df['date'], format='%Y%m%d')后续数据分析需要按月来操作，因此需要读取date（用户消费时间列）中的月份，In [6]: df['month'] = df['date'].astype('datetime64[M]')注意看这里这里从date（用户消费时间列）中获取月份的方式，并没有使用：In [7]: df['date'].dt.month为什么呢？使用.dt.month虽然可以获取月份，但是遇到多个年份中提取月份的话就会出现问题，做个简单小对比：In [8]: demo = pd.DataFrame({ ...: 'date': ['2022-07-23', '2021-07-23'] ...: }) ...: demo['date'] = pd.to_datetime(demo['date']) ...: demo['astype'] = demo['date'].astype('datetime64[M]') ...: demo['dt.month'] = demo['date'].dt.month ...: demo Out[8]: date astype dt.month 0 2022-07-23 2022-07-01 7 1 2021-07-23 2021-07-01 7由此可以看出，当遇到多个年份有相同月份的时候，使用.dt.month仅仅可以提取出来月份，但是具体是哪一年的7月，确提取不出来，因此其存在着弊端。数据分析月统计量分析按月份统计产品购买数量、消费金额、消费次数以及消费人数。df.pivot_table( index='month', aggfunc={ 'user': 'count', # 每个月的顾客数量（当同一个顾客下多次订单时，都按照新顾客统计） 'product': 'count', 'amount': 'sum' })注意：这种聚合分析方法存在的问题就是，如果一个用户一个月内有多条消费记录，那么其每次都会被当作一个新的消费者记录，即通过'user': 'count'聚合得到的是一个月所有消费记录的数量，并不是本月内有多少不同的顾客来过该店里。那么如何实现统计每个月内有多少顾客（无重复，比如一个顾客一个月内来了10次店里消费，也按成一次处理）来过店里呢，这就需要分组后做一次去重操作。df.pivot_table( index='month', aggfunc={ 'user': lambda x: pd.Series.nunique(x), 'product': lambda x: x.nunique(), 'amount': 'sum' })使用groupby()方法，然后配合聚合函数也可以达到相同的效果：df.groupby('month')'user', 'product', 'amount'.agg({ 'user': lambda x: pd.Series.nunique(x), 'product': lambda x: pd.Series.nunique(x), 'amount': 'sum' })订单商品数量与订单金额关系分析df'product', 'amount'.plot(kind='scatter', x='product', y='amount')由图中可以看出，大部分订单都有趋向于购买价格少、数量少的特征。用户消费金额占比（贡献度）重点使用.cumsum()函数。user_cumsum = df.pivot_table(index='user', values='amount', aggfunc='sum') user_cumsum = user_cumsum.sort_values(by='amount').reset_index() user_cumsum['cumsum_per'] = (user_cumsum['amount'].cumsum() / np.sum(user_cumsum['amount'])).apply( lambda x: format(x, '.2%')) user_cumsum.tail(10) Out[9]: # 消费累计中贡献度最大的十位用户 user amount cumsum_per 23560 15162 4234.45 97.56% 23561 3049 4262.85 97.73% 23562 499 4378.55 97.90% 23563 22279 4490.64 98.08% 23564 6569 4968.00 98.28% 23565 7931 6497.18 98.54% 23566 19339 6552.70 98.80% 23567 7983 6973.07 99.08% 23568 14048 8976.33 99.44% 23569 7592 13990.93 100.00%用户购买行为分析每天客流量分析df.pivot_table(index='user', values='date', aggfunc='min').value_counts().sort_index().plot()用户最后一次消费时间df.pivot_table(index='user', values='date', aggfunc='max').value_counts().sort_index().plot()由此图可以发现，在门店开业初期吸引了大量用户，但是大部分用户仅在前期参与活动后，后期便再也没来过。RFM模型最近一次消费 (Recency)——date列最大值消费频率 (Frequency)——product购买商品数量代替消费金额 (Monetary)——amount消费金额rfm_model = df.pivot_table( index='user', # values=['product', 'amount', 'date'], aggfunc={ 'product': 'size', 'amount': 'sum', 'date': 'max' } ) # 求出R rfm_model['R'] = df['date'].max() - rfm_model['date'] rfm_model['R'] = rfm_model['R'].dt.days # rfm_model.rename(columns={ 'amount': 'M', 'product': 'F' }, inplace=True) rfm_model.drop('date', axis=1, inplace=True)根据RFM模型对用户进行分类：def check_cus_grade(s): label = s.apply(lambda x: '1' if x >= 1 else '0') label = label['R'] + label['F'] + label['M'] grade = { '111': '重要价值客户', '011': '重要保持客户', '101': '重要发展客户', '001': '重要挽留客户', '110': '一般价值客户', '010': '一般保持客户', '100': '一般发展客户', '000': '一般挽留客户' } return grade[label] rfm_model['label'] = rfm_model.apply(lambda x: x - x.mean()).apply(check_cus_grade, axis=1) rfm_modelRFM模型的可视化：用户生命周期分析所谓用户生命周期是指用户第一次消费与最后一次消费的时间间隔。用户最后一次购买日期\=\=第一次购买的日期，说明用户仅仅购买了一次或者用户在同一天内购买了两次。这个时候就需要排除掉仅由一条消费记录的顾客，如果用户仅有一条购买数据，那还谈何生命周期，直接赋值为NaN。三种不同的实现方法：# 1、 df.groupby('user').apply(lambda X: X['date'].max() - X['date'].min() if X.shape[0] > 1 else np.NaN) # 2、 df.groupby('user')['date'].agg(lambda s: (s.max() - s.min()) if len(s) > 1 else np.NaN)是否多次消费df.groupby('user').apply(lambda X: 1 if X.shape[0] > 1 else 0).value_counts(normalize=True).plot.pie(autopct='%1.1f%%') plt.legend(['仅消费一次', '多次消费'])生命周期分析user_life_cycle = df.pivot_table(index='user', values='date', aggfunc=(lambda s: (s.max() - s.min()) if len(s) > 1 else np.NaN)) user_life_cycle.div(np.timedelta64(1, 'D')).plot.hist() plt.title('所有用户生命周期直方图') plt.xlabel('生命周期天数') plt.ylabel('用户人数')复购率与回购率分析复购率计算方式：在自然月内，购买多次的用户在总消费人数中的占比（若客户在同一天消费了多次，也称之复购用户）。user_repurchase_rate = df.pivot_table(index='user', columns='month', values='product', aggfunc='count') #复购用户：1；非复购的消费用户: 0；自然月没有消费记录的用户: NAN user_repurchase_rate = user_repurchase_rate.applymap(lambda num: 1 if num > 1 else 0 if num == 1 else np.NaN)回购率计算方式：在一个时间窗口（一个月）内进行了消费，在下一个窗口内又进行了消费。回购用户：本月消费了，上个月也消费了，记为1；非回购用户：本月消费了，但是上个月未消费，记为0；本月未消费，记为NaN。def back_purchase(s): every_month_status = [] # 第一个月不会有回购用户 if s[0] == 1: # 第一个月消费了 every_month_status.append(0) # 记为非回购用户 else: # 第一个月未消费 every_month_status.append(np.NaN) # 记为NaN for i in range(1, len(s)): if s[i] == 1: # 当前月份消费了 if s[i - 1] == 1: # 上个月份也消费了 every_month_status.append(1) # 回购用户 else: # 下个月份未消费 every_month_status.append(0) # 非回购用户 else: # 当前月份未消费 every_month_status.append(np.NaN) return pd.Series(every_month_status, user_repurchase_rate.columns) purchase_return = user_repurchase_rate.apply(back_purchase, axis=1)复购率与回购率的可视化purchase_return.sum().div(purchase_return.count()).plot(label='回购率') user_repurchase_rate.sum().div(user_repurchase_rate.count()).plot(label='复购率') plt.legend() plt.ylabel('百分比(%)') plt.title('用户复购率与回购率对比图')plt.plot(purchase_return.sum(), label='回购人数') plt.plot(user_repurchase_rate.sum(), label='复购人数') plt.legend() plt.ylabel('百分比(%)') plt.xlabel('month') plt.title('用户复购人数与回购人数对比图')总结1、用户个体特征：每笔订单的金额和商品购买量都集中在区间的低段水平，都是小金额小批量进行购买，此类交易群体，可在丰富产品线和增加促销活动提高转换率和购买率。2、大部分用户的消费总额和购买总量都集中刚在低段，长尾分布，这个跟用户需求有关，可以对商品进行多元文化价值的赋予，增强其社交价值属性，提高用户的价值需求。3、用户的消费周期：有二次以上消费的用户，平均68天，所以在50天到60天期间，应该对这批用户进行刺激召回，细致点，比如10天回复满意度，30天发放优惠券，55天的时候提醒优惠券的使用。4、用户的生命周期：有二次及以上消费的用户的平均生命周期是276天。用户的生命周期分别在20天内与400至500天间，应该在20天内对客户进行引导，促进其再次消费并形成消费习惯，延长其生命周期；在100至400天的用户，也要根据其特点推出有针对性的营销活动，引导其持续消费。5、新客户的复购率约为12%，老客户的复购率在20%左右；新客户的回购率在15%左右，老客户的回购率在30%左右，需要营销策略积极引导其再次消费及持续消费。6、用户质量：用户个体消费有一定规律性，大部分用户的消费集中在2000以下，用户消费反应了2/8法则，消费排名前20%的用户贡献了80%的消费额。所以说，狠抓高质量用户是万古不变的道理，这些高质量客户都是“会员”类型，需要专门]为会员优化购物体验，比如专线接听、特殊优惠等等。
在Python中，我们可以使用map()函数对list对象中的每一个元素进行循环迭代操作，例如：In [1]: a = [i for i in range(10)] In [2]: a Out[2]: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] In [3]: list(map(lambda x: x**2 ,a)) # 对list对象a中的每一个元素都进行计算平方值。map() 函数生成的是一个map对象，需要使用list()函数对其强制转换为list对象才可以。 Out[3]: [0, 1, 4, 9, 16, 25, 36, 49, 64, 81] 那么在Pandas操作中，有没有类似的功能可以实现对矩阵或者向量进行操作呢？当时是有的，这篇笔记来汇总下自己了解的几种方法。apply()在Pandas中，无论是矩阵（DataFrame）或者是向量（Series）对象都是有apply()方法的。对DataFrame对象使用该方法的话就是对矩阵中的每一行或者每一列进行遍历操作（通过axis参数来确定是行遍历还是列遍历）；对Series对象使用该方法的话，就是对Series中的每一个元素进行循环遍历操作。对DataFrame对象使用apply()方法：In [4]: import pandas as pd In [5]: df = pd.DataFrame({'a': [10, 20, 30], 'b': [20, 30, 40]}) In [6]: df Out[6]: a b 0 10 20 1 20 30 2 30 40 In [7]: df.apply(sum) # 对df中每一列的Series使用sum函数 Out[7]: a 60 b 90 dtype: int64 In [10]: df.apply(lambda s: s.min(), axis=1) # 对df中的每一行Series使用.min()方法，axis=1设置对df中的行进行操作 Out[10]: 0 10 1 20 2 30 dtype: int64对Series对象使用apply()方法：In [12]: type(df['a']) Out[12]: pandas.core.series.Series In [14]: df['a'].apply(lambda x: x * 20) Out[14]: 0 200 1 400 2 600 Name: a, dtype: int64当apply()配合有参数的函数时，也可以为函数传递参数：In [15]: def calc_n(x, n): ...: return x ** n ...: In [16]: df['a'].apply(calc_n, n=10) Out[16]: 0 10000000000 1 10240000000000 2 590490000000000 Name: a, dtype: int64applymap()矩阵（DataFrame）的applymap()方法可以对矩阵中每一个元素进行遍历迭代操作：In [18]: df.applymap(lambda x: x * 2) Out[18]: a b 0 20 40 1 40 60 2 60 80 In [19]: df*2 Out[19]: a b 0 20 40 1 40 60 2 60 80行、列的迭代除了对矩阵使用apply()方法进行迭代外，还可以.iteritems()、.iterrows()与.itertuples()方法进行行、列的迭代，以便进行更复杂的操作。.iteritems()列迭代每次取出的i是一个元组，在元组中，第[0]项是原来的列名称，第[1]列是由原来该列的元素构成的一个Series：In [20]: for i in df.iteritems(): ...: print(type(i), i) ...: <class 'tuple'> ('a', 0 10 1 20 2 30 Name: a, dtype: int64) <class 'tuple'> ('b', 0 20 1 30 2 40 Name: b, dtype: int64) In [21]: for i in df.iteritems(): ...: print(type(i[1]), i[1]) ...: <class 'pandas.core.series.Series'> 0 10 1 20 2 30 Name: a, dtype: int64 <class 'pandas.core.series.Series'> 0 20 1 30 2 40 Name: b, dtype: int64行迭代与列迭代的形式一样：In [22]: for i in df.iterrows(): ...: print(type(i), i , end='\n --------- \n') ...: <class 'tuple'> (0, a 10 b 20 Name: 0, dtype: int64) --------- <class 'tuple'> (1, a 20 b 30 Name: 1, dtype: int64) --------- <class 'tuple'> (2, a 30 b 40 Name: 2, dtype: int64) --------- In [23]: for i in df.iterrows(): ...: print(type(i[1]), i[1] , end='\n --------- \n') ...: <class 'pandas.core.series.Series'> a 10 b 20 Name: 0, dtype: int64 --------- <class 'pandas.core.series.Series'> a 20 b 30 Name: 1, dtype: int64 --------- <class 'pandas.core.series.Series'> a 30 b 40 Name: 2, dtype: int64 --------- .itertuples()方法取出的每一行是一个Pandas对象：In [24]: for i in df.itertuples(): ...: print(type(i), i) ...: <class 'pandas.core.frame.Pandas'> Pandas(Index=0, a=10, b=20) <class 'pandas.core.frame.Pandas'> Pandas(Index=1, a=20, b=30) <class 'pandas.core.frame.Pandas'> Pandas(Index=2, a=30, b=40)函数向量化Series是一个向量，但是其中的元素却是一个个数值，如何将两个Series像两个数值元素一样进行使用？In [27]: import numpy as np In [28]: def vs_num(a, b): ...: if a > b: ...: return 'a获胜！' ...: elif a == b: ...: return '平局~' ...: else: ...: return 'b获胜！' ...: In [29]: np.vectorize(vs_num)(df['a'], df['b']) Out[29]: array(['b获胜！', 'b获胜！', 'b获胜！'], dtype='<U4') In [30]: df Out[30]: a b 0 10 20 1 20 30 2 30 40在函数前加上@np.vectorize语法糖也是同样的效果。
产生需求任何一个功能的出现，必定是有其对应的需求的。使用Python获取斐波拉契数列，斐波拉契数列（Fibonacci）：除第一个数以外，所有的数字都是由前面两项相加得到。def fib(n): """ 获取斐波拉契数列 :param n: 需要获取自然数n前的所有斐波拉契数列中的数字 :return: result，以列表的形式返回数列 """ a, b = 0, 1 result = [a, b] while True: c = a + b if c <= n: result.append(c) a, b = b, c else: break return result a = fib(100) # [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89]问题通过上面的函数我们可以获取数字n前面的所有斐波拉契数列，但是如果我使用函数fib()获取了100前的所有斐波拉契数列，但是我再使用的过程中却发现我只需要10以前的斐波拉契数列就行了。这个时候，就产生了内存浪费，制造了100前的所有斐波拉契数列，但是你实际需求只用了10以前的。抛出疑问能不能写出这样一个函数，就是我所产生的斐波拉契数列个数是小于等于100前的，但是至于我用到哪些，现在是还不确定的，但是我可以随用随造，不用10以后的斐波拉契数列的话，那么这个函数就不会制造后面的呢？解决办法这个时候，就用到了生成器。结合列表推导式，把其两端的[]换成()即可得到生成器，一个简单的举例：a = (i for i in range(3))与列表推导式不同的是，列表生成器在创建的时候直接生成了一个列表：[0, 1, 2]，而生成器则是创建了一个生成器对象，只有在后续调用的时候才会生成其中的元素。b = [i for i in range(3)] print(a) # <generator object <genexpr> at 0x000001A7767A5690> print(b) # [0, 1, 2]如何使用生成器呢？生成器每次使用都需要进行调用，且生成的值只能使用一次，之后就产生下一个元素：print(a.__next__()) # 0 print(a.__next__()) # 1 print(next(a)) # 2 print(a.__next__()) # 报错：StopIteration如果超出生成器中元素个数，就会报错。上述逐个调用生成器中的元素的话，如果超出则会报错，但是使用for循环的话就不会报错：for i in a: print(i) #0, 1, 2其实Python中很多地方都有生成器的使用场景，例如使用open()函数对文件执行逐行读取的时候，就是每次读取一行，不读就不取。yield关键字回到最初提出的问题，如何实现函数中也实现生成器的效果呢？这个时候就需要使用yield关键字。def fib(n): a, b = 0, 1 i = 0 while True: i += 1 print(f'开始第{i}次循环...') c = a + b if c <= n: yield c a, b = b, c else: break print(f'第{i}次循环结束...') # 正常函数在创建的时候就会运行完毕，但是加入yield关键字的函数，在创建的时候，运行到yield关键字的位置就会停止，并在第一次调用.__next__()方法的时候将yield关键字后面的变量返回出来，然后继续执行后续代码，直至再次遇见yield。 a = fib(100) # [1, 2, 3, 5, 8, 13, 21, 34, 55, 89] while True: try: print(a.__next__()) except Exception as e: print('已打印完所有元素...', e) break
什么是列表推导式？当一个for循环中代码比较简单的时候，往往都可以写成列表推导式，这样代码更加简洁、高效。举个栗子：a = list(range(1, 11)) # [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] # 智障版 for index, value in enumerate(a): # a[index] = value + 1 a[index] += 1 a # [2, 3, 4, 5, 6, 7, 8, 9, 10, 11] # map()函数 list(map(lambda x: x + 1, a)) # [3, 4, 5, 6, 7, 8, 9, 10, 11, 12] # 列表生成式 [i+1 for i in a] #[3, 4, 5, 6, 7, 8, 9, 10, 11, 12]配合if else的高级使用仅配合if语句：[i ** 2 for i in a if i % 2 == 0] Out[5]: [0, 4, 16, 36, 64, 100]配合if else语句：首先复习一个三元表达式：b = 2 c = 1 if b >= 3 else 0 c Out[8]: 0然后再来看：[i ** 2 if i % 2 == 0 else i ** 3 for i in a ] Out[9]: [0, 1, 4, 27, 16, 125, 36, 343, 64, 729, 100]注意看，这个时候if else语句的位置已经在for语句之前了。配合if elif else语句：if elif else这个语句其实是不能实现的，但是可以绕一下，通过if else (if else)来实现（不加括号也是可以的，但是为了看起来更容易理解就加上了）：d = 2 if b >= 3 else (0 if b >= 0 else None) d Out[12]: 0理解了上述语句之后，再来看列表推导式就不难了：[i ** 2 if i % 2 == 0 else (i ** 3 if i % 3 == 0 else None) for i in a] Out[13]: [0, None, 4, 27, 16, None, 36, None, 64, 729, 100]配合函数的高级使用当需要对for寻循环中的每一个i都需要进行比较复杂操作的时候，可以通过封装函数的方法来实现。lambda匿名函数：[(lambda x: x ** x)(i) for i in a] # 计算列表a中每一个x数的x次方 Out[15]: [1, 1, 4, 27, 256, 3125, 46656, 823543, 16777216, 387420489, 10000000000]配合lambda匿名函数加上if语句：[(lambda x: x ** x)(i) for i in a if i % 2 == 0] # 计算列表a中每一个x数的x次方 Out[16]: [1, 4, 256, 46656, 16777216, 10000000000]配合自定义函数：def calc(n): return n ** n [(lambda x: x ** 0)(i) if i % 2 == 0 else calc(i) for i in a] Out[19]: [1, 1, 1, 27, 1, 3125, 1, 823543, 1, 387420489, 1]除了使用列表推导式+函数的方法来对列表中的每一个元素进行操作，还有一种简便的方法就是使用map()函数。# map对象 map(lambda x: x ** 2, [i for i in a if i % 2 == 0]) Out[21]: # 配合lambda匿名函数 list(map(lambda x: x ** 2, [i for i in a if i % 2 == 0])) Out[22]: [0, 4, 16, 36, 64, 100] # 配合一般函数 list(map(calc, [i for i in a if i % 2 == 0])) Out[24]: [1, 4, 256, 46656, 16777216, 10000000000]需要注意的是，map()函数返回值是一个map对象，需要使用list()强制转化。这里有个小疑问：Pandas里的apply()函数可以为其函数传递参数，但是map()函数好像没有这个功能。import pandas as pd def calc_n(x, n): return x ** n a_s = pd.Series(a) a_s.apply(calc_n, n = 10) # 函数calc_的参数也可以通过apply()函数进行传递 Out[32]: 0 0 1 1 2 1024 3 59049 4 1048576 5 9765625 6 60466176 7 282475249 8 1073741824 9 3486784401 10 10000000000 dtype: int64
今天在B站看了个视频，终于让我搞懂了什么是装饰器，解决了心里多年的一个困惑，记录下。产生的需求一个新闻门户网站，分为sports、news、user_center三个页面，想要实现用户必须登录才能看片的效果。源代码：user_info = { # 用户信息 'user_name': 'Yacan', 'password': '123' } def sports(): print("体育页面...") def news(): print("新闻页面...") def user_center(): print('用户中心...')简单粗暴的实现方法# 为用户增加一个is_login的属性，当用户登陆后，将其赋值为True user_info = { 'user_name': 'yacan', 'password': '123', 'is_login': False } # 增加一个check_login函数，用来检查用户的登录状态 def check_login(): name = input('请输入用户名：') password = input('请输入密码：') if name == user_info['user_name'] and password == user_info['password']: user_info['is_login'] = True print('登录成功~') else: print('账号或者密码错误，登陆失败！') # 在加载每一个板块内容前都先执行登录函数 def sports(): check_login() print("体育页面...") def news(): check_login() print("新闻页面...") def user_center(): check_login() print('用户中心...')问题：简单粗暴的方式违背了代码的“开放-封闭”原则，即当我们增加新功能时，要尽量地不在原代码上产生改动，尽量把所有的改动都写到函数外。升级版如果我们可以将需求函数sports()、news()、user_center()放入check_login()函数来执行的话，就可以不在源代码上进行修改，从而把全部的新增代码都放到了外面，通过Python函数的高级用法——把一个函数作为参数传入另一个函数内：def check_login(fun): name = input('请输入用户名：') password = input('请输入密码：') if name == user_info['user_name'] and password == user_info['password']: user_info['is_login'] = True print('登录成功~') fun() # 执行传入的函数 else: print('账号或者密码错误，登陆失败！') sports = check_login(sports) sports() news = check_login(news) news() user_center = check_login(user_center) user_center()问题：在执行sports = check_login(sports)的时候，由于check_login(fun)没有return（返回值），导致sports = None，后面执行sports()的时候，就是None()，直接报错。小知识点在Python中定义一个函数def fun()后，直接调用fun返回的是函数的内存地址，调用fun()返回的是函数的return值。def fun(): return 1 print(fun(), ' | ', fun) # 1 | <function fun at 0x000001FBF182EB90>完善想要解决上述问题，就是想办法让执行sports = check_login(sports)的时候返回的是函数的内存ID，就可以解决这个问题。这个时候就用到了一个伟大的概念——闭包。def check_login(fun): def inner(): # 所谓闭包，就是在函数里面加上一个嵌套函数，然后将内部嵌套函数返回出去。 name = input('请输入用户名：') password = input('请输入密码：') if name == user_info['user_name'] and password == user_info['password']: user_info['is_login'] = True print('登录成功~') fun() else: print('账号或者密码错误，登陆失败！') return inner news = check_login(news) sports = check_login(sports) news() sports()所谓闭包，就是在函数里面加上一个嵌套函数，然后将内部嵌套函数返回出去。执行sports = check_login(sports)的时候返回的是函数的内存ID，而不是None了。至此，装饰器起始已经产生了，但是为了书写方便，Python有一个专门的写法：@check_login def sports(): print("体育页面...") @check_login def news(): print("新闻页面...") @check_login def user_center(): print('用户中心...')就是在每一个函数前面的一行加上@标识符，在sports()函数前面一行加上@check_login的作用就是相当于执行了sports = check_login(sports)代码。完美，Python太优雅了。
举个🌰来说明通过一个例子来说明什么是中心化与标准化，并且如何进行标准化与中心化。1、Python包的加载与画图设置（不懂可以不看）：# 包的加载 import pandas as pd import numpy as np import seaborn as sns import matplotlib.pyplot as plt import palettable #python颜色库 # 用于显示中文 plt.rcParams['font.sans-serif'] = ['SimHei'] plt.rcParams['axes.unicode_minus'] = False2、加载数据：不懂代码的话可以理解为这里是读取一个表格数据。# 这里做一个对矩阵进行标准化与中心化的对比，阐述其意义 state_data_0 = pd.read_csv('./state.csv') state_data_0.set_index('Unnamed: 0', inplace=True) state_data_0.index._name = 'state' state_data = state_data_0.copy() state_data这个数据表是美国几个州的统计数据，每一行代表一个州，每一列分别是人口（Population）、收入（Income）、受教育程度（Illiteracy）、平均寿命（Life Exp）等等。表格数据如图所示3、使用表格数据画图：plt.figure(dpi=120) # 在进行标准化与中心化处理之前绘制热图，由于不同列的数据差距过大，在绘制热图的时候就会导致反应出来的效果不好（例如：某一列的数据都在100~900，而另一列的数据都在10~99） # sns.heatmap(state_data) sns.heatmap(data=state_data, cmap=sns.cubehelix_palette(as_cmap=True), cbar=True, cbar_kws={'label': 'ColorBar', #color bar的名称 'orientation': 'vertical', #color bar的方向设置，默认为'vertical'，可水平显示'horizontal' # "ticks": np.arange(4.5, 8, 0.5), #color bar中刻度值范围和间隔 # "format": "%.2f", #格式化输出color bar中刻度值 "pad": 0.05, #color bar与热图之间距离，距离变大热图会被压缩 }, # annot=True, #默认为False，当为True时，在每个格子写入data中数据 # fmt=".2f",#设置每个格子中数据的格式，参考之前的文章，此处保留两位小数 # annot_kws={'size':4,'weight':'normal', 'color':'blue'},#设置格子中数据的大小、粗细、颜色 linewidths=1, #格子与格子，默认为0 # linecolor='red', #每个格子边框颜色,默认为白色 ) plt.title('未经过中心化与标准化的数据')出图4、遇到了问题：如上图所示，使用表格数据绘制了一张热图，但是发现由于表格中的数据相差过大（Population列都在千、万级别，Forst都在10~1000级别）导致绘制出来的图并不能很好地表达出来一定的区分度。这个时候就在想有没有一种方法能够解决这个问题呢？——数据的标准化与中心化，也被称为均值方差归一化（Standardization），把所有数据归一到均值为0方差为1的分布中，适用于数据分布没有明显边界的情况中。标准化与中心化其实就是对一组数据求平均值和方差，然后计算：$$ 处理后的结果=\frac{(源数据-源数据平均值)}{源数据的方差} $$Python代码实现：def norm_(pd_raw): """ 定义一个可以对DataFrame进行中心化与标准化的函数 :param x: DataFrame :return: 经过标准化的DataFrame """ pd_mean = np.mean(pd_raw, 0) # 求DataFrame每一列的平均值 pd_std = np.std(pd_raw, 0) # 求DataFrame每一列的标准差 return (pd_raw - pd_mean) / pd_std上面这个函数就是定义一个可以用于将数据（Python中的DataFrame对象）进行标准化与中心化的函数，不懂代码的话可以理解为这一步就是如何将数据进行标准化与中心化。5、对源数据进行标准化与中心化，然后重新出图：state_data_norm = norm_(state_data) # 这里对数据进行标准化与中心化处理，处理后各个列的数据都向中间靠拢 sns.heatmap(data=state_data_norm, cmap=sns.cubehelix_palette(as_cmap=True), cbar=True, cbar_kws={'label': 'ColorBar', #color bar的名称 'orientation': 'vertical', #color bar的方向设置，默认为'vertical'，可水平显示'horizontal' # "ticks": np.arange(4.5, 8, 0.5), #color bar中刻度值范围和间隔 # "format": "%.2f", #格式化输出color bar中刻度值 "pad": 0.05, #color bar与热图之间距离，距离变大热图会被压缩 }, # annot=True, #默认为False，当为True时，在每个格子写入data中数据 # fmt=".2f",#设置每个格子中数据的格式，参考之前的文章，此处保留两位小数 # annot_kws={'size':4,'weight':'normal', 'color':'blue'},#设置格子中数据的大小、粗细、颜色 linewidths=1, #格子与格子，默认为0 # linecolor='red', #每个格子边框颜色,默认为白色 ) plt.title('经过中心化与标准化的数据')经过中心化与标准化的数据可以看到，经过中心化与标准化后，可以在热图上很好地对级别相差较大的数据进行区分。这就是数据中心化与标准化的流程。最值归一化在有些文章中还会提到一个概念叫做最值归一化（Normalization），所谓最值归一化就是把一组数据映射到0~1之间：$$ x_{scale} = \frac{x-x{min}}{x_{max} - x_{min}} $$这个方法适用于分布有明显边界的情况，且受outlier的影响较大。scikit-learn中的封装方法MinMaxScaler在sklearn.preprocessing.MinMaxScaler中封装着最值归一化的方法，具体的使用方法可以参考官网链接。StandardScaler均值方差归一化的调用对象为sklearn.preprocessing.StandardScaler，下面是一个举例：# Import package from sklearn.preprocessing import StandardScaler # Create object scale_scaler = StandardScaler() # Input train data set scale_scaler.fit(x_train) # Get mean and std of train data set scale_scaler.mean_, scale_scaler.scale_ # Get standard train data and test data x_train_std = scale_scaler.transform(x_train) x_test_std = scale_scaler.transform(x_test)注意 / WARNING这里有一个值得注意的问题：在进行机器学习的过程中，将特征矩阵x分为训练集x_train和测试集x_test后，需要分别进行归一化操作。但是切记，这里并不是分别求出x_train和x_test的均值mean和方差std，然后再分别进行归一化；而是对x_test进行归一化也要采用x_train的均值和方差。还有一个非常重要的点是，如果仅对x_train进行了归一化，而x_test却没有归一化，那么机器学习预测结果将会非常非常差😡😡😡！！！别问我怎么知道的！都是泪！DIY一个import numpy as np class StandardScalerDiy: """ Do it yourself: sklearn.preprocessing.StandardScaler """ def __init__(self): self.mean_ = None self.scale_ = None def fit(self, x): assert x.ndim == 2, 'The dimension of x must be two.' self.mean_ = np.mean(x, axis=0) self.scale_ = np.std(x, axis=0) return self def transform(self, x): assert x.ndim == 2, 'The dimension of x must be two.' assert x.shape[1] == len(self.mean_), 'x must have same number of column with self.mean_.' return (x - self.mean_) / self.scale_