前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >开源图书《Python完全自学教程》12.4科学计算

开源图书《Python完全自学教程》12.4科学计算

作者头像
老齐
发布2022-12-09 20:23:44
1.4K0
发布2022-12-09 20:23:44
举报
文章被收录于专栏:老齐教室

12.4 科学计算

科学计算是科学、工程等项目中必不可少的,MATLAB 曾风光一时,但它是收费的,并且有“被禁”的风险——坚决反对用盗版软件,“被禁”不是盗版的理由。其实,Python ——开源、免费——是做科学计算的选择之一,它不仅能做 MATLAB 所能做的一切,还能做它不能做的。所以隆重推荐,在科学计算上选用 Python 。

12.4.1 Jupyter

Jupyter 是一款基于浏览器的开源的交互开发环境,常用于科学计算、数据科学、机器学习等业务中。其官方网站是:https://jupyter.org/ ,目前有 JupyterLab 和 Jupyter Notebook ,一般认为 JupyterLab 更趋近于 IDE,功能多于 Jupyter Notebook。下面以 JupyterLab 为例,演示安装和使用方法。

代码语言:javascript
复制
% pip install jupyterlab

使用 pip 安装,安装之后,执行如下指令启动 JupyterLab :

代码语言:javascript
复制
% jupyter-lab
... (省略部分显示信息)
    Or copy and paste one of these URLs:
        http://localhost:8888/?token=549de63a2fcd53325c5d17b8dd7a00ce4b63766aa368efb1
     or http://127.0.0.1:8888/?token=549de63a2fcd53325c5d17b8dd7a00ce4b63766aa368efb1

一般情况下,会自动打开本地计算机操作系统默认的浏览器,并显示图12-4-1所示效果。

图12-4-1 JupyterLab 默认页面

如果没有打开图12-4-1所示页面,或者打算用其他浏览器,可以复制终端所显示的地址(如:http://localhost:8888/?token=549de63a2fcd53325c5d17b8dd7a00ce4b63766aa368efb1 )到指定浏览器地址栏,同样能打开图12-4-1所示界面(注意,地址中的 token 项不能缺少)。

JupyterLab 所提供的各项功能,与通常 IDE 类似,此处不做详细说明,有意逐项了解的读者可以自行查阅专门资料。下面仅就如何使用它编写程序和执行所编写的代码给予说明,这是最基本的应用。

在图12-4-1所示界面的 Launcher 标签下,选择 Notebook 中的第一项“Python 3”(读者的开发环境可能与图中所示不同,只要选择“Python 3 ”作为程序执行的驱动即可。然后进入图12-4-2所示页面。

图12-4-2 编辑页面

参照图12-4-3所示的功能,将当前文件改名为 scicomputing.ipynb ,并将左侧文件栏折叠收起。此外,从图示中还可看到其他关于文件和目录的基本操作,供读者参阅。

图12-4-3 基本操作

然后点击菜单 View ,并选中如图12-4-4所示的项目,即可在代码块中显示行号。

图12-4-4 显示代码块中的行号

将鼠标移动到代码块中并单击,如图12-4-5所示,开始输入一行代码,然后回车,输入第二行——注意,这里与 Python 交互模式不同,回车意味着换行,而不是执行当前行代码。输入完图12-4-5所示的所有代码。

图12-4-5 编写代码

按照图12-4-5所示,点击该按钮,可以执行当前代码块;或者按组合键 shift + Return/Enter 执行,效果如图12-4-6所示。

图12-4-6 执行代码块

执行了代码块 [1] 之后,执行结果显示在代码之后,并自动进入下一个代码块。在编辑界面顶部,有各种常见的编辑、控制按钮,把鼠标滑动到按钮上,会显示该按钮的作用。

以上是 JupyterLab 的最基本使用方法,其他复杂功能,可以在应用过程中学习。

12.4.2 第三方库

Python 生态中拥有非常丰富的支持科学计算的第三方库,常用的有 NumPy 、Pandas 、SciPy 、Matplotlib 、SymPy 等,建议读者将这些库依次安装。

代码语言:javascript
复制
% pip install numpy pandas sympy scipy matplotlib

除了能在终端执行安装指令之外,在 JupyterLab 的代码块中也可以执行终端指令,如图12-4-7所示,在代码块中输入 !pip install numpy 并执行——注意前面的 ! 符号。其效果等同于以往在终端执行安装 pip install numpy 指令。

图12-4-7 在代码块中执行安装指令

安装好之后,在代码块中输入如下代码,并执行,即可查看所安装的 NumPy 的版本。

代码语言:javascript
复制
[3]: import numpy as np
     np.__version__
[3]: '1.19.4'       # 上述代码块的输出结果

在数据科学中,引入一些常用的第三方库时,习惯于再命名一个别称(或简称),例如以 np 作为 numpy 的别称,并且这是一种习惯命名,如果非要以 ny 为别称,语法上没有问题,但不是大多数人的习惯——代码更多时候是给人看的。

安装好基础库之后,再列举几个示例(随后几个小节内容),体会 Python 在科学计算中的应用。

12.4.3 矩阵

矩阵不仅在线性代数中占据重要地位,也是科学计算的主角。如果读者还忌惮于当初用纸笔完成有关矩阵计算的痛苦,现在使用 Python 中科学计算的工具包则会体验到无比的畅快。

在 JupyterLab 代码块中输入如下代码(如无特别声明,本节的代码均在 JupyterLab 中调试)。

代码语言:javascript
复制
[4]: import numpy as np
     I_m = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
     I_m
[4]: array([[1, 0, 0],
            [0, 1, 0],
            [0, 0, 1]])

变量 I_m 引用的对象是用 np.array() 创建的二维数组,可以用二维数组表示矩阵—— I_m 表示一个

3\times3

的单位矩阵。

此外,用 NumPy 中的 np.mat() 也能生成矩阵对象:

代码语言:javascript
复制
[5]: np.mat("1 2 3; 4 5 6")
[5]: matrix([[1, 2, 3],
             [4, 5, 6]])

以这两种形式表示的矩阵,在矩阵加法和数量乘法上没有任何差别,例如:

代码语言:javascript
复制
# 二维数组表示的矩阵相加
[6]: a = np.array([[1, 4], [8, 9]])
     b = np.array([[2, 3], [7, 11]])
     a + b
[6]: array([[ 3,  7],
            [15, 20]])
代码语言:javascript
复制
# 矩阵对象相加
[7]: ma = np.mat("1 4; 8 9")
     mb = np.mat("2 3; 7 11")
     ma + mb
[7]: matrix([[ 3,  7],
             [15, 20]])
代码语言:javascript
复制
# 矩阵数量乘法
[8]: print('array: ', 7 * a)
     print("")
     print('mat: ', 7 * ma)
[8]: array:  [[ 7 28]
              [56 63]]

     mat:  [[ 7 28]
            [56 63]]

但是,矩阵乘法会有所不同。首先来看数组 ab ,如果计算 a * b ,其所得并不是它们所表征的矩阵乘法结果。

代码语言:javascript
复制
[9]: a * b    # 数组相乘,不是矩阵乘法
[9]: array([[  2,  12],
            [ 56,  99]])

而矩阵对象 mamb 直接相乘 ma * mb 的结果是根据矩阵相乘规则所得。

代码语言:javascript
复制
[10]: ma * mb
[10]: matrix([[ 30,  47],
              [ 79, 123]])

如果用数组表示矩阵,且计算矩阵乘法,可以通过 np.dot() 函数或者数组的方法 dot() 实现。

代码语言:javascript
复制
[11]: np.dot(a, b)   # 或者 a.dot(b)
[11]: array([[ 30,  47],
             [ 79, 123]])

在实际的项目中,经常会遇到稀疏矩阵——详细内容请参阅拙作《机器学习数学基础》(电子工业出版社)。例如:

代码语言:javascript
复制
[12]: sm = np.array([[1, 0, 2, 0, 0, 0, 0, 0],
                     [0, 0, 0, 0, 0, 0, 0, 0],
                     [0, 0, 3, 0, 0, 0, 0, 0]])
      sm
[12]: array([[1, 0, 2, 0, 0, 0, 0, 0],
             [0, 0, 0, 0, 0, 0, 0, 0],
             [0, 0, 3, 0, 0, 0, 0, 0]])

sm 表示一个稀疏矩阵,由于行列数还不是很多,像上面那样写出来还能忍受。稀疏矩阵包含了很多零,虽然耗费了时间和精力,但它们事实上没有什么意义,还占用内存空间。为此可以对此类矩阵进行压缩,并直接生成压缩矩阵对象。

代码语言:javascript
复制
[13]: from scipy.sparse import csr_matrix
      row = np.array([0, 0, 2])
      col = np.array([0, 2, 2])
      data = np.array([1, 2, 3])
      csrm = csr_matrix((data, (row, col)), shape=(3, 8))
      csrm
[13]: <3x8 sparse matrix of type '<class 'numpy.longlong'>'
             with 3 stored elements in Compressed Sparse Row format>

现在得到的压缩矩阵 csrm 如果转化为数组,与前面创建的 sm 一样。

代码语言:javascript
复制
[14]: csrm.toarray()
[14]: array([[1, 0, 2, 0, 0, 0, 0, 0],
             [0, 0, 0, 0, 0, 0, 0, 0],
             [0, 0, 3, 0, 0, 0, 0, 0]], dtype=int64)

本书作者在《机器学习数学基础》(http://math.itdiffer.com)一书中,专门介绍了矩阵及其有关计算的 Python 实现方法,推荐有兴趣的读者深入学习参考。

12.4.4 解线性方程组

最一般的解线性方程组的方法是高斯消元法,在传统的数学教材中,还会列出其他巧妙的方法。这个工作如果教给 Python ,就是找到适合的函数——背后的事情教给了发明者,如果读者立志创新研究,可以自己成为发明者。

例如,求方程组的解:

\begin{cases}-x_1 + 3x_2 - 5x_3 &= -3 \\2x_1 -2x_2 + 4x_3 &= 8 \\ x_1 + 3x_2 &= 6\end{cases}

用矩阵的方式,可以将方程组表示为:

\begin{bmatrix}-1&3&-5\\2&-2&4\\1&3&0\end{bmatrix}\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}=\begin{bmatrix}-3\\8\\6\end{bmatrix}

然后,编写如下代码:

代码语言:javascript
复制
[15]: A = np.mat("-1 3 -5; 2 -2 4; 1 3 0")    # 系数矩阵
      b = np.mat('-3 8 6').T                  # 常数项
      r = np.linalg.solve(A, b)               # 求解
      print(r)
[15]: [[ 4.5]
       [ 0.5]
       [-0. ]]

这里使用 np.linalg.solve() 函数求解线性方程组,从输出结果可知,其解为

x_1 = 4.5, x_2 = 0.5, x_3 =0

但是,如果遇到这样的方程组:

\begin{cases}x_1+3x_2-4x_3+2x_4&=0\\3x_1-x_2+2x_3-x_4&=0\\-2x_1+4x_2-x_3+3x_4&=0\\3x_1+9x_2-7x_3+6x_4&=0\end{cases}

还是使用前面的函数,对此方程组求解。

代码语言:javascript
复制
[17]: A = np.mat("1 3 -4 2;3 -1 2 -1;-2 4 -1 3;3 0 -7 6")
      b = np.mat("0 0 0 0").T
      r = np.linalg.solve(A, b)
      print(r)
[17]: [[ 0.]
       [ 0.]
       [-0.]
       [ 0.]]

当所有变量的解都是 0 ,原线性方程组成立,但这仅仅是一个特解。根据线性代数的知识可以判断,此方程组有无穷多个解(参阅《机器学习数学基础》2.4.2节),还能用程序计算吗?

代码语言:javascript
复制
[18]: from sympy import *
      from sympy.solvers.solveset import linsolve
      x1, x2, x3, x4 = symbols("x1 x2 x3 x4")
      linsolve([x1 + 3*x2 - 4*x3 + 2*x4, 
                3*x1 - x2 + 2*x3 - x4, 
                -2*x1 + 4*x2 - x3 + 3*x4, 
                3*x1 +9*x2 - 7*x3 + 6*x4], 
               (x1, x2, x3, x4))

输出:

这就是该线性方程组的通解,如果对应到未知量,即:

\begin{cases}x_1=\frac{x_4}{10}\\x_2=-\frac{7}{10}x_4\\x_3=0\\x_4=x_4\end{cases}

其中

x_4

是自由变量。

12.4.5 假设检验

假设检验是统计学中的重要内容,也广泛地应用在科学研究和工程项目中。下面选用《机器学习数学基础》6.5.1节的一个案例,读者从中感受 Python 在这方面的应用。

假设有人制造了一个骰子,他声称是均匀的,也就是假设分布律为:

H_0:P(X=i)=\frac{1}{6}, (i=1,2,3,4,5,6)

为了证明自己的判断,他做了

n=6\times10^{10}

次投掷试验,并将各个点数出现的次数记录下来(为了便于观察,出现次数用

10^{10}

加或减一个数表示。

点数

1

2

3

4

5

6

次数

接下来利用

\chi^2

检验法,通过如下程序对此人的结论——假设——进行检验。

代码语言:javascript
复制
[19]: from scipy.stats import chisquare
      chisquare([1e10-1e6, 1e10+1.5e6, 
                 1e10-2e6, 1e10+4e6, 1e10-3e6, 1e10+0.5e6])
[19]: Power_divergenceResult(statistic=3250.0, pvalue=0.0)

输出结果显示,检验统计量的值

\chi^2 =3250.0

。根据如下公式(对此公式详细讲解,请参阅《机器学习数学基础》6.5.1节):

\chi^2 = \sum_{i=1}^k\frac{(np_i-f_i)^2}{np_i}

此处

𝑘=6

,则在

\alpha = 0.05

的显著水平下,

\chi^2_{0.05}(6-1)

的值为:

代码语言:javascript
复制
[20]: from scipy.stats import chi2 
      chi2.isf(0.05, (6-1))
[20]: 11.070497693516355

由输出值可知,在显著水平

\alpha = 0.05

下,

\chi^2\gt\chi^2_{0.05}(6-1)=11.07

,故试验数据不支持“骰子均匀”这个假设。

还可以根据 p 值检验此假设。

代码语言:javascript
复制
[21]: p_value = 1 - chi2.cdf(3250.0, (6-1)) 
      print(p_value)
[21]: 0.0

得到的 p 值结果说明拒绝原假设犯错误的概率是 0.0% 。

自学建议 科学计算不仅仅是科研和工程项目的必备技能,也是后面所列举的数据分析、机器学习的基础。所以,有意在数据科学方向发展的读者务必要学习科学计算。首先,要具备相关的数学基础,其次要掌握相关计算的第三方包的引用,最后要将前两者应用到实际问题之中。也正是基于这些思考,我出版了《机器学习数学基础》(网址:http://math.itdiffer.com),在这本书中,不强调传统数学教材中的“纸笔计算”,重点是在理解有关数学原理之后,用程序工具完成计算,并以贴近真实的问题为案例。有兴趣的读者不妨参阅。”

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2022-08-17,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 老齐教室 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 12.4 科学计算
    • 12.4.1 Jupyter
      • 12.4.2 第三方库
        • 12.4.3 矩阵
          • 12.4.4 解线性方程组
            • 12.4.5 假设检验
            领券
            问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档