高等数学整理(二)

原创
2021/04/06 11:33
阅读数 5.4K

高等数学整理

多元复合函数求导法则

多元复合函数是用在bp神经网络或者叫做神经网络的bp算法当中。深度学习是基于深度神经网络的。多元复合函数在神经网络算法当中有很大的用处。习惯性当中,把多元复合函数求导法则称为链式法则。

多元函数全增量

全增量为
△z=f(x0+△x,y0+△y)-f(x0,y0)

=f(x0+△x,y0+△y)-f(x0+△x,y0)
+f(x0+△x,y0)-f(x0,y0)

=fy'(x0+△x,ξ1)·△y+fx'(ξ2,y0)·△x
此处用到的为拉格朗日中值定理,其中ξ1在y0与y0+△y之间,
ξ2在x0与x0+△x之间】

然后利用一阶偏导数的连续性,
【ξ1→y0,ξ2→x0】

fy'(x0+△x,ξ1)=fy'(x0,y0)+α
fx'(ξ2,y0)=fx'(x0,y0)+β
【其中,α,β是无穷小】

从而,
△z=fy'(x0,y0)△x+fx'(x0,y0)△y
+α△x+β△y

其中,α△x+β△y=o(ρ)

全增量等于两个偏增量的和

  • 一元函数与多元函数复合

若函数u=å(t)及v=ß(t)都在t点可导,函数z=f(u,v)在对应点(u,v)具有连续偏导数,那么复合函数z=f[å(t),ß(t)]在点t可导,且

证明:设t获得增量,此时u=å(t)、v=ß(t)的对应增量为

函数的全增量

->0,->0时,->0,->0

等式两边同时除以

  • 多元函数与多元函数复合

如果函数u=å(x,y),v=ß(x,y)都在点(x,y)具有对x及对y的偏导数,函数z=f(u,v)在对应点(u,v)具有连续偏导数,那么复合函数z=f[å(x,y),ß(x,y)]在点(x,y)的两个偏导数都存在,且有

证明的过程跟一元函数跟多元函数复合一样,只不过求对x的偏导的时候把y看成常数,求对y的偏导的时候把x看成常数。

类似的,设u=å(x,y),v=ß(x,y)及w=(x,y)都在点(x,y)具有对x及对y的偏导数,函数z=f(u,v,w)在对应点(u,v,w)具有连续偏导数,则复合函数z=f[å(x,y),ß(x,y),(x,y)]在点(x,y)的两个偏导数都存在,且可用下列公式计算:

如果函数u=å(x,y)在点(x,y)具有对x及对y的偏导数,函数v=ß(y)在点y可导,函数z=f(u,v)在对应点(u,v)具有连续偏导数,那么复合函数z=f[å(x,y),ß(y)]在点(x,y)的两个偏导数都存在,且有

设z=f(u,x,y)具有连续偏导数,u=å(x,y)具有偏导数,则复合函数z=f[å(x,y),x,y]具有对自变量x及y的偏导数,且有

多元复合函数求导法则示例

  • ,求

令u=x+y,v=xy,则

隐函数求导公式

隐函数:一般地,如果变量x和y满足一个方程g(x,y)=0,在一定条件下,当x取某区间内的任一值时,相应地总有满足这方程的唯一的y值存在,那么就说方程g(x,y)=0在该区间内确定了一个隐函数。

例如

这里y=3√(2-x),所以当x取任意值时,y的值唯一,所以这里确定了一个隐函数。大致图形如下

,这个方程能否确定一个隐函数,它的图形大致如下

由图形可知,在(0,1)范围内随意取一个x值,y的值不唯一,所以这个方程不能确定一个隐函数。

在点(0,1)的某个邻域内存在隐函数吗?

我们上面一题中其实是从一元函数的角度来看待的,而这里其实是从二元函数的角度来看待的,(0,1)这个点就是A点,A点点邻域内,它是可以存在隐函数的,只要这个邻域不过大。

隐函数求导公式:设函数g(x,y)在点P(x0,y0)的某一邻域内具有连续偏导数,且g(x,y)=0,,则方程g(x,y)=0在点(x0,y0)的某一邻域内恒能唯一确定一个连续导数的函数y=f(x),他满足条件y0=f(x0),并有,同时也被称为隐函数存在定理1。

证明:将y=f(x)代入g(x,y)=0得恒等式g(x,f(x))0,左端看作是x的一个复合函数,两端求导得,将其变型后就可得

设函数g(x,y,z)在点的某一邻域内具有连续偏导数,且,则方程g(x,y,z)=0在点的某个邻域内恒能唯一确定一个连续且有连续偏导数的函数z=f(x,y),它满足条件,并有,同时也被称为隐函数存在定理2。

隐函数求导公式示例

验证方程在点(0,1)的某一个邻域内能唯一确定一个有连续导数,当x=0,y=1时的隐函数y=f(x),并求这函数的一阶导数在x=0的值。

,则

根据隐函数存在定理1,可以确定隐函数y=f(x)存在,且一阶导数x=0的值为0

多元函数极值

  • 二元函数极值
  1. 函数z=f(x,y)定义域为D,为D内点,若存在的某个邻域,使得对于该邻域内异于的任何点(x,y)都有,则称为函数的极大值。极大值可视为局部最大值。
  2. 函数z=f(x,y)定义域为D,为D内点,若存在的某个邻域,使得对于该邻域内异于的任何点(x,y)都有,则称为函数的极小值。极小值可视为局部最小值。

极大值与极小值均称为极值。

极值点示例

  • 函数在(0,0)处取得极小值

函数图像如下

其中绿色坐标为x轴,红色坐标为y轴,蓝色坐标为z轴。根据图像可以看到,(0,0)处取得极小值,同时也是最小值。

  • 函数在(0,0)处取得极大值

函数图像如下

根据图像可以看出(0,0)处取得极大值,同时也是最大值。

  • 函数z=xy在(0,0)处不取得极值

函数图像如下

这是双曲抛物面,形状类似于马鞍。我们可以看到在(0,0)处取不到极值,因为在(0,0)的邻域内,有比它大的函数值也有比它小的函数值。

极值必要条件:设函数z=f(x,y)在点具有偏导数,且在该点取得极值,则有

证明:设函数在点处取得极大值,则在的某个邻域内异于的点都满足,故,既是一元函数处取得极大值,根据一元函数极值的必要条件,则,同理

极值必要条件推广:若三元函数f(x,y,z)在点取得极值,则有,可推广至更多元函数。

极值充分条件:设函数z=f(x,y)在点的某邻域内连续,且一阶二阶偏导数连续,又满足,记

  1. 时有极值,且当A<0时有极大值,A>0时有极小值
  2. 时没有极值
  3. 时,可能存在极值,也可能不存在极值

求函数极值示例

  • 求函数极值

解方程组(一阶偏导数)

这是每一个偏导数方程都是一个一元二次方程,对x的偏导方程中,根据一元二次方程求根公式,将a=3、b=2、c=-2代入得x1=(√7-1)/3,x2=-(√7+1)/3。对y的偏导方程中,将a=-3、b=2、c=0代入得y1=0,y2=2/3。则得到4组解((√7-1)/3,0),((√7-1)/3,2/3),(-(√7+1)/3,0),(-(√7+1)/3,2/3)。

二阶偏导数

  1. 对于第一组解((√7-1)/3,0),A=2√7,B=0,C=2,当√7取正时,则AC=4√7>=0有极值,且A>0有极小值;当√7取负时,则AC=4√7<=0,没有极值
  2. 对于第二组解((√7-1)/3,2/3),A=2√7,B=0,C=-2,当√7取正时,则AC=-4√7<=0,没有极值;当√7取负时,则AC=-4√7>=0有极值,且A<0有极大值
  3. 对于第三组解(-(√7+1)/3,0),A=4-2√7,B=0,C=2,当√7取正时,则AC=8-4√7<=0,没有极值;当√7取负时,则AC=8-4√7>=0有极值,且A>0有极小值
  4. 对于第四组解(-(√7+1)/3,2/3),A=4-2√7,B=0,C=-2,当√7取正时,则AC=4√7-8>=0有极值,且A<0有极大值;当√7取负时,则AC=4√7-8<=0,没有极值

条件极值

  • 一个长方体表面积为,求该长方体体积的最大值

设长宽高分别为x、y、z,在条件下,求V=xyz最大值。

这道题的解题思路可以通过条件式将x写成y、z的表示形式,代入V=xyz中,再求一阶偏导数方程以及二阶偏导数。通过解来看是否有极值,再求极值。也就是采用消元法,转化成二元函数最大值。

但有的时候,有些约束条件,并不好将一个未知数表示成其他的未知数的形式,比如,这又该如何处理呢?

拉格朗日乘数法:考虑函数z=f(x,y)在条件g(x,y)=0下的极限。作辅助函数,将其视为三元函数,根据极值存在必要条件

我们继续求长方体体积极值问题。

作辅助函数

由上面的求解过程,我们可以看到其实求最左边的方程组其实并不是那么容易的。

拉格朗日乘数法的局限

  1. 方程变量过多时,即使是简单线性方程组,求解困难。
  2. 有复杂方程存在时很难求解。
  3. 这只是必要条件,而非充要条件,所以即便我们正确求解出来之后,我们无法判断求出来的是极大值还是极小值。
  4. 多变量情况下,需根据二阶偏导数构成的海森矩阵的正定性判断是否属于极值点。

全微分

之前我们在讲一元函数的微分的时候,指的是函数变化量的线性部分∆y=A•∆x+o(∆x)

设函数z=f(x,y)在点(x,y)的某邻域有定义,如果函数在点(x,y)的全增量∆z=f(x+∆x,y+∆y)-f(x,y)可表示为,A、B与∆x、∆y无关,称函数在点(x,y)可微分,A•∆x + B•∆y称函数在点(x,y)处的全微分,记为dz=A•∆x + B•∆y

全微分必要条件

如果函数z=f(x,y)在点(x,y)可微分,则该点处偏导数必存在,且函数在点(x,y)的全微分为。因为这里是必要条件而不是充分条件,所以偏导数存在,不一定可微分

全微分充分条件

如果函数z=f(x,y)的偏导数在点(x,y)连续,则在该点处可微分

理解全微分

之前在一元函数中说过如果在某点微分存在,在该点附近可以用该点切线近似表示这个函数。

在二元函数中,如果在某点微分存在,在该点附近可以用该点的切平面近似表示这个函数。

全微分计算示例

  • 计算函数的全微分

首先计算偏导数

全微分为

全微分的用处

  • 近似计算

根据f(x+∆x,y+∆y)-f(x,y),则f(x+∆x,y+∆y)=f(x,y)+,由于是一个无穷小量,我们可以将其舍弃,则=1+3*0.01+0*0.02=1.03(将x=1、y=3代入,∆x=0.01,∆y=0.02),则≈1.03

连续、偏导数、可微分三者之间的关系

  • 一元函数的情况

  • 二元函数的情况

二元函数相对一元函数要复杂的多,它对我们在理解上的困难是呈几何倍数增加的。但是在理解了二元函数的基础上再去理解多元函数(三元或者三元以上)就没有那么困难了。现在我们来看一下一些我们已知的情况。

  • 现在我们来看一下连续能否推导出偏导数存在吗?答案是不能。以一个反例来说明

根据偏导数的定义,在(0,0)点对x的偏导数为

当∆x趋近于0+的时候,=1,而∆x趋近于0-的时候,=-1;也就是说这个极限不存在。所以虽然函数连续,但在(0,0)处偏导数不存在。在(0,0)处对y的偏导数也一样。

  • 偏导数存在能推出连续吗?答案是不能,以一个反例来说明

根据偏导数的定义,在(0,0)点对x的偏导数为=0,同理,而我们前面也说过,该二元函数的其他方向上的趋近并不为0,故虽然偏导数存在,但并不连续,既然不连续则不可微。

则最终三者之间的关系如下图

方向导数与梯度下降算法

方向导数

偏导数反映函数沿坐标轴方向的变化率,如何表示任意方向的变化率,我们将这个任意方向的变化率称为方向导数。

方向导数的意义:可用于搜索函数极值。

假设有一个人位于山坡的A点,现在他要迅速到达山顶的位置,那么他就要选择变化率最大的方向搜索。他可以向四周迈出一小步,当他迈出之后,他可以比较迈出之后的位置的海拔高度差,他就选择那个海拔高度差最大的位置的方向走,那么就可以较快地到达山顶。这同样也是一种搜索函数极值的方式

物理学、气象学应用,温度、气压在某些方向的变化率。

射线的参数方程

单位方向向量,射线以为起点,且于同向,求射线的参数方程。

设射线上的点为P(x,y),则

,即=(tcos å,tcos ß),t≥0

转化为

由于为单位向量,所以它的模为1,即

函数在某方向的增量

函数z=f(x,y)在某邻域有定义,到P函数增量

到P的距离为:

方向导数:比值,P沿l趋近时,比值极限即为方向导数

理解方向导数

  1. 是关于方向的函数
  2. 通过表示方向,是的函数,记为
  3. 的最大值表示原函数在某个方向上的具有最大变化率。

方向导数化简

根据全微分,A为函数对x的偏导数,B为函数对y的偏导数,则有

这里∆x为,∆y为

代入化简

因为o(t)是t的高阶无穷小,根据高阶无穷小的定义,b是比a高阶的无穷小,则o(t)/t=0

方向导数的最值

记向量,则。其中为两个向量的夹角,表示两个向量的点乘,又等于两个向量的模相乘再乘以两个向量的夹角余弦。具体可以参考线性代数整理

要想让取得最大值,由于=1,=,只有=1的时候,才能取得最大值。故=的时候,即两个向量同向时,取得最大值

又只有=-1的时候,才能取得最小值。故=180°的时候,即两个向量反向时,取得最小值

梯度:向量称为梯度,记为,梯度指示了函数变化率最大的方向。梯度反方向表示函数减小最快的方向。所以可利用梯度搜索函数极值。

如何实现按方向移动

初始位置,按照向量指示的方向移动,移动后的位置为,则,这里因为同向,既是

梯度下降算法

现在我们回到多元函数,这里是一个n元函数

梯度:表示各个变量在某个点的偏导数所构成的向量,这个向量表示移动的方向在函数值的变化上是否最大化。

极大值迭代:这里的i代表的是变量,如i=1代表第一个变量,i≤n,最多为n个变量。k代表迭代的次数。m称为迭代步长,又称为学习率

极小值迭代:极小值迭代又叫梯度下降算法。

梯度下降算法示例

  • 求函数最小值

这个函数比较简单,最小值很明显就是x1=0、x2=0的情况下,函数值最小值为0,现在我们来看看用梯度下降法怎么算

  1. 求偏导数:Lx1(x1,x2)=2•x1                       Lx2(x1,x2)=4•x2
  2. 赋初值:,这里我们把迭代步长设置为0.01
  3. 迭代:第一次迭代,根据,则有,将x1=0.98、x2=2.88代入原函数得

通过第一次迭代,我们发现函数值是下降的。那么后续的迭代也是一样的,这里就不继续算下去了。

整个算法流程如下

代码实现

class GradientDescent:

    def __init__(self, m, x1, x2):
        # 利用梯度下降算法,求函数最小值:L = x1^2 + 2*x2^2
        self._m = m         #迭代步长,学习率
        self._x1 = x1
        self._x2 = x2

    def fun_min(self):
        L = self._x1**2 + 2 * self._x2**2
        # 迭代次数
        k = 0
        # 是否继续循环
        err = 1
        # 误差值
        threshold = 0.000000001
        while err > threshold:
            # 2*x1为函数对x1的偏导数
            self._x1 = self._x1 - self._m * 2 * self._x1
            # 4*x2为函数对x2的偏导数
            self._x2 = self._x2 - self._m * 4 * self._x2
            # 计算前后两次迭代后的函数差值的绝对值
            err = abs(self._x1**2 + 2 * self._x2**2 - L)
            L = self._x1 ** 2 + 2 * self._x2 ** 2
            # 更新迭代次数
            k += 1
        return (self._x1, self._x2, L, k)
from playLA.GradientDescent import GradientDescent

if __name__ == "__main__":
    grad = GradientDescent(0.01, 1, 3)
    print(grad.fun_min())

运行结果

(0.00015563843334256782, 6.065078047552769e-08, 2.4223329290363272e-08, 434)

结果第一个数为第一个变量x1经过434此迭代后的结果,第二个数是第二个变量x2经过434次迭代后的结果,第三个数是经过434次迭代后的函数值,第四个数为迭代次数。

现在我们来把迭代前后差值的绝对值用图形表示出来

import matplotlib.pyplot as plt

class GradientDescent:

    def __init__(self, m, x1, x2):
        # 利用梯度下降算法,求函数最小值:L = x1^2 + 2*x2^2
        self._m = m         #迭代步长,学习率
        self._x1 = x1
        self._x2 = x2

    def fun_min(self):
        L = self._x1**2 + 2 * self._x2**2
        # 迭代次数
        k = 0
        # 是否继续循环
        err = 1
        # 误差值
        threshold = 0.000000001
        value = []
        while err > threshold:
            # 2*x1为函数对x1的偏导数
            self._x1 = self._x1 - self._m * 2 * self._x1
            # 4*x2为函数对x2的偏导数
            self._x2 = self._x2 - self._m * 4 * self._x2
            # 计算前后两次迭代后的函数差值的绝对值
            err = abs(self._x1**2 + 2 * self._x2**2 - L)
            value.append(err)
            L = self._x1 ** 2 + 2 * self._x2 ** 2
            # 更新迭代次数
            k += 1
        plt.plot(value)
        plt.show()
        return (self._x1, self._x2, L, k)
from playLA.GradientDescent import GradientDescent

if __name__ == "__main__":
    grad = GradientDescent(0.01, 1, 3)
    print(grad.fun_min())

运行结果

(0.00015563843334256782, 6.065078047552769e-08, 2.4223329290363272e-08, 434)

如果我们修改了函数表达式或者是参数之后,有可能出现不收敛的情况,导致我们的程序一直运行。为了避免这种情况发生,我们可以加一些异常处理机制,比如将迭代次数限制在一个范围之内

import matplotlib.pyplot as plt

class GradientDescent:

    def __init__(self, m, x1, x2):
        # 利用梯度下降算法,求函数最小值:L = x1^2 + 2*x2^2
        self._m = m         #迭代步长,学习率
        self._x1 = x1
        self._x2 = x2

    def fun_min(self):
        L = self._x1**2 + 2 * self._x2**2
        # 迭代次数
        k = 0
        # 是否继续循环
        err = 1
        # 误差值
        threshold = 0.000000001
        value = []
        while err > threshold and k < 10000:
            # 2*x1为函数对x1的偏导数
            self._x1 = self._x1 - self._m * 2 * self._x1
            # 4*x2为函数对x2的偏导数
            self._x2 = self._x2 - self._m * 4 * self._x2
            # 计算前后两次迭代后的函数差值的绝对值
            err = abs(self._x1**2 + 2 * self._x2**2 - L)
            value.append(err)
            L = self._x1 ** 2 + 2 * self._x2 ** 2
            # 更新迭代次数
            k += 1
        plt.plot(value)
        plt.show()
        return (self._x1, self._x2, L, k)

梯度下降算法的局限性

  1. 初值敏感,变量越多越敏感,不同的初值可能收敛到不同的极值点
  2. 学习率敏感,可动态改变学习率。(学习率特别大的话,会导致我们迭代过程不收敛;学习率特别小的话,会导致最后效果不是特别好,并且迭代次数会增加)
  3. 无法完全保证取得的是全局最值。

偏导数与方向导数之间的关系

  • 求函数在点(1,2)处沿从点(1,2)到点(3,5)的方向的方向导数

首先,方向向量为(3,5)-(1,2)=(2,3),方向余弦为,这里方向余弦的分母即为√(2^2+3^2)=√13,分子分别为2、3。而也就为单位方向向量。求偏导数,最后代入方向导数公式,得,这里x=1,y=2。

现在我们想知道的是方向导数存在能推出偏导数存在吗?

答案是不可以,上面的例题中,之所以能通过公式求出方向导数是在函数可微分的前提下推导出来的。现在我们来举出一个反例

这个函数我们之前在连续能否推出偏导存在的时候已经说明过,它的偏导数是不存在的。

根据方向导数的定义,设单位方向向量为

说明函数在(0,0)处方向导数存在,但偏导数不存在。

同样偏导数存在,不一定存在方向导数。因为对于多元函数来说,它是0~360度逼近的,任何一个方向的导数存在都不能说明其他方向的导数存在,而偏导数其实也就是其中的两个方向而已。

方向导数存在也不能推出连续,而连续也不能推出方向导数存在。同时方向导数存在不能推出可微分,但是可微分是一定可以推出方向导数存在的。

二元函数的泰勒公式

z=f(x,y)在点(x0,y0)的某一邻域内连续且有(n+1)阶连续偏导数,(x0+h,y0+k)为此邻域内任一点,则有

对比一下一元函数的泰勒展式,我们可以看到它们在写法上有一些不同。

以一阶导数项为例:一元函数f'(x0)(x-x0)其实就是f'(x0)∆x,而二元函数的h和k就相当于∆x和∆y,那么,等于就是一元函数的一阶导数乘以增量∆x变成了二元函数的两个偏导数分别乘以x,y的增量∆x和∆y再相加。

而二阶导数项一元函数为即为∆x^2f''(x0)/2!,二元函数撇开1/2!,则,,它比一元函数多了一个,即为对x求偏导后再对y求偏导项,其他两项跟一元函数类似。

被称为拉格朗日余项。这里不做证明。对于二元函数的泰勒公式,我们有一个了解就好了。

积分定律

导函数与原函数

在某个定义区间上,如果F'(x)=f(x)或dF(x)=f(x)dx,则F(x)称为f(x)的一个原函数。

比如sin(x)是cos(x)的一个原函数,因为sin'(x)=cos(x)。sin(x)+C(C是常数)也是cos(x)的一个原函数,因为(sin x + C)'=cos x。

连续函数一定有原函数

不定积分:函数f(x)带有任意常数项的原函数称为f(x)的不定积分。它的表示为,积分号由莱布尼茨发明,Sum首字母拉长,它代表求和的意思。

  • 求不定积分,求1/x的不定积分

因为ln'x=1/x,根据lnx的图像我们看到

所以x的定义域为大于0的,则

当x>0时,

当x<0时,,则

综上,

常见积分表

积分的性质

求和公式

数量乘公式

第一类换元积分

(函数的两端求微分)

则有

我们来看一个例子

令u=1+2x

这里1+2x对应上面的就是,常数2可以提到积分的前面,而就是,因为我们要构造的导数即(1+2x)',因为(1+2x)'=2,所以要乘以一个1/2。这是的由来。整体构成u的微分,则有,最后根据积分公式,求得

第二类换元积分

         (表示t为的反函数)

条件:单调可导(只有单调的,才能求反函数),具有原函数

例子

(目的是为了消除根号),(在这个范围之内是单调的,存在反函数),(t是的反函数),(求微分)

这里第一步将代入可得;第二步,将4a^2提出,根据余弦的平方公式,则有;第三步,根据积分的求和公式,则有;第四步,根据可得+C,根据,又中t是自变量,则2dt=d2t,可得=;第五步,我们要将t换回x,将代入,则有=a^2sin(2arcsin(x/2a)),根据sin2x=2sinxcosx可得a^2sin(2arcsin(x/2a))=2a^2sin(arcsin(x/2a))cos(arcsin(x/2a)),根据sin(arcsinx)=x以及cos(arcsinx)=√(1-x^2),可得2a^2sin(arcsin(x/2a))cos(arcsin(x/2a))=2a^2(x/2a)√(1-(x/2a)^2)=ax√[(4a^2-x^2)/4a^2]=x/2√(4a^2-x^2),则=

分部积分

之前已知,则

两边求不定积分,得

这里u、v都是关于x的函数

例子

令u=x,dv=sinxdx,则du=dx,v=-cosx

u与dv的选取原则

  1. v易求
  2. ∫vdu比∫udv更容易求

令u=x,,则du=dx,

不定积分是微分的逆运算。是否所有函数的原函数都有解析式?答案是并不是所有函数的原函数都有解析式。

定积分的定义

我们来看一个面积问题

  • 与x轴围成的封闭图形的面积

微元法:-1到1平均分成n段,每段长度(-1到1的总长度为2)

第i个矩形的高(右侧)为,这里其实就是找出第i个矩形右侧的横坐标x=-1+(2/n)*i(起始点+总距离)

封闭面积

根据平方和公式以及1+2+3+...+n=n(n+1)/2,将展开后,再代入这两个公式,最终可得,当n越大,越接近真实面积。现在我们来看一下n在几个不同的值的情况下的图形

n=6

n=32

n=88

由图中我们可以看到,当n在逐渐的增大,矩形在逐渐的减小。

现在我们用的是矩形的右侧高来计算面积的,那么我们使用左侧高来计算会如何呢,将微元部分使用梯形来计算结果是否一致。微元面积是否可以任意划分。

现在我们将n=20来看一下各种情况的封闭面积,首先是右侧高

此时计算出来的封闭面积为1.33

左侧高

此时算出来的封闭面积也是1.33

梯形微元

此时算出来的封闭面积也是1.33。但这里要说明的是这个函数比较特殊,一般大部分函数在使用微元法计算面积的时候只有当n趋近于∞的时候,各种微元法的最终结果才相同,n越小,差距越大。

之前我们都是等分矩形,矩形的底边都是相等的,若任意划分矩形,矩形的底边不一定相等。区间固定,小矩形最大的底边趋于0,面积是否为定值。

数学描述:函数f(x)在[a,b]上有界,在[a,b]上任意插入若干分点,把区间[a,b]分成了n个小区间,每个小区间的长度,记

定积分:每个小区间上任取一点,若当时,极限总存在,且与闭区间[a,b]的分法及点的取法无关,则称该极限为函数f(x)在区间[a,b]上的定积分(简称积分),记作,即

定积分的意义在于面积,它的几何含义如下

这里需要注意的是,第一个函数的定积分指的是阴影部分的面积,第二个函数的定积分指的是x轴上方的阴影部分面积之和减去x轴下方的阴影面积。

定积分的物理意义

路程        ,有一个物体在一条直线上做变速运动,那么它在t1到t2这段时间内走过的路程就等于v(t)在t1到t2这段时间内的定积分。

做功           ,有一个变力推动一个物体直线运动,力与运动方向一致,那么这个变力在s1到s2这段距离内做的功就等于F(s)在s1到s2这段距离内的定积分。

定积分存在的充分条件

  1. 函数在闭区间上连续
  2. 函数在闭区间上有界,且只有有限个间断点

性质

  1. 为常数
  2. 若在[a,b]上f(x)≥0,则
  3. 若在[a,b]上f(x)≤0,则
  4. 若在[a,b]上f(x)最大值为M,最小值为m,则

定积分中值定理:若f(x)在[a,b]上连续,则在[a,b]上至少存在一个点,满足

几何解释

定积分的面积=矩形面积

可以找到一个矩形,它的面积等于定积分的面积(阴影部分的面积),主要是这个矩形的高是一定存在的

物理解释:总有一个时刻的瞬时速率等于平均速率,这里v(å)为平均速率,表示总路程除以总时间,这也是区间测速原理。

牛顿-莱布尼茨公式

计算定积分

[0,1]平均分成n个小区间,分点为,取,和式为

这里是将代入得i^3/n^4(i=1,2,...,n)在根据1到n的立方和公式可得

最后n趋近于∞求极限,这里给出一个近似图,因为我们做图是不可能真的做出n为∞的图

[0,π]平均分成n个小区间,分点为,取,和式为

我们算到这里已经没办法计算下去了,因为这里没有上面的那个1到n的立方和公式之类的东西。所以一般来说,这种先求和,再求极限,非常困难。那么是否有便捷之法呢,答案是肯定的。

我们来看一个特例,考虑变速直线运动,[t1,t2]经过的路程可用定积分表示,[t1,t2]经过的路程可用位置函数S(t)表示为,所以

定理:如果函数f(x)在区间[a,b]上连续,则积分上限函数在[a,b]上可导,且

证明:

这里是根据定积分的性质2进行了区间拆分。

根据积分中值定理,将该式进行转化,得(介于x与x+∆x之间)

当∆x->0时,->x,则有,最终可得,得证

原函数存在定理:根据原函数定义,是f(x)的一个原函数

微积分基本定理(牛顿-莱布尼茨公式):如果函数F(x)是连续函数f(x)在区间[a,b]上的一个原函数,则有

简要证明:F(x)是f(x)的一个原函数,也是f(x)的一个原函数,于是

令x=a,得,这里是因为f(t)dt的定积分从[a,a]的面积为0,所以

令x=b,得,即,所以

例子

的原函数为,因为()'=

根据牛顿-莱布尼茨公式=1^3/4-0^3/4=1/4

sin x的原函数为-cos x

根据牛顿-莱布尼茨公式=-cos π+cos 0=1+1=2

定积分与和式极限

  • 计算

求解这一类问题,一般都是先展开,再利用自然数平方和公式,先求和再求极限

这种方法有很大的局限性,那就是自然数的次方公式,利用定积分的定义同样可以解这道题

像这种先求和再求极限的问题,我们称为和式极限。

  • 计算函数的图像围成的封闭图形的面积

将[0,1]区间分成n个长度相等的小区间,令,则面积为

[0,1]区间的定积分                            

  • 计算

这里的重点就在于我们需要将一般的表达式转换称和式极限表达式,再转换成左边的形式,再转换成右边的定积分,最后求出结果,这里的原函数为x^(3/2+1)/(3/2+1)=2x^(5/2)/5,再根据牛顿-莱布尼茨公式,得2/5-0=2/5

  • 计算

[a,b]区间的定积分                                

现在我们回到

这里我们设定(b-a)/n=2/n,2i/n-1=a+(b-a)i/n,则设定a=-1,b=1;1-x^2的原函数为x-x^3/3,根据牛顿-莱布尼茨公式,得(1-1/3)-(-1+1/3)=4/3

  • 计算

这里我们设定(b-a)/n=3/n,a+(b-a)i/n=2+3i/n,则设定a=2,b=5;1/x的原函数为lnx,根据牛顿-莱布尼茨公式,得ln5-ln2。

定积分应用-求平面曲线的弧长

元素法:我们之前所说的微元法又称为元素法。

用定积分表示所求量Q的条件

  1. Q是一个变量x的变化区间[a,b]有关的量
  2. Q对于区间[a,b]具有可加性
  3. 部分量∆Q的近似值可表示为

步骤

选取积分变量x,确定其变化区间[a,b]。把区间[a,b]分成n个小区间,取其中任一小区间[x,x+dx],求出相应于这个小区间的部分量∆Q的近似值,如果∆Q=f(x)dx+o(dx)(高阶无穷小),就把f(x)dx称为量Q的元素,记为dQ,dQ=f(x)dx,则

  • 计算平面曲线长,曲线的参数方程

相应于[a,b]上任一小区间[t,t+dt]的小弧度的长度(弦长),,∆s的近似值即弧长元素为ds,,则,这里我们舍弃了验证o(dx),o(dy),这个验证过程比较麻烦,我们就直接略过了,但是结果是正确的。

对于y=f(x)的情况,此时的参数方程为,弧长为

  • 举例,求函数在[-1,1]上的弧长

这里√(1+x^2)的原函数,我们可以采用分部积分法来求解,令u=√(1+x^2),dv=dx,则du=(√(1+x^2))'dx,v=x。

(√(1+x^2))'是一个复合函数求导,根据复合函数求导公式,令du=1+x^2,=1/[2√(1+x^2)],=2x,则(√(1+x^2))'=x/√(1+x^2)

移项化简得

这里,设x=tant,则dx=d(tant)=tan'tdt=sec^2tdt

∫(1/√(1+x^2))dx=∫(1/√(1+tan^2t))sec^2tdt=∫costsec^2tdt=∫sectdt=∫(cost/cos^2t)dt=∫(1/cos^2t)dsint=∫(1/1-sin^2t)dsint

令sint=θ化为∫(1/(1-θ^2))dθ=(1/2)∫(2/[(1+θ)(1-θ)])dθ=(1/2)∫[1/(1+θ)-1/(1-θ)]dθ=(1/2)∫1/(1+θ)dθ-(1/2)∫1/(1-θ)dθ,根据,则有(1/2)∫1/(1+θ)dθ-(1/2)∫1/(1-θ)dθ=(ln|1+θ|-ln|1-θ|)/2+C。根据对数公式以及,则(ln|1+θ|-ln|1-θ|)/2+C=ln√[(1+θ)/(1-θ)]+C=ln√[(1+sint)/(1-sint)]+C=ln|(1+sint)/cost|+C=ln|sect+tant|+C。因为sect=√(1+tan^2t),则ln|sect+tant|+C=lnl√(1+tan^2t)+tantl+C=lnl√(1+x^2)+xl+C

定积分的应用-用元素法求平面图形的面积

夹逼准则:如果数列 满足下列条件:

  1. 从某项起,即,当时,有

那么数列的极限存在,且

如果

  1. (去心邻域)(或|x|>M)时,m(x)≤f(x)≤n(x);

那么存在,且等于A

  • 求平面图形的面积,计算两条抛物线:所围成的图形的面积

其中红色曲线为,黄色曲线为。它们有两个交点O(0,0)和A(2,4)。其中小矩形BCIF的面积为为矩形的高BC,dx为矩形的宽CI。(考虑dx>0的情况)。这里我们知道,dS是不等于的,如果我们要使用元素法,就要判定它是否满足元素法的条件,∆S和dS是不是相差一个比dx高阶的无穷小

要证(考虑dx>0的情况)

因为

所以

易知

所以,只需证,只要证明了这两个式子成立,根据夹逼准则,就成立。

令C点的横坐标为

即证,我们将dx=0代入这两个式子,结果很明显它们都等于0,得证,满足元素法使用条件

,A=(2,4)

这里的原函数为4√2x^(3/2)/3-x^3/3

  • 计算抛物线与直线y=3x-6所围成的图形的面积

这里A=(1,-3),B=(4,6)

过A点垂直于x轴的直线将该图形分成了两半部分,我们需要将这两半部分的面积相加就可以了

其中的原函数为12x^(3/2)/3;的原函数为2x^(3/2)-3x^2/2+6x

定积分的应用-求连续型随机变量的概率

分布函数

我们翻开一本讲概率的书,会发现很多公式都是带有积分∫这个符号的,高等数学是一门基础课,后续的很多课程都会用到。虽然我们这里讲的是微积分的应用,但微积分的应用往往都会结合其他的学科,这样才能发挥出它实际的意义。一个实际的问题往往是多学科交叉的东西,很少有一种单学科的知识能够在实际当中应用的特别的好。

随机变量X的分布函数为F(x),若存在非负函数f(x),使得对于任意实数x,有,则X为连续型随机变量,f(x)为X的概率密度函数。

性质

  1. 首先概率密度分布函数f(x)≥0。
  2. 然后概率密度分布函数在-∞到+∞的积分是等于1的。这里是因为总的概率之和为1。
  3. 随机变量在x1和x2之间(x1,x2]的概率就为
  4. 如果概率密度分布函数f(x)在点x连续,则有F'(x)=f(x)

注意:P{X=a}=0,随机变量取某一个点的时候,它的概率为0。比如身高,通常来说我们都会认为它是连续型的随机变量,小明的身高精确值为1.75米,则P{X=1.75}=0,这似乎存在矛盾,其实概率为0的事件,它是有可能发生的。

正太分布

正太分布,其中表示均值,表示方差。

概率分布函数:

概率密度函数:时取最大值

由于计算过于复杂,我们会把概率分布函数做成一个表,通过查表计算

正太分布概率密度函数图

蓝色曲线表示均值为0,方差为1的正太分布;红色曲线表示均值为1,方差为4的正太分布。

如果要求蓝色曲线中(-1,1]的概率,就是将概率密度函数在(-1,1]上做定积分。

  • 计算概率:某连续型随机变量X的概率密度函数,求P{-1<X≤1}

这里由于e的指数-|x|一定是小于等于0的,所以我们要分成两部分来计算

微分方程

微分方程的意义

客观世界的描述:函数是客观事物内部联系在数量方面的反映。通过函数研究事物规律。那么如何获得函数?

  • 根据观测数据拟合

比如说运动轨迹,我要扔一个东西,就建立一个三维的空间坐标系,以我所在的地方为坐标原点。坐标系建立之后,我们就可以在每一个时间点上记录这个物体在坐标系中的位置。经过了一系列的记录之后,就会得到物体的运动轨迹。有了这个运动轨迹,我们就可以对它进行拟合,得到这个函数。

某产品价格与利润,我们可以对历史数据进行采样分析。

这里横坐标是产品的价格,纵坐标是利润。我们可以将这些历史数据点连接起来,从连接起来的平面曲线,我们可以看到它具有一定的趋势,或者说它和我们的一些函数长的比较像。这个时候我们可以溢出任何函数来拟合它。比如我们觉得它比较像二次函数ax^2+bx+c,虽然这里a、b、c我们还不知道,我们可以根据这些数据来找出二次函数的参数。这个过程就是一个拟合的过程。

  • 建立微分方程,求解微分方程,得到函数

我们来看一个刹车距离的问题,汽车在平直路上行驶,速度10m/s,突然发现前方20米处有一障碍物,司机立即刹车,汽车获得的加速度,问汽车能否撞到障碍物?

我们将这个题目用数学化的语言描述一下,这也是一个数学建模的过程,数学建模的第一步就是重述问题,也就是说我们会将原来的问题尽量转化为数学符号的描述。为了方便后续建模和求解,我们将这道题目进行这样的描述。我们知道速度是位移对时间的导数,同时我们又知道加速度是速度对时间的导数。这个时候,我们知道加速度a=-2,初始速度v0=10,要求位移s。即为已知,求s。

,这里v(t)是速度对时间的函数,亦是加速度的积分,这里-2的原函数为-2t

代入初始值v0=10,即为v(0)=-2t+C=10,此时t为0,求得C=10,则v(t)=-2t+10,由此可以看出,t=5时,速度减为0,汽车停止。

在这5秒中,汽车行驶的距离,,这里s(t)是位移对时间的函数,亦为速度的积分,这里-2t+10的原函数为。代入初始值s(0)=0,则,则,求得,所以会撞上障碍物。

总群自然生长模型:在没有限制的条件下,总群增长速度和数量呈一次函数关系。

我们用N表示总群的数量,则N对时间求导,表示一次函数关系,如何求得N(t)

微分方程应用领域:物理学、生物学、天文学、经济学、文化、人口、资源。如果我们所研究的问题是和变化率相关的问题,可以考虑使用微分方程。

求几种特定形式的微分方程的通解

现在我们来看一下什么叫微分方程

未知函数、未知函数的导数、自变量。例如(二阶微分方程)、(三阶微分方程),这里y是未知函数,x是自变量。我们主要探讨的是一阶微分方程。微分方程的解:函数。通解:解中含任意常数,常数个数等于阶数。特解:通过初始条件确定了常数值。

可分离变量的一阶微分方程

形如变量处于等号两端,两端积分。

  • 例子:求y'=3xy通解

y'=3xy    因为y'=dy/dx,分离变量,得

两边积分,得

齐次方程

形如,令,这里是根据两个函数的积的微分公式,再两边分别除以dx。

,分离变量,得,这里我们可以看到等式的左边不包含x,右边不包含u,变量处于等号两端

两端积分,得

  • 求解微分方程

化简,得,令

,化简,得,得    u-ln|u|+C2=ln|x|+C3,根据对数公式,这里C1代表常数,不分正负。即,即

可化为齐次的方程

令x=X+m,y=Y+n,则dx=dX,dy=dY

若方程组有解(),则

,令,则

令新变量h=ax+by,则

,这是可分离变量的方程

  • 例子,求微分方程(x+y)dx+(3x+3y-4)dy=0的通解

,显然

令h=x+y,,分离变量,得,即,两端积分并整理,得,这里化简得[3/2+1/(2-h)]dh=dx,3/2+1/(2-h)的原函数为3h/2+ln|2-h|,等式右边1的原函数为x。

  • 求微分方程(2x-5y+3)dx-(2x+4y-6)dy=0的通解

令x=X+m,y=Y+n,(2X-5Y+3+2m-5n)dX-(2X+4Y-6+2m+4n)dY=0

由于2/2≠-4/5,则

代入上式,得dY/dX=(2-5u)/(2+4u),由于Y=uX,则Y'=dY/dX=(uX)'=u'X+u=Xdu/dX+u=(2-5u)/(2+4u),则,即,此时变量处于等号两端,两端积分

代回变量,得方程通解为

一阶线性微分方程

形如

,此时变量处于等号两边,两端积分,得

当Q(x)≠0,将常数C换成未知函数u(x),,则,代入标准方程,得,所以,两端积分,得,故的通解为,即

  • 求微分方程的通解

这里面

思考:

  1. 求解过程技巧性高
  2. 不同形式对应不同的解法,难度高
  3. 是否存在一种容易掌握的方法,答案是有的。

利用Python求微分方程的通解

我们先来看一个最简单的微分方程f'(x)=1,结果我们已经知道为f(x)=x+C

# 该库支持符号运算,高精度计算,模式匹配,绘图,解方程,微积分计算,离散数学,概率统计
from sympy import *

if __name__ == "__main__":

    # 定义自变量
    x = symbols("x")
    # 定义函数
    f = symbols("f", cls=Function)
    # 描述微分方程f'(x)=1
    diffeq = Eq(f(x).diff(x), 1)
    # 输出结果
    print(dsolve(diffeq, f(x)))

运行结果

Eq(f(x), C1 + x)

现在我们来看一下上一小节的几个例子

  • y'=3xy
# 该库支持符号运算,高精度计算,模式匹配,绘图,解方程,微积分计算,离散数学,概率统计
from sympy import *

if __name__ == "__main__":

    # 定义自变量
    x = symbols("x")
    # 定义函数
    y = symbols("y", cls=Function)
    # 描述微分方程y'=3xy
    diffeq = Eq(y(x).diff(x), 3 * x * y(x))
    # 输出结果
    print(dsolve(diffeq, y(x)))

运行结果

Eq(y(x), C1*exp(3*x**2/2))

这里exp(3*x**2/2))表示

# 该库支持符号运算,高精度计算,模式匹配,绘图,解方程,微积分计算,离散数学,概率统计
from sympy import *

if __name__ == "__main__":

    # 定义自变量
    x = symbols("x")
    # 定义函数
    y = symbols("y", cls=Function)
    # 描述微分方程y^2+x^2y'=xyy'
    diffeq = Eq(y(x)**2 + x**2 * y(x).diff(x), x * y(x) * y(x).diff(x))
    # 输出结果
    print(dsolve(diffeq, y(x)))

运行结果

Eq(y(x), -x*LambertW(C1/x))

这里LambertW函数又称欧米茄函数或乘积对数,是f(w)=we^w的反函数,而且是多值的。

  • (x+y)dx+(3x+3y-4)dy=0

dy/dx=-(x+y)/(3x+3y-4)

# 该库支持符号运算,高精度计算,模式匹配,绘图,解方程,微积分计算,离散数学,概率统计
from sympy import *

if __name__ == "__main__":

    # 定义自变量
    x = symbols("x")
    # 定义函数
    y = symbols("y", cls=Function)
    # 描述微分方程(x+y)dx+(3x+3y-4)dy=0,即dy/dx=-(x+y)/(3x+3y-4)
    diffeq = Eq(y(x).diff(x), -(x + y(x)) / (3 * x + 3 * y(x) - 4))
    # 输出结果
    print(dsolve(diffeq, y(x)))

运行结果

Eq(y(x), x**2*(-C1*(3*C1/(3*C1 - 4) - 1)/(3*C1 - 4) + 3*C1/(3*C1 - 4) - 1)/(2*(3*C1 - 4)) + x**3*(-6*C1*(-3*C1/(3*C1 - 4) + 1)/(3*C1 - 4) - C1*(-6*C1*(-3*C1/(3*C1 - 4) + 1)/(3*C1 - 4) - 18*C1/(3*C1 - 4) + (3*C1/(3*C1 - 4) - 1)**2 + 6)/(3*C1 - 4) - 18*C1/(3*C1 - 4) + (3*C1/(3*C1 - 4) - 1)**2 + 6)/(6*(3*C1 - 4)**2) + x**4*(-54*C1*(3*C1/(3*C1 - 4) - 1)/(3*C1 - 4) - 18*C1*(-3*C1*(3*C1/(3*C1 - 4) - 1)/(3*C1 - 4) + 9*C1/(3*C1 - 4) + (-3*C1/(3*C1 - 4) + 1)*(3*C1/(3*C1 - 4) - 1) - 3)/(3*C1 - 4) - C1*(-54*C1*(3*C1/(3*C1 - 4) - 1)/(3*C1 - 4) - 18*C1*(-3*C1*(3*C1/(3*C1 - 4) - 1)/(3*C1 - 4) + 9*C1/(3*C1 - 4) + (-3*C1/(3*C1 - 4) + 1)*(3*C1/(3*C1 - 4) - 1) - 3)/(3*C1 - 4) + 162*C1/(3*C1 - 4) + 18*(-3*C1/(3*C1 - 4) + 1)*(3*C1/(3*C1 - 4) - 1) + (3*C1/(3*C1 - 4) - 1)*(-6*C1*(-3*C1/(3*C1 - 4) + 1)/(3*C1 - 4) - 18*C1/(3*C1 - 4) + (3*C1/(3*C1 - 4) - 1)**2 + 6) - 54)/(3*C1 - 4) + 162*C1/(3*C1 - 4) + 18*(-3*C1/(3*C1 - 4) + 1)*(3*C1/(3*C1 - 4) - 1) + (3*C1/(3*C1 - 4) - 1)*(-6*C1*(-3*C1/(3*C1 - 4) + 1)/(3*C1 - 4) - 18*C1/(3*C1 - 4) + (3*C1/(3*C1 - 4) - 1)**2 + 6) - 54)/(24*(3*C1 - 4)**3) + x**5*(-648*C1*(-3*C1/(3*C1 - 4) + 1)/(3*C1 - 4) - 108*C1*(-6*C1*(-3*C1/(3*C1 - 4) + 1)/(3*C1 - 4) - 18*C1/(3*C1 - 4) + (-3*C1/(3*C1 - 4) + 1)**2 + 2*(3*C1/(3*C1 - 4) - 1)**2 + 6)/(3*C1 - 4) - 6*C1*(-108*C1*(-3*C1/(3*C1 - 4) + 1)/(3*C1 - 4) - 18*C1*(-6*C1*(-3*C1/(3*C1 - 4) + 1)/(3*C1 - 4) - 18*C1/(3*C1 - 4) + (-3*C1/(3*C1 - 4) + 1)**2 + 2*(3*C1/(3*C1 - 4) - 1)**2 + 6)/(3*C1 - 4) - 324*C1/(3*C1 - 4) + 18*(-3*C1/(3*C1 - 4) + 1)**2 + (-3*C1/(3*C1 - 4) + 1)*(-6*C1*(-3*C1/(3*C1 - 4) + 1)/(3*C1 - 4) - 18*C1/(3*C1 - 4) + (3*C1/(3*C1 - 4) - 1)**2 + 6) + 36*(3*C1/(3*C1 - 4) - 1)**2 + 6*(3*C1/(3*C1 - 4) - 1)*(-3*C1*(3*C1/(3*C1 - 4) - 1)/(3*C1 - 4) + 9*C1/(3*C1 - 4) + (-3*C1/(3*C1 - 4) + 1)*(3*C1/(3*C1 - 4) - 1) - 3) + 108)/(3*C1 - 4) - C1*(-648*C1*(-3*C1/(3*C1 - 4) + 1)/(3*C1 - 4) - 108*C1*(-6*C1*(-3*C1/(3*C1 - 4) + 1)/(3*C1 - 4) - 18*C1/(3*C1 - 4) + (-3*C1/(3*C1 - 4) + 1)**2 + 2*(3*C1/(3*C1 - 4) - 1)**2 + 6)/(3*C1 - 4) - 6*C1*(-108*C1*(-3*C1/(3*C1 - 4) + 1)/(3*C1 - 4) - 18*C1*(-6*C1*(-3*C1/(3*C1 - 4) + 1)/(3*C1 - 4) - 18*C1/(3*C1 - 4) + (-3*C1/(3*C1 - 4) + 1)**2 + 2*(3*C1/(3*C1 - 4) - 1)**2 + 6)/(3*C1 - 4) - 324*C1/(3*C1 - 4) + 18*(-3*C1/(3*C1 - 4) + 1)**2 + (-3*C1/(3*C1 - 4) + 1)*(-6*C1*(-3*C1/(3*C1 - 4) + 1)/(3*C1 - 4) - 18*C1/(3*C1 - 4) + (3*C1/(3*C1 - 4) - 1)**2 + 6) + 36*(3*C1/(3*C1 - 4) - 1)**2 + 6*(3*C1/(3*C1 - 4) - 1)*(-3*C1*(3*C1/(3*C1 - 4) - 1)/(3*C1 - 4) + 9*C1/(3*C1 - 4) + (-3*C1/(3*C1 - 4) + 1)*(3*C1/(3*C1 - 4) - 1) - 3) + 108)/(3*C1 - 4) - 1944*C1/(3*C1 - 4) + 108*(-3*C1/(3*C1 - 4) + 1)**2 + 6*(-3*C1/(3*C1 - 4) + 1)*(-6*C1*(-3*C1/(3*C1 - 4) + 1)/(3*C1 - 4) - 18*C1/(3*C1 - 4) + (3*C1/(3*C1 - 4) - 1)**2 + 6) + 216*(3*C1/(3*C1 - 4) - 1)**2 + 36*(3*C1/(3*C1 - 4) - 1)*(-3*C1*(3*C1/(3*C1 - 4) - 1)/(3*C1 - 4) + 9*C1/(3*C1 - 4) + (-3*C1/(3*C1 - 4) + 1)*(3*C1/(3*C1 - 4) - 1) - 3) + (3*C1/(3*C1 - 4) - 1)*(-54*C1*(3*C1/(3*C1 - 4) - 1)/(3*C1 - 4) - 18*C1*(-3*C1*(3*C1/(3*C1 - 4) - 1)/(3*C1 - 4) + 9*C1/(3*C1 - 4) + (-3*C1/(3*C1 - 4) + 1)*(3*C1/(3*C1 - 4) - 1) - 3)/(3*C1 - 4) + 162*C1/(3*C1 - 4) + 18*(-3*C1/(3*C1 - 4) + 1)*(3*C1/(3*C1 - 4) - 1) + (3*C1/(3*C1 - 4) - 1)*(-6*C1*(-3*C1/(3*C1 - 4) + 1)/(3*C1 - 4) - 18*C1/(3*C1 - 4) + (3*C1/(3*C1 - 4) - 1)**2 + 6) - 54) + 648)/(3*C1 - 4) - 1944*C1/(3*C1 - 4) + 108*(-3*C1/(3*C1 - 4) + 1)**2 + 6*(-3*C1/(3*C1 - 4) + 1)*(-6*C1*(-3*C1/(3*C1 - 4) + 1)/(3*C1 - 4) - 18*C1/(3*C1 - 4) + (3*C1/(3*C1 - 4) - 1)**2 + 6) + 216*(3*C1/(3*C1 - 4) - 1)**2 + 36*(3*C1/(3*C1 - 4) - 1)*(-3*C1*(3*C1/(3*C1 - 4) - 1)/(3*C1 - 4) + 9*C1/(3*C1 - 4) + (-3*C1/(3*C1 - 4) + 1)*(3*C1/(3*C1 - 4) - 1) - 3) + (3*C1/(3*C1 - 4) - 1)*(-54*C1*(3*C1/(3*C1 - 4) - 1)/(3*C1 - 4) - 18*C1*(-3*C1*(3*C1/(3*C1 - 4) - 1)/(3*C1 - 4) + 9*C1/(3*C1 - 4) + (-3*C1/(3*C1 - 4) + 1)*(3*C1/(3*C1 - 4) - 1) - 3)/(3*C1 - 4) + 162*C1/(3*C1 - 4) + 18*(-3*C1/(3*C1 - 4) + 1)*(3*C1/(3*C1 - 4) - 1) + (3*C1/(3*C1 - 4) - 1)*(-6*C1*(-3*C1/(3*C1 - 4) + 1)/(3*C1 - 4) - 18*C1/(3*C1 - 4) + (3*C1/(3*C1 - 4) - 1)**2 + 6) - 54) + 648)/(120*(3*C1 - 4)**4) + C1 - C1*x/(3*C1 - 4) + O(x**6))
  • (2x-5y+3)dx-(2x+4y-6)dy=0

dy/dx=(2x-5y+3)/(2x+4y-6)

运算时间较长,这里只给出代码

# 该库支持符号运算,高精度计算,模式匹配,绘图,解方程,微积分计算,离散数学,概率统计
from sympy import *

if __name__ == "__main__":

    # 定义自变量
    x = symbols("x")
    # 定义函数
    y = symbols("y", cls=Function)
    # 描述微分方程(2x-5y+3)dx-(2x+4y-6)dy=0,即dy/dx=(2x-5y+3)/(2x+4y-6)
    diffeq = Eq(y(x).diff(x), (2 * x - 5 * y(x) + 3) / (2 * x + 4 * y(x) - 6))
    # 输出结果
    print(dsolve(diffeq, y(x)))
# 该库支持符号运算,高精度计算,模式匹配,绘图,解方程,微积分计算,离散数学,概率统计
from sympy import *

if __name__ == "__main__":

    # 定义自变量
    x = symbols("x")
    # 定义函数
    y = symbols("y", cls=Function)
    # 描述微分方程dy/dx-y(x+1)=x+1
    diffeq = Eq(y(x).diff(x) - y(x) * (x + 1), x + 1)
    # 输出结果
    print(dsolve(diffeq, y(x)))

运行结果

Eq(y(x), C1*exp(x*(x/2 + 1)) - 1)

微分方程的数值解——欧拉法

无法求解的微分方程:其实有大量的微分方程,它们是没有代数解析式的,比如

我们可以用Python来验证一个,比如y'=2xy+1

# 该库支持符号运算,高精度计算,模式匹配,绘图,解方程,微积分计算,离散数学,概率统计
from sympy import *

if __name__ == "__main__":

    # 定义自变量
    x = symbols("x")
    # 定义函数
    y = symbols("y", cls=Function)
    # 描述微分方程y'=2xy+1
    diffeq = Eq(y(x).diff(x), 2 * x * y(x) + 1)
    # 输出结果
    print(dsolve(diffeq, y(x)))

运行结果

Eq(y(x), (C1 + sqrt(pi)*erf(x)/2)*exp(x**2))

这里erf是一个误差函数,也被称为高斯误差函数,通常是用一个积分的形式来表示,也就是说对这个微分方程来说,最终的通解是没有解析式的。

微分方程的数值解:实际应用中,知道函数的值即可。

比如标准正太分布概率密度函数

计算随机变量的取值在(-1,1]之间的概率P(-1<X≤1),这其实是在求概率密度函数在(-1,1]的定积分

但是这个定积分是没有表达式的,f(x)的一个原函数记为g(x),则P(-1<X≤1)=g(1)-g(-1)可以获取数值解,方程的数值解通常可以满足实际需求。

欧拉法

考虑一般问题

变量离散化,自变量x的变化量∆x=,这个n是可以大也可以小的,是我们自己设定的(太小的话,可能精度达不到要求;如果太大的话,可能计算量会比较大,计算时间比较长)

泰勒公式:

忽略及以上的项,这里实际上就相当于微分

递推计算

这里y(x1)是用y1来代替的

  • 例子:

我们令n=50,则有

由此递推下去,可以用程序来实现。

这里的蓝色曲线是使用欧拉方法递推出来的值,黄色曲线是准确值。这里明显是有误差的,我们可以通过增大n的取值来改善

我们令n=100,则有                

y1=y0-2x0y0h=1/e-2*1*1*0.04/e=0.3384

y2=y1-2x1y1h=0.3384-2*(1+0.04)*0.3384*0.04=0.3102

欧拉法误差分析

泰勒公式        

截断误差项        表示多项式的最低次数,通过截断误差项,我们可以看出,把h取的越小(即n越大),误差越小。我们这里考虑的是截断误差,没有考虑累计误差。

若某算法的局部截断误差为,则称该算法有p阶精度。欧拉法具有1阶精度。p越大,精度越高,如何提高精度。是有办法的,这是后面的龙格-库塔法。

使用Python实现欧拉法

我们还是来实现上一小节的第一个例子

import numpy as np
import matplotlib.pyplot as plt

if __name__ == "__main__":

    # y'=x,x0=-2,y0=3
    a = -2
    b = 2
    x0 = -2
    y0 = 3
    n = 500
    h = (b - a) / n
    y = []
    y.append(y0)
    # 自变量离散化
    # np.arange返回一个有终点和起点的固定步长的排列,如[1,2,3,4,5],起点是1,终点是6,步长为1
    # 这里也就是说最后一个值不是b,而是b-h
    x = np.arange(a, b, h)
    for i in x:    # 迭代过程
        tmpy = y0 + i * h
        y.append(tmpy)
        y0 = tmpy
    # 这里我们需要最后一个值是b
    x1 = np.arange(a, b + h, h)
    # 函数表达式
    y1 = 0.5 * x1**2 + 1
    plt.plot(x1, y, lw=3, label="euler_value")
    plt.plot(x1, y1, lw=3, label="exact_value")
    plt.legend()
    plt.show()

运行结果

现在我们再来实现第二个例子

import numpy as np
import matplotlib.pyplot as plt
import math

if __name__ == "__main__":

    # y'=-2xy,x0=1,y0=1/e, y=exp(-x**2)
    a = 1
    b = 5
    x0 = 1
    y0 = 1 / math.e
    n = 100
    h = (b - a) / n
    y = []
    y.append(y0)
    # 自变量离散化
    # np.arange返回一个有终点和起点的固定步长的排列,如[1,2,3,4,5],起点是1,终点是6,步长为1
    # 这里也就是说最后一个值不是b,而是b-h
    x = np.arange(a, b, h)
    for i in x:    # 迭代过程
        tmpy = y0 - 2 * i * y0 * h
        y.append(tmpy)
        y0 = tmpy
    # 这里我们需要最后一个值是b
    x1 = np.arange(a, b + h, h)
    # 函数表达式
    y1 = math.e**(-x1**2)
    plt.plot(x1, y, lw=3, label="euler_value")
    plt.plot(x1, y1, lw=3, label="exact_value")
    plt.legend()
    plt.show()

运行结果

微分方程的数值解——龙格-库塔法

欧拉法单步迭代,只用了端点信息(我们预测n+1个点的时候,我们只用了第n个点的信息),现在我们考虑增加数据点数,使误差多项式阶数提高(提高精度)。

递推公式

这里k1跟欧拉法是一样的,k1、k2可认为是两个不同点处的增量。下一个值的增量部分是两个增量的加权平均。a1和a2是我们要求的,代表权重。

求解参数a1、a2、b、c

将k2在点处进行二元函数的泰勒展开

二阶龙格-库塔法

在点处按一元函数泰勒展开

根据多元复合函数求导法则

要使截断误差达到2阶,即,则有

三个方程,4个未知数,解不唯一,我们可以任取一组解

三阶龙格-库塔法

8个未知系数,6个方程,用了三个点的数据,达到三阶精度,一种常见递推式

四阶龙格-库塔法

用4个点的信息,达到四阶精度

13个未知系数,11个方程,一种常见递推式

利用Python实现龙格-库塔法

我们这里要实现的是四阶龙格-库塔法

import numpy as np
import matplotlib.pyplot as plt
import math

if __name__ == "__main__":

    # y'=-2xy,x0=1,y0=1/e,y=exp(-x**2)
    a = 1
    b = 5
    x0 = 1
    y0 = 1 / math.e
    n = 50
    h = (b - a) / n
    y = []
    y.append(y0)
    # 自变量离散化
    # np.arange返回一个有终点和起点的固定步长的排列,如[1,2,3,4,5],起点是1,终点是6,步长为1
    # 这里也就是说最后一个值不是b,而是b-h
    x = np.arange(a, b, h)
    for i in x:         # 迭代过程
        k1 = h * (-2 * i * y0)
        k2 = h * (-2 * (i + 0.5 * h) * (y0 + 0.5 * k1))
        k3 = h * (-2 * (i + 0.5 * h) * (y0 + 0.5 * k2))
        k4 = h * (-2 * (i + h) * (y0 + k3))
        tmpy= y0 + 1 / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
        y.append(tmpy)
        y0 = tmpy
    # 这里我们需要最后一个值是b
    x1 = np.arange(a, b + h, h)
    # 函数表达式
    y1 = math.e**(-x1**2)
    plt.plot(x1, y, lw=3, label="runge-kutta_value")
    plt.plot(x1, y1, lw=3, label="exact_value")
    plt.legend()
    plt.show()

运行结果

对比同样n=50的时候,欧拉法的图像

我们可以看到欧拉法的误差是非常明显的。

使用Python库函数求微分方程的数值解

import numpy as np
import matplotlib.pyplot as plt
import math
from scipy.integrate import odeint

if __name__ == "__main__":

    # y'=-2xy,x0=1,y0=1/e,y=exp(-x**2)
    a = 1
    b = 5
    x0 = 1
    y0 = 1 / math.e
    n = 50
    h = (b - a) / n
    y = []
    y.append(y0)
    # 自变量离散化
    # np.arange返回一个有终点和起点的固定步长的排列,如[1,2,3,4,5],起点是1,终点是6,步长为1
    # 这里也就是说最后一个值不是b,而是b-h
    x = np.arange(a, b, h)
    for i in x:         # 迭代过程
        k1 = h * (-2 * i * y0)
        k2 = h * (-2 * (i + 0.5 * h) * (y0 + 0.5 * k1))
        k3 = h * (-2 * (i + 0.5 * h) * (y0 + 0.5 * k2))
        k4 = h * (-2 * (i + h) * (y0 + k3))
        tmpy= y0 + 1 / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
        y.append(tmpy)
        y0 = tmpy
    # 这里我们需要最后一个值是b
    x1 = np.arange(a, b + h, h)
    # 函数表达式
    y1 = math.e**(-x1**2)
    def f(x, y):
        # 微分方程表达式
        return -2 * x * y
    y2 = odeint(f, 1 / math.e, x1)
    plt.plot(x1, y, lw=3, label="runge-kutta_value")
    plt.plot(x1, y1, lw=3, label="exact_value")
    plt.plot(x1, y2, lw=3, label="ode_value")
    plt.legend()
    plt.show()

运行结果

常见微分方程数学建模

传染病的微分方程模型

数学建模:根据实际问题,建立数学模型。然后对这个模型进行求解。达到解决问题的目的。一般过程为:问题重述(将问题用数学化的语言表述)->合理假设(一个实际的问题影响因素非常多,我们在研究一个问题,往往会抓住主要因素,忽略掉次要因素)->建立模型(处理方法,什么方法处理什么问题,需要积累)->模型求解(不同的模型对应不同的算法)->模型检验(根据常识或者是专业知识来判定它的结果的正确性和合理性或者根据模型的结果进行一些预测(预测类模型))

传染病的自由增长模型

模型假设:我们的主要目的是为了得到传染病患者的人数和时间t的一个关系。

  1.  单位时间内每个人有效接触(足以传染)人数为a
  2. 病人数为x(t)
  3. 在∆t内,病人数增加量为x(t)∆ta,则有x(t+∆t)-x(t)=x(t)∆ta,两边同时除以∆t,∆t->0取极限,得到,即x'(t)=x(t)a

令初始值为x(0),解该方程得,这里采用分离变量法,dx/dt=xa得dx/x=adt,两边同时积分∫dx/x=∫adt,得ln|x|=at+C1,得x=(+/-)e^(at+C1),最终x=Ce^(at)

,我们可以看出,随着时间增加,病人数将无限增长,不符合实际情况。

传染病SI模型

因为传染病的自由增长模型不符合实际情况,所以我们需要改进我们的方法,改变假设条件,不断的把各种因素考虑的更周全一些,让最终结果符合实际情况。当我们考虑的因素越来越多的时候,我们的模型会越来越复杂。所以我们采取的策略是从简单到复杂,层层递进的方式。

SI模型,S代表的是易感染者,也可以理解成健康人,I表示感染者,也叫病人。

模型假设:

  1. 每个病人每天有效接触人数为a,称为日有效接触率
  2. 总人数不变,设为N
  3. 所有人分为易感染者(健康人)和感染者(病人)两类
  4. 符号定义
N 总人数
s(t) t时刻健康人所占比例
i(t) t时刻病人所占比例
a 每个病人每天有效接触人数

我们的目标就是求出s(t)和i(t),这里s(t)+i(t)=1

构建模型:

  1. 每个病人每天可感染人数为as(t)
  2. Ni(t)个病人每天可传染人数为Ni(t)as(t)
  3. 由2可得Ni'(t)=Ni(t)as(t),即i'(t)=i(t)as(t),则有i'(t)=i(t)a(1-i(t)),假定初始值
  4. 求解该方程得

用Python直接求解该微分方程i'(t)=i(t)a(1-i(t))

# 该库支持符号运算,高精度计算,模式匹配,绘图,解方程,微积分计算,离散数学,概率统计
from sympy import *

if __name__ == "__main__":

    # 定义自变量
    t = symbols("t")
    # 定义函数
    i = symbols("i", cls=Function)
    a = 3
    # 描述微分方程i'(t)=i(t)a(1-i(t))
    diffeq = Eq(i(t).diff(t), i(t) * a * (1 - i(t)))
    # 输出结果
    print(dsolve(diffeq, i(t)))

运行结果

Eq(i(t), 1/(C1*exp(-3*t) + 1))

传染病初期取,我们可以这么来考虑,比如一个城市的人口是1000000,对应的是1个,表示第一个病人刚刚出现,图形如下

这里横坐标表示时间t,纵坐标表示感染病人所占比例i。由图可以看出,i在很长时间内,变化是很小的;但是过了一个时间点之后就会爆发式增长。所以能够在尽早的初期发现传染病对于传染病的控制来说,有非常重大的意义。这里SI模型没有考虑到病人可治愈,还需要做进一步的改进。

传染病SIS模型

模型假设

  1. 新增假设:每天治愈人数比例为u
  2. u代表平均治愈率,代表平均治愈时间
  3. a  每个患者每天有效接触人数
  4. 感染期内每个患者感染的平均健康人数

微分方程:

变形得

对于这个微分方程,一方面我们可以求出这个微分方程的通解,再根据初始条件,求它的特解。另外一种就是微分方程的数值求解,必须要给定a、b一个具体的值才可以求解。a、b是可变的,对于不同的a、b,对应的求解微分方程的函数图像有很大的差别。

现在我们对b进行分析,当b≤1时,i(t)单调递减。取a=0.2,b=0.5,

取a=0.2,b=0.5,

上面两个图是利用龙格-库塔法或者Python库函数做出的。

关键数据,它极有可能对应i(t)的极小值或者极大值。我们先取<

取a=0.2,b=2,

取a=0.2,b=1.2,

现在我们取>,取a=0.2,b=1.2,

我们可以看到无论>(递增)还是<(递减),它都是趋近于的。

从图像中可以看到,减小b的值有利于我们抗击传染病,为了减小b的值有两种方案,一是减小a(隔离传染源切断传染途径),另一个是增大u(提高医疗水平)。SIS模型没有考虑感染者会死亡以及感染者治愈之后产生抗体不会再被感染这两个因素。

传染病的SIR模型

模型假设

  1. 修改假设:每天治愈人数和死亡人数之和的比例为u
  2. 新增假设:治愈者都会产生抗体
  3. r(t)~t时刻治愈和死亡患者之和所占的比例

微分方程组

                    初始值

这个方程组是没有解析解的,所以我们需要数值解

,r0=0

第68天的时候,感染比例最大约为38%

,r0=0

第38天,感染比例最大约为55%

由图中我们可以看到,增大了a,将感染比例最大值提前了,感染比例的最大值也增大了。

,r0=0

第92天,感染比例最大,约为15%。

传染病控制方法

  1. 减小a:隔离、阻断传染途径
  2. 增加u:提高医疗水平

利用Python实现求微分方程组的数值解

龙格-库塔法

之前我们已经知道了龙格-库塔法解微分方程

变形得

龙格-库塔法求微分方程组的数值解

递推式

传染病SIR模型的微分方程组

对应

Python实现

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

if __name__ == "__main__":

    # 求微分方程组的数值解
    a = 0.5
    u = 0.15
    b = 200
    i0 = 0.000001
    s0 = 1 - i0
    n = b * 100
    i = i0
    s = s0
    h = b / n
    idata = []
    sdata = []
    idata.append(i0)
    sdata.append(s0)
    t = np.arange(0, b, h)
    for j in t:        #迭代过程
        k1 = a * s * i - u * i
        m1 = -1 * a * s * i
        k2 = a * (s + h / 2 * m1) * (i + h / 2 * k1) - u * (i + h / 2 * k1)
        m2 = -1 * a * (s + h / 2 * m1) * (i + h / 2 * k1)
        k3 = a * (s + h / 2 * m2) * (i + h / 2 * k2) - u * (i + h / 2 * k2)
        m3 = -1 * a * (s + h / 2 * m2) * (i + h / 2 * k2)
        k4 = a * (s + h * m3) * (i + h * k3) - u * (i + h * k3)
        m4 = -1 * a * (s + h * m3) * (i + h * k3)
        tmpi = i + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
        tmps = s + h / 6 * (m1 + 2 * m2 + 2 * m3 + m4)
        idata.append(tmpi)
        sdata.append(tmps)
        i = tmpi
        s = tmps
    t1 = np.arange(0, b + h, h)
    # 描述微分方程
    def func(Iiitial_value, var, para):
        ii, ss, rr = Iiitial_value
        pa, pu = para
        return np.array([pa * ss * ii - pu * ii, -1 * pa * ss * ii, pu * ii])
    p = odeint(func, (i0, s0, 1 - i0 - s0), t1, args=([a, u],))
    idata_1 = p[:, 0]
    sdata_1 = p[:, 1]
    rdata_1 = p[:, 2]
    rdata = 1 - np.array(idata) - np.array(sdata)
    plt.plot(t1, idata, lw=3, label="i(t)")
    plt.plot(t1, sdata, lw=3, label="s(t)")
    plt.plot(t1, rdata, lw=3, label="r(t)")
    plt.plot(t1, idata_1, lw=1, label="i(t)_1")
    plt.plot(t1, sdata_1, lw=1, label="s(t)_1")
    plt.plot(t1, rdata_1, lw=1, label="r(t)_1")
    plt.legend()
    plt.show()

运行结果

SIR模型待改进的地方,一般来说传染病有一个潜伏期,但是SIR模型没有考虑潜伏期,这是可改进的地方,另外我们在模型当中使用的两个重要参数a和u,我们是设为固定值的,但实际上这两个参数很可能是变化的。这些因素都可能会影响到我们最终利用这个模型来预测实际情况时的准确性。

线性回归

最小二乘法

客观世界的描述:函数是客观事物内部联系在数量方面的反映,通过函数来研究事物的规律。

如何获得函数:

  1. 微分方程模型
  2. 数据拟合

拟合与回归

  1. 数据之间可由某种函数或方程近似表示,但具体参数未知。
  2. 可以用拟合的方式求得这些参数。
  3. 回归是求得这些参数的一种具体方法或模型。
  4. 其他常见的拟合方法:ARIMA(差分整合移动平均自回归模型,又称为整合移动平均自回归模型)、VAR(向量自回归模型)、神经网络等。

线性回归

数据呈线性关系,通过线性回归模型求得具体参数,比如某视频课程阅读量与上架时间关系

横坐标表示上架时间,纵坐标表示阅读量(累计阅读量)。我们可以很明显的发现,这些点近似分布在一条直线上,这就是一种典型的数据呈线性关系。对这样的数据,我们可以用线性回归来拟合。

 

分析:根据散点图观察,存在线性关系,则f(t)=kt + b,     k、b为待求参数。找到合适的k、b,使点尽可能均匀分布在直线的两侧。通常来说,我们的数据里面都是带噪声的,数据有噪声,几乎不可能让数据点都在一条直线上,不能在直线上,我们就尽可能的让这些点分布在直线的两侧。

最小二乘法:找到合适的k、b,使残差平方和最小。

残差:在上面的公式中,为拟合的数据,为实际的数据,这两个数据对应的差值,我们就称为残差。我们要求的就是min(M)

M是k、b的二元函数,

要求M的极值的必要条件,即M对k、b的偏导数要等于0

            即        

这是一个复合函数求导,,u=

整理合并,分离k、b后,得

其中

由此我们可以看出,上面就是一个二元线性方程组,最终得

最终函数为,图像如下

这里阅读量的拟合值是在一条直线上的

多变量情况

这就是求多元函数的极值,如果变量非常多,那么我们使用这种方法计算起来就会非常麻烦,有没有什么简单一点的方法呢,答案是有的。

使用线性代数实现最小二乘法

关于用到的线性代数的内容请参考线性代数整理 ,但这里需要补充一点,

向量

  • 例子,3个自变量,5个样本

这是一个三元线性函数,我们要求的是a0、a1、a2、a3

该式子中,第一部分是5个样本的函数值y向量(即真实数据向量),是拟合数据向量(矩阵乘以向量为向量),则结果为残差m。我们的目标是使向量m的平方最小(即m中的元素平方再求和的值最小),则等价于最小。

即函数最小。

类比极限存在的必要条件,函数L对向量a的导数为0.这里

Python实现

  • 函数关系

数据

from numpy.linalg import inv   #求矩阵的逆
from numpy import dot
import pandas
import numpy as np

if __name__ == "__main__":

    # 最小二乘法 y=1.5+0.5*x1+0.8*x2-0.2*x3
    filedata = pandas.read_csv("/Users/admin/Downloads/data.csv")
    Xdata = np.array(filedata)
    print(Xdata)
    X1 = Xdata[:, [0, 1, 2]]
    y = Xdata[:, 3]
    m, n = X1.shape
    # 生成数据全为1的向量
    b = np.ones(m)
    X = np.c_[b, X1]
    print(X)
    a = dot(dot(inv(dot(X.T, X)), X.T), y)
    print(a)

运行结果

[[ 4.  10.  24.   6.7]
 [10.  14.  26.  12.5]
 [ 3.  19.  25.  13.2]
 [10.  14.  20.  13.7]
 [ 3.  13.  22.   9. ]]
[[ 1.  4. 10. 24.]
 [ 1. 10. 14. 26.]
 [ 1.  3. 19. 25.]
 [ 1. 10. 14. 20.]
 [ 1.  3. 13. 22.]]
[ 1.5  0.5  0.8 -0.2]

这里我们需要注意的是,这里的a0、a1、a2、a3是未知的,是我们拿到数据的时候需要求的,这里只是为了方便比较结果才给出了多元线性函数。

线性回归的假设与检验

  1. 回归方程的形式应该是准确的。如果因变量和自变量的关系不是线性的,而是二次函数关系,如果使用线性回归的方法来做,肯定是有问题的。
  2. 全部的影响因素都应该考虑到。假设我们的回归方程的影响因素有5个,而我们只找到了3个,采集数据也只拿到了这3个影响因素的数据,那么得到的回归方程,显然它的解释力度是不够的。
  3. 不相干的因素不被纳入模型。
  4. 采集的样本数据应当准确。通常来说我们所能获取的数据都是存在噪声的,不可能是纯净的数据,所以我们只能尽量去让它准确。、

如何使回归方程更准确

  1. 一致性。随着样本量的不断增加,估计的值就会越来越接近被估计的总体的参数。也就是说随着样本量的不断增加,咱们估计的值就不断的接近真实值。一般一致性是在大样本的条件下去说的,小样本一般不会谈一致性。
  2. 无偏性。通常来说我们对整体的某一个参数的估计是通过采样的方式来完成的。由于采样是存在着偶然因素或者随机性的,每一次抽出样本都是有差别的,所以我们根据样本得到的估计值也不会完全相同,单凭某一次抽样的样本是不太具有说服力的。所以我们往往会进行多次抽样,因为每一次抽样,我们都会有一个估计值,在多次抽样之后,对这些估计值取一个平均,我们可以理解为求期望。如果这个期望值它总体参数是一致的,那么这种情况我们就说是满足无偏性。
  3. 有效性。有效性指的是对同一个总体参数,如果有多个无偏估计量,那么标准差最小的估计量更有效,实际上这里说的是估计量和总体参数的离散程度。离散程度越小的估计量,我们说它更有效。

假设我们要估计一个城市的20-50岁的人的平均身高,这个城市的人口比较多,我们要做这个事情不太可能把这个城市的20-50岁的人都来测量一下他们的身高,然后取平均值。我们会采取抽样的方法,我们会随机去抽取200个满足20-50岁条件的人,通过测量他们的身高,然后得到这200个人的平均值。如果我们将这个平均值作为对整体的一个估计,一般来说是不太合适的。因为单凭一次的抽样,偶然性太大,偏差会很大,所以我们会采取多次抽样的方式,假设我们进行了100次的抽样,那么我们每一次抽样之后都会得到一个平均值,就会得到100个平均值。那么我们将这100个平均值再求一次平均值,那么最后得到的这个平均值,如果说它就等于我们的总体的参数,就等于这个城市20-50岁的人的平均身高的话,那我们就说估计量满足无偏性。同时这100个数存在着方差,方差代表的是对均值的离散程度。我们假设有2个团队,他们分别都在做这个事,现在我们假设A、B两个团队,他们的估计量都满足无偏性,然而A团队估计量的方差比B团队估计量的方差更小,那么这个时候我们就说A团队的估计量更有效。

模型假设

  1. 线性性。比如因变量和自变量满足这样一个关系y=0.5*x1+0.6*x2^2,很明显它不满足线性性。再比如y=0.5*x1+0.6*x2+0.5*x1*x2,这样的关系它也是不满足线性性的。不满足线性性,我们就不能直接用线性回归来处理。
  2. 解释变量与随机误差项不相关。如果解释变量与随机误差项相关的话,我们就说模型存在内生性问题。内生性问题是线性回归中比较复杂的问题。比如说我们的模型没有考虑到全部的影响因素,而模型之外的影响因素又和模型之内的影响因素存在着一定的相关关系,那么这种情况就会造成解释变量与随机误差项存在相关。又比如数据的测量方面存在误差,这也可能导致解释变量与随机误差项存在相关。
  3. 随机误差项无序列相关性。指的是随机误差项本身是不存在自相关的。如果随机误差项本身存在自相关,那么使用线性回归来处理效果就不会很好。随机误差项存在自相关很可能是我们遗漏了时间这个因素,当然也还可能有其他因素。
  4. 随机误差项的方差为常数。如果随机误差项的方差不为常数的话,会影响模型的有效性。
  5. 随机误差项服从0均值的正太分布的。这条假设的意义就在于我们可以在此基础之上做检验。

线性性

  1. 可以通过散点图初步判断。
  2. F检验。
  3. t检验。
  4. 不满足线性性可尝试的处理方法:对数变换、倒数变换、平方根变换。做了这些变换后可能依然存在不满足线性性,我们就要考虑一下是不是采集数据存在问题,或者要考虑一下是否要更换模型。如果不满足线性性,我们要考虑不用线性回归来处理。

F检验

  1. 判定模型整体是否显著
  2. 假设各个系数都为0,计算F统计值,得出概率,根据概率判断。
  3. F值服从自由度为(p,n-p-1)的F分布。

这里SSR叫回归平方和,表示是因变量的变异转移回归模型中所含的p个自变量所能够解释的部分,SSE称为误差平方和,又叫残差平方和。

t检验

  1. 判定具体某个解释变量与因变量之间是否存在线性关系。
  2. 根据某个系数为0,计算t值,得出概率,根据概率判断是第i个变量的偏回归系数,是其标准误,标准误不是标准差,标准误=标准差/√n,n代表样本量。

模型的解释能力评价

解释能力指的是因变量它总的变异当中,可以由回归模型中自变量解释的部分所占的比例,那么这个比例越大,显然模型的解释能力越强

  1. 决定系数。决定系数可以用来评价解释能力。不过决定系数有不好的地方,如果我向模型当中再添加一些其他的变量,即使这个变量和我们的因变量y是没有什么关系的,但是决定系数还是会增大。所以决定系数不够客观,不够准确。
  2. 校正的决定系数。校正的决定系数克服了决定系数存在的问题。我们在看这个校正的决定系数的时候,是看它越大,对模型的解释能力是越强的,但最大不会超过1。
  3. 剩余标准差。
  4. 赤池信息准则。

决定系数

因变量的总变异中可由回归模型中自变量解释部分所占比例,越大越好,不足:变量增加,就会增加,即使增加的变量没有统计学意义。

校正的决定系数

残差分析

  1. 残差与自变量的散点图判断该自变量与因变量的线性关系。
  2. 残差间是否相互独立——D-W检验
  3. 自变量个数小于4,统计量大于2,基本可以判定残差相互独立。
  4. 判断方差齐性:绘制标准差与因变量的标化预测值散点图。
  5. 判断正太分布:绘制标准化残差的直方图、茎叶图、正太概率分布图

多重共线性问题

多重共线性指的是自变量间存在明显的共线性,这个时候我们是不能直接使用最小二乘法进行回归分析的

  1. 矩阵的行列式值接近0.如果多重共线性越严重,那么这个矩阵的行列式就会越接近0,特别严重的情况下,我们是无法求出系数向量的。但不太严重的情况下,我们还是可以求出系数向量a,但是这个时候是不稳定的。比如当我们改变X中的数据做微小的改变,那么系数向量a可能就会有很大的变化,这就是叫做不稳定。这一点对应的就是有效性不够好。
  2. 判定:容忍度、方差膨胀因子(VIF)、特根值、条件指数
  3. 处理方法:
    1. 逐步回归
    2. 岭回归
    3. 主成分分析
展开阅读全文
加载中

作者的其它热门文章

打赏
0
1 收藏
分享
打赏
0 评论
1 收藏
0
分享
返回顶部
顶部