Python基础语法学习总结

实验目的

学习Python基本语法

实验场地与设备

线上

实验方式

阅读教程与程序设计

实验设计

Python语言基础

\[ 图1.1 Python基础语法学习实验设计 \]

实验内容

1. Python语法总结

1.1 Python基本语法

（1） 基本语句

① 首先是输入输出语句，输入语句比较简单为name=input()，基本输出语句为print(),拼接输出使用逗号。

② 注释采用# 进行书写

③ 代码风格：Python采用的是缩进式代码风格，所以对于复制粘贴比较不友好

④ 条件判断语句：if 条件1 :...elif 条件2 : ... else : ...

⑤ 循环语句：

第一种是for循环：for x in []: for x in ...: 循环就是把每个元素代入变量x，然后执行缩进块的语句

第二种是while循环：while 条件判断语句 : break、continue和java中用法相同

（2） 数据类型

① 整数：对于很大的数，很难数清楚0的个数。Python允许在数字中间以_分隔。

② 浮点数：允许使用科学计数法定义

③ 字符串：在Python没有严格要求''和""的区别在，也就是说没有区分字符和字符串使用二者没有任何区别。

  • 转义符和Java中保持一致
  • Python允许用r''表示''内部的字符串默认不转义

④ 布尔值：

在Python中要注意：True、False要注意开头首字母大写。 可以进行与、或、非的运算，运算符分别为：and，or，not

⑤ 空值：空值用None表示，意义与Java中的null相同。

⑥ list：

list是Python内置的一种数据类型，list是一种有序的集合，可以随时添加和删除其中的元素。此数据类型在Java的实用类中有封装。list和数组很像，声明方式如下：

1
classname = ['老六','老八','老九']

想要调取其中的某个元素也和数组一致，赋值修改等也相同
下面列举一下list的ADT

1
2
3
4
5
6
list:
append('Elem')  # 在末尾添加新的元素
insert(i,'Elem') # 将元素插入指定位置
pop() # 删除末尾元素
pop(i) # 删除i处的元素
len(list) # list列表的长度

list允许混合类型，也允许list嵌套，从而出现多维数组。

⑦ tuple

tuple被称为元组，其最大的特点就是不可修改，声明方式如下：

1
classname = ('老六','老八','老九')

tuple在定义时要确定元素个数，这里有一个问题，在定义只有一个元素的tuple时，Python语法会认为这是一个小括号，因此在定义一个元组的tuple时，要加一个,避免歧义。

1
t=(1,)

⑧ 字典（dict）

字典全称为dictionary，在Java实用类中叫hash map。其由键值对（key-value）组成，查找速度快。 下面是一种初始化方法：

1
d = {'Michael': 95, 'Bob': 75, 'Tracy': 85}

也可以放入指定的key中：

1
d['Adam'] = 67

查找value:

1
d['Adam']

key与value是多对一的关系，key需要是一个不可变对象保证key做hash运算后的唯一性。如果多次对某个key赋值，后边的value会覆盖前面的value 提供了几个函数：

  1. 通过in来判断key是否在dict中，返回值为布尔值，格式为：key in dict
  2. get()方法，dict.get('key',空返回值)key不存在时返回空返回值，空返回值可自定义，如果没有定义的话返回None
  3. pop()方法，删除key，如果有value也一并删除，格式为pop('key')

⑨ 集合（set）

set是一组key的集合,集合特点；无序性、确定性、互异性 要创建一个set，需要提供一个list作为输入集合：

1
s = set([1, 2, 3])
  • 方法： add(key)添加一个新的元素 remove(key)删除一个元素
  • 两个set可以做交运算和并运算： 交运算：s1&s2 并运算：s1|s2

（3） 理解变量

在Python中变量仅仅是一个一个字母，变量与所对应的值之间的关系靠指针联系起来的。所以很重要的一点就是：当我们使用变量时，更多的要关注变量指向的东西，他可能是值，也可能是一个函数，也可能是一个变量

1.2 模块

（1） 模块导入

1
import numpy as np

（2） 模块下载

模块下载有比较复杂的方法，也有比较傻瓜式的。先说复杂的，使用Python中自带的pip包管理工具，用命令：

1
pip install numpy

但是使用pip需要事先了解要导的包的名字，而且不能批量导入，而且在Python编程里也有编程一分钟，导包一小时的说法。pip下载第三方库的源可能会很慢或者失效，需要会自己添加国内的高速镜像。

傻瓜式的导包，例如在pycharm中可以直接在代码中写出自己需要的包，然后交给pycharm自己去下载，或者用Anaconda提前构建好的Python的库环境。

1.3 函数式编程

（1） 函数

① 函数定义

在Python中定义函数为，def 函数名(参数):然后，在缩进块中编写函数体，函数的返回值用return语句返回。
如果没有return语句，函数执行完毕后也会返回结果，只是结果为None。return None可以简写为return。

1）空函数：

1
2
def nop():
    pass

在这里pass作为占位符，表示跳过，也可以用在if的缩进块。

2）参数限制：

1
2
if not isinstance(x, (int, float)):
      raise TypeError('bad operand type')

实际上参数限制就是定义一个报错，isinstance()判断数据类型，如果不是就提出一个错误。 作为一个弱类型语言，定义这一步是很有必要的，有助于读懂代码。

3）返回值：

Python允许返回多个值，其返回的实际上是一个tuple元组，但是也可以用两个变量接收。

② 参数定义

在Python中函数参数的定义也比较灵活，提供位置参数、默认参数、可变参数、关键字（key）参数等

1）位置参数：位置参数指的是参数在传入时，实参和形参有着严格的位置对应关系，为常用参数形式。

2）默认参数：默认参数是指在位置参数的基础上为其添加默认值，有默认值的参数为默认参数，没有默认值的参数为必选参数 基本定义形式如下：

1
2
3
def my_def(a,b=1):
    a=b+1
    return

需要注意的是：

  • 默认参数必须在必选参数后边，否则会无法辨认是否输入必选参数，从而报错。
  • 默认参数的默认值一定是不变对象，由于Python中的变量定义为指针指向，会导致可变对象值发生变化

3）不可变对象有：数值类型、字符串、tuple元组、None等

4）可变参数：可变参数指的是参数的数目不固定，定义形式：

1
2
3
4
5
def my_function(*v):
    sum = 0
    for vi in v:
        sum+=vi
    return sum

在可变参数中传入的所有参数将作为一个tuple被接收，该tuple的变量名为函数在定义时的形参名，定义时的需要在参数名前加一个*。

5）关键字（key）参数

此处的关键字和c语言中的关键字并不是一个意义，而是在dict中的key的意义。即在传递参数时，同时传递键（key）和值(value),Python会自动封装为一个dict。

1
2
3
def my_function(**v):
    print(v)
    return

6）命名关键字参数

在关键字参数上，进一步限制传入的key的命名，就有了命名关键词参数：

1
2
def person(name, age, *, city, job):
    print(name, age, city, job)

这里需要一个*区分位置参数与命名关键字参数，如果在这之前有可变参数，那么就不需要加*。
命名关键字参数必须传入参数名，这和位置参数不同。如果没有传入参数名，调用将报错：

7）参数组合

在一个函数中使用多个参数要保证其中的顺序，依次为：必选参数、默认参数、可变参数、命名关键字参数和关键字参数。

1
2
def onefunction(a,b,c=0,*args,job,city,**kw):
    pass

tips：

  • 使用*args和**kw是Python的习惯写法。
  • 可变参数和关键字参数有一点层级的感觉，中间包裹的是命名关键字参数这个比较尴尬的参数。

③ 递归函数

写法与Java相同。

（2） 实用方法

① 切片

切片是一个针对tuple和list方便地取元素的方法，语法规则：

1
L[起始坐标:终止坐标:步长]

当起始坐标为0时可以省略；步长为1时可以省略。

② 迭代

迭代是循环的增强，但是想要弄清迭代，需要知道两件事：一个是能不能迭代，一个是迭代出的数据是什么

想要知道一个数据能否迭代可以通过一个函数来完成：

1
2
3
from collections.abc import Iterable
L=[1,2,3]
isinstance(L,Iterable)

迭代出的是什么，和要迭代的对象的储存方式，要特殊记忆一下dic。

③ 列表生成器

一种快捷生成list的方式，一个例子如下：

1
[x * x for x in range(1, 11)]

如果想要筛选生成的值，可以在for后加上if作为筛选条件，注意这里是筛选条件， 因此这里和平时的if else并不是一个东西。

1
[x * x for x in range(1, 11) if x % 2 == 0]

④ 生成器

生成器是一种惰性的计算方式。包含yield关键字，当一个函数包含yield关键字时，他就成了一个generator函数。yield在generator函数中起到了一个return的作用，即到yield便返回。 在调用时，使用一个变量接受一个generator对象。使用next()函数依次获得下一个返回值。

⑤ 迭代器

区分Iterable和Iterator

Iterable是可迭代的，是直接可用于for循环的。包括dict、list、tuple、set、str、grenerator。 Iterator是迭代器，是直接可用于next()函数的，生成器都是Iterator对象，集合数据类型可以通过iter()获取Interator对象。

（3） 函数式编程

函数式编程是一种面向过程的编程思想，实际上是将复杂问题转化为一个个函数。

在Java的函数定义中，除去void类型不返回值，其余的都需要返回值。因此也就经常存在，使用一个变量接受函数值：

1
2
3
4
int function(x,y){
return;
}
int a=function(x,y);

那么是不是存在一种可能，我们可以将函数嵌套，让函数调用函数，让函数返回函数，彻底抛弃变量？

抛弃变量、只有函数就是彻底的函数式编程

① 理解高阶函数

之前有过变量名和值的理解，在Python中变量名和值是一个指针指向的关系。同理，函数名和函数也是这样的，函数名也是一个变量。也就是说，我们可以通过函数名，拿到函数体。也就是说函数名是什么并不重要，我们看中的是函数体。

绘图1

那么设想一种情况，现在我们定义了函数f2，那么我可以随便写一个函数，然后返回一个变量f2，那么实际上我就拿到了函数体。

1
2
3
4
5
def f2(a,b):
    return a+b
def f3():
    return f2
print(f3()(1,2))
image-20220909173741530

然后我们在设想另一种情况，现在我们定义了另一种情况，我们在一个函数中写了一个f1作为局部变量，那么我就可以传入变量f2，然后就相当于传入了函数体。

1
2
3
4
5
def f2(a,b):
    return a+b
def f1(a,b,f):
    return f(a,b)
print(f1(1,2,f2))

现在就可以进行一个区分：

  • f代表函数名，是变量
  • f()代表数值，是函数的返回值，返回值是一个量

高阶函数，就是让函数的参数能够接收别的函数。

实用的几个函数，有必要查表即可

② 返回函数

同上文理解，只不过是将一个函数嵌套入了另一个函数

③ lambda表达式

与Java中语法相同，目的是为了简化返回函数嵌套

1.4 面向对象编程

（1）类和对象

创建类：语法如下

1
class 类名(继承的类):

python的类非常随意，几乎可以不定义就能用。在类中自带有一个构造函数__init__(),此函数可以重新定义

生成对象：

1
a=A()

（2）访问权限

如果要让内部属性不被外部访问，可以把属性的名称前加上两个下划线__，在Python中，实例的变量名如果以__开头，就变成了一个私有变量（private），只有内部可以访问，外部不能访问。

此外，__ __这种变量都是特殊变量，在不清楚的时候不要随便乱改

（3）继承和多态

和Java中的思想完全相同

（4）常用变量和方法

① __slots__

用这个变量可以起到参数列表的功能，可以在一定程度上限制参数的变量名，用turple进行限定

② @property

注解编程，可以起到一个简化定义setter和getter函数的作用。@property注解在getter方法上，然后会自动生成 @函数名.setter 的注解，但是要注意的一点是，在getter中就不能使用函数名作为自身的调用值，否则会出现无限的调用，产生爆栈。

③ 多继承

与Java相同

⑤ __str__:和Java中的toString方法相同

1.5 错误调试

（1）错误处理

参照Java中，对比来学习即可：

两种方法，一是尝试，二是抛出，尝试采用：

1
2
3
4
5
6
try:
	pass
except baseexception  :
	pass
finally:
	pass

抛出采用raise关键字

（2）测试

① 断言：assert的意思是，表达式n != 0应该是True，否则，根据程序运行的逻辑，后面的代码肯定会出错。

如果断言失败，assert语句本身就会抛出AssertionError

② 断点：在强大IDE的辅助下，使用断点调试应该是最简单的。

2.实践

2.1 石头剪子布

使用random包中的random函数和条件控制语句，模拟两个电脑互相猜拳：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
import random

def win(pc,cc):
    if (cc==1 and pc==2) or (cc==2 and pc==3)or(cc==3 and pc==1):
        print("电脑一输")
    elif pc==cc:
        print("平")
    else:
        print("电脑二输")

def computer_choice():
    cc=random.randint(1,3)
    return cc
def show(pc,cc):
    print("电脑一的出招为",pc)
    print("电脑二的出招为",cc)

if __name__=="__main__":
    while(True):
        pc=computer_choice()
        cc=computer_choice()
        show(pc,cc)
        win(pc,cc)
image-20220914212607801

改进提升一下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
import random

def win(pc,cc):
    if (cc==1 and pc==2) or (cc==2 and pc==3)or(cc==3 and pc==1):
        print("玩家输")
    elif pc==cc:
        print("平")
    else:
        print("玩家赢")

def computer_choice():
    cc=random.randint(1,3)
    return cc
def str(cc):
    if cc==1:
        return '石头'
    elif cc==2:
        return '剪刀'
    else:
        return '布'
def show(f,pc,cc):
    print("电脑一的出招为",f(pc))
    print("电脑二的出招为",f(cc))

if __name__=="__main__":
    while(True):
        cc=computer_choice()
        print("请输入：1.石头 2.剪刀 3.布")
        pc=input()
        show(str,pc,cc)
        win(pc,cc)
image-20220914213324805

2.2 ATM模拟

通过类和对象简单的设计了一个ATM取钱模拟器

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
import Account


class ATM(object):
    def __init__(self,money,accounts):
        self.money=money
        self.accounts=accounts
    @property
    def money(self):
        return self._money;
    @money.setter
    def money(self,value):
        self._money=value
    @property
    def accounts(self):
        return self._accounts
    @accounts.setter
    def accounts(self,value):
        self._accounts=value

    def searchId(self,id):
        for account in self.accounts:
            if account.id==id:
                return account
    def lode(self):
        print('请输入账号id')
        id = input()
        account1 = self.searchId(id)
        print('请输入密码')
        password = input()
        if password == account1.password:
            print("欢迎", account1.name)
        return account1
    # 存钱
    def save_money(self):
        account=self.lode();
        print("请输入要存入的数目")
        saveMneyValue=input()
        print('存款成功')
        account.remain=int(account.remain)+int(saveMneyValue)
        print('您的账户余额为',account.remain)
        self.money=self.money+int(saveMneyValue)

    # 取钱
    def withdraw_money(self):
        account=self.lode()
        print('请输入要取出的数目')
        withdrawMoneyValue=input()
        if account.remain > withdrawMoneyValue:
            account.remain=int(account.remain)-int(withdrawMoneyValue)
            print('取款成功，您的账户余额为',account.remain)
        else:
            print('您的账户余额不足')
        self.money=int(self.money)-int(withdrawMoneyValue)

    def __str__(self):
        print("当前ATM中有金额",self.money,"元")

if __name__=="__main__":
    # atm1=ATM(1000)
    # atm1.__str__()
    # atm1.ave_money(200)
    # atm1.__str__()
    # atm1.withdraw_money(200)
    # atm1.__str__()
    accounts=[]
    for i in range(5):
        name=input()
        id = input()
        password=input()
        remain=input()
        accounts.append(Account.account(name, id, password, remain))
    atm2=ATM(10000,accounts)
    atm2.save_money()
    atm2.withdraw_money()
1
2
3
4
5
6
7
8
class account(object):
    def __init__(self,name,id,password,remain):
        self.name=name
        self.remain=remain
        self.password=password
        self.id=id

    __slots__ = ('name','remain','password','id')
image-20220914214759256

2.3 圣诞树画图

使用Python自带的turtle包，进行圣诞树绘制：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
import turtle

screen = turtle.Screen()
screen.setup(375, 700)

circle = turtle.Turtle()
circle.shape('circle')
circle.color('red')
circle.speed('fastest')
circle.up()

square = turtle.Turtle()
square.shape('square')
square.color('green')
square.speed('fastest')
square.up()

circle.goto(0, 280)
circle.stamp()

k = 0
for i in range(1, 13):
    y = 30 * i
    for j in range(i - k):
        x = 30 * j
        square.goto(x, -y + 280)
        square.stamp()
        square.goto(-x, -y + 280)
        square.stamp()

    if i % 4 == 0:
        x = 30 * (j + 1)
        circle.color('red')
        circle.goto(-x, -y + 280)
        circle.stamp()
        circle.goto(x, -y + 280)
        circle.stamp()
        k += 3

    if i % 4 == 3:
        x = 30 * (j + 1)
        circle.color('yellow')
        circle.goto(-x, -y + 280)
        circle.stamp()
        circle.goto(x, -y + 280)
        circle.stamp()

square.color('brown')
for i in range(13, 17):
    y = 30 * i
    for j in range(2):
        x = 30 * j
        square.goto(x, -y + 280)
        square.stamp()
        square.goto(-x, -y + 280)
        square.stamp()
turtle.mainloop()
image-20220914215352995

3.总结

Python作为一个弱类型语言，是有他的弊端的，在一些需要数据类型转换和严格控制数据类型的情况下，会非常难受。而Python最大的优势在于有大量的库，这些库在特定的编程领域会非常便利。Python本身的语言具有极强的灵活性，而灵活性的言外之意就是规范性很难确定。因此，Python的重点是将第三方包为我所用，在数值计算中发挥他最大的作用。

实验二 numpy和matplotlib包学习

实验目的

  1. 掌握基本的numpy对象及其对应方法
  2. 掌握常用的numpy数学函数，学习查找numpy帮助文档
  3. 重点学习numpy线性代数方法
  4. 掌握matplotlib的绘图对象关系
  5. 掌握基本的绘制图形的方法，包括绘制、属性设置、子图
  6. 能够通过查阅文档、示例，画出复杂图像

实验场地与设备

实验室4074

实验方式

阅读教程与程序设计

实验设计

常用计算库

实验内容

概述

1.numpy概述

NumPy是使用Python进行科学计算的基础包。它包含如下的内容：

  • 一个强大的N维数组对象。
  • 复杂的（广播）功能。
  • 用于集成C / C ++和Fortran代码的工具。
  • 有用的线性代数，傅里叶变换和随机数功能。

除了明显的科学用途外，NumPy还可以用作通用数据的高效多维容器。可以定义任意数据类型。这使NumPy能够无缝快速地与各种数据库集成。

API文档：Array objects — NumPy v1.23 Manual

2.matplotlib概述

Matplotlib是一个Python 2D绘图库，它以多种硬拷贝格式和跨平台的交互式环境生成出版物质量的图形。

API文档：API Reference — Matplotlib 3.6.0 documentation

具体实验

1.数组存取以及关于内存位置相同的思考

通过切片、整数列表、列表生成式等方法尝试数组取出，并对切片内存相同提出解决方案。联想到线性代数，取出子数组对于判断海森矩阵是否正定有很大意义。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
# 实现对数组元素的存取
import numpy as np #导入numpy包


print('-----------数组存取----------')
a=np.arange(9).reshape(3,3) #建立0-8的一个3*3的二位数组
print('输出数组a：')
print(a)

print('-----------取出第一行----------')
# 取出第一行所有元素
b=a[0,:]
print("输出数组b：")
print(b)
# 取子数组
print('-----------取出子数组----------')
x=a[:2,:2]
print(x)

print('-----------验证内存相同----------')
# 修改数组b内一元素，并输出a,b
b[0]=11
print('修改后的数组a，b：')
print(a)
print(b)

print('-----------如何处理内存相同----------')
# 可见通过切片存取数组，只是一个可修改的视图，二者共用同一内存空间。我考虑到两种方法：一是使用列表生成器，二是进行深拷贝
# 列表生成器取出a第一列
print('-----------列表存取法----------')
l=[i for i in range(a.shape[1])]
print('生成的list：')
print(l)
b=a[0,l]
print("输出数组b：")
print(b)
# 修改数组b内一元素，并输出a,b
b[0]=3
print('修改后的数组a，b：')
print(a)
print(b)

# 深拷贝方法（我个人认为这个更简单）
# 先用切片取出第一行所有元素
print('-----------深拷贝方法----------')
b=a[0,:]
print("输出数组b：")
print(b)
# 进行拷贝
c=b.copy()
# 修改数组b内一元素，并输出b,c
c[0]=5
print('修改后的数组b,c：')
print(b)
print(c)

输出结果：

image-20220925191202086

而在内存存储中，ndarray是以头地址加地址步长构成的，所以可以在处理连续的存取时，就可以通过指针直接映射到一个视图上。而整数序列不能保证连续分布，就只能拷贝一份。

2.结构数组练习

新建一个dtype:student。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
import numpy as np
# 结构数组类似c中的strut，和面向对象中的对象数组也比较相似
# 定义结构数组
studentType=np.dtype({
    'names':['name','age','id','grades'],
    'formats':['S20','i8','S20','f8']
})
# 创建studentType类型的对象
a=np.array([('lisi',11,'202011000',90.2)],
            dtype=studentType)
# 输出a
print(a)
# 这时出现了一个问题：每个字符串类型的数据都自动补齐了一个b，观察一下属性列表
print(a.dtype)
# 发现没有什么问题
# 查阅资料得知，b表示字符串为bytes类型，是字节字符串，在字符的前面会加一个b。numpy中没办法将bytes直接转化为str类型，会出现如下报错：
'''
print(a['name'],decode())
AttributeError: 'numpy.ndarray' object has no attribute 'decode'
'''
# 结构数组的优点在于赋予了数组中各个数值意义，那么也就有相应的存取方案
# 例如，取出姓名
print(a['name'])
# 判断name对应的numpy属性
print(a['name'].dtype)

输出结果： image-20221002130955868

3.自定义ufunc函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
import numpy as np

def triangle_wave(x,c,c0,hc):
    x=x-int(x)
    if x>=c:r=0.0
    elif x <c0: r=x/c0*hc
    else:r=(c-x)/(c-c0)*hc
    return r

# 调用一个高阶函数frompyfunc
triangle_ufunc=np.frompyfunc(triangle_wave,4,1)
# 在这里第二个参数指的是输入到高阶函数的返回值中的函数的参数个数，第三个参数是返回函数的返回参数个数
# 新建一个数组
a=np.random.rand(9).reshape(3,3)
print(a)
# 尝试调用
y=triangle_ufunc(a,0.6,0.4,1.0)
# 输出y
print(y)
image-20221002192935954

4.numpy实现梯度下降法求最小值

使用numpy实现梯度下降，具体见下注释：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
import numpy as np

# 函数对多项式函数求导函数
'''
分析：多项式函数求导数公式：
y=a*x^n
y‘=n*a*x^(n-1)
思路：在numpy中的有多项式函数构造器，本质上其利用的就是一个list列表，那么操作list列表就达到了求导的目的；
算法流程：
1.获取np.array对象的shape，并将shape元组中的的长度取出，设为count
2.声明一个临时变量list
3.对np.array的内容进行迭代，按照顺序，ndarray第1个元素乘count-1（count-1减去的是常数项那一位，刚好得到是首个元素的次数），将结果添加到list中，count每次循环-1
4.返回ndarray，其中要将list进行切片（目的是将常数项去除，保证使用构造器时次数正确）
'''
def dxfun1(x):
    m=x.shape
    count=int(m[0])
    list=[]
    for i in np.nditer(x):
        list.append(i*(count-1))
        count=count-1
    return np.array(list[0:len(list) - 1])
# 构造多项式函数
def dxfun(x):
    return np.poly1d(dxfun1(x))
# 理想最速梯度下降算法：
'''
理想最速梯度下降：
梯度下降需要解决的两个任务：寻找搜索方向和搜索步长，理想状态下搜索方向就是最快的沿着函数曲线下降，一阶导数表明了运动的趋势，二阶导数表明了运动的速度
所以当搜寻最小值时，那么一阶导数就要为负，而二阶导数要为负梯度，才能保证其下降速度是最快的
搜索步长：最佳步长可对lambda求导数，得到一个表达式

算法步骤：
1.给定一个varepsilon作为可容忍的误差，选定初始点x
2.根据梯度下降和最优步长移动点x，更新x位置
3.判断x位置是否满足可容忍的误差，如满足则输出最小值，如不满足则继续进行迭代
'''
def GD(x,varepsilon,a):
    min1=0
    while(True):
        if np.poly1d(dxfun1(a))(x)<=varepsilon:
            min1 = np.poly1d(a)(x)
            return min1
        else:
            ddx=dxfun1(dxfun1(a))
            lambda_1=1/np.poly1d(ddx)(x)**2
            x=x-lambda_1*np.poly1d(ddx)(x)


if __name__=='__main__':
    # a=np.linspace(1,5,5)
    a=np.array([1,0,4])
    # 输出原函数：
    print('原函数：')
    print(np.poly1d(a))
    # 输出导函数
    print('导函数')
    print(dxfun(a))
    print('最小值为：')
    print(GD(1,0.001,a))

    # 拟合多项式
    a_1 = np.array([1, 0, 0])
    f = np.poly1d(a)
    x=np.linspace(0,4*np.pi,40)
    y=f(x)*np.sin(x)+5*np.cos(x)
    z = np.polyfit(x, y, 6)
    f1=np.poly1d(z)
    z1 = f1(x)
    print(GD(20, 0.001, z))

输出结果： image-20221002130924057

5.梯度下降法可视化

使用散点图记录中间过程，用折线图画出函数轮廓得到如图：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
import numpy as np
import matplotlib.pyplot as plt

# 函数对多项式函数求导函数
'''
分析：多项式函数求导数公式：
y=a*x^n
y‘=n*a*x^(n-1)
思路：在numpy中的有多项式函数构造器，本质上其利用的就是一个list列表，那么操作list列表就达到了求导的目的；
算法流程：
1.获取np.array对象的shape，并将shape元组中的的长度取出，设为count
2.声明一个临时变量list
3.对np.array的内容进行迭代，按照顺序，ndarray第1个元素乘count-1（count-1减去的是常数项那一位，刚好得到是首个元素的次数），将结果添加到list中，count每次循环-1
4.返回ndarray，其中要将list进行切片（目的是将常数项去除，保证使用构造器时次数正确）
'''
def dxfun1(x):
    m=x.shape
    count=int(m[0])
    list=[]
    for i in np.nditer(x):
        list.append(i*(count-1))
        count=count-1
    return np.array(list[0:len(list) - 1])
# 构造多项式函数
def dxfun(x):
    return np.poly1d(dxfun1(x))
# 理想最速梯度下降算法：
'''
理想最速梯度下降：
梯度下降需要解决的两个任务：寻找搜索方向和搜索步长，理想状态下搜索方向就是最快的沿着函数曲线下降，一阶导数表明了运动的趋势，二阶导数表明了运动的速度
所以当搜寻最小值时，那么一阶导数就要为负，而二阶导数要为负梯度，才能保证其下降速度是最快的
搜索步长：最佳步长可对lambda求导数，得到一个表达式

算法步骤：
1.给定一个varepsilon作为可容忍的误差，选定初始点x
2.根据梯度下降和最优步长移动点x，更新x位置
3.判断x位置是否满足可容忍的误差，如满足则输出最小值，如不满足则继续进行迭代
'''
def GD(x,varepsilon,a):
    min1=0
    while(True):
        if np.poly1d(dxfun1(a))(x)<=varepsilon:
            min1 = np.poly1d(a)(x)
            plt.scatter(x, min1)
            return min1
        else:
            ddx=dxfun1(dxfun1(a))
            lambda_1=1/np.poly1d(ddx)(x)**2
            x=x-lambda_1*np.poly1d(ddx)(x)
            y=np.poly1d(a)(x)
            plt.scatter(x,y)


if __name__=='__main__':
    # a=np.linspace(1,5,5)
    a=np.array([1,0,4])
    # 输出原函数：
    print('原函数：')
    y1=np.poly1d(a)
    print(y1)
    # 输出导函数
    print('导函数')
    print(dxfun(a))
    print('最小值为：')
    print(GD(10,0.001,a))
    # 绘制二次函数y1图像
    x1=np.linspace(-10,10,20)
    plt.plot(x1,y1(x1))
    # 拟合多项式
    '''
    a_1 = np.array([1, 0, 0])
    f = np.poly1d(a)
    x=np.linspace(0,4*np.pi,40)
    y=f(x)*np.sin(x)+5*np.cos(x)
    z = np.polyfit(x, y, 6)
    f1=np.poly1d(z)
    z1 = f1(x)
    print('拟合计算最小值')
    print(GD(13, 0.001, z))
    plt.figure(x,z1)'''
    plt.show()

二次函数梯度下降可视化，如图所示：

image-20221002152415012

自定义函数f1：

image-20221002154352417

但是，如果出的试点选择为0和5，那么得到的结果为：

image-20221002160351783image-20221002160300703

经过测试可知，选取初始点极其重要，不好的选点在非凸优化中很容易引起陷入局部最优，从而找不到全局最优。

6.折线图

采用鸢尾花数据集，但是仅观察数据集中各个数据的特征，不做聚类分析，首先是三种鸢尾花的不同特征的比较：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
import matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets

iris =datasets.load_iris()
'''
通过查看load_iris源代码，可以看出，一共是三种花：'setosa', 'versicolor', 'virginica'，各有四种属性：'sepal length (cm)', 'sepal width (cm)','petal length (cm)', 'petal width (cm)'
'''
print('查看品类名称：')
print(iris.target_names)
print('查看特征名称：')
print(iris.feature_names)
# 整理各种数据，分为三类画图
y1=[]
y2=[]
y3=[]
for i in range(iris.data.shape[0]):
    if iris.target[i]==0:
        y1.append(iris.data[i,:])
    elif iris.target[i]==1:
        y2.append(iris.data[i, :])
    else:
        y3.append(iris.data[i, :])
# 建立一个和list长度相同的数组x，便于画图
x1=np.linspace(1, len(y1), len(y1))
x2=np.linspace(1, len(y2), len(y2))
x3=np.linspace(1, len(y3), len(y3))
# 将list y转化为ndarray
y1=np.array(y1)
y2=np.array(y2)
y3=np.array(y3)

# 针对不同类型分别绘制折线图
# sepal length (cm)
plt.figure(1)
plt.subplot(2,2,1)
plt.plot(x1,y1[:,0],label="setosa")
plt.plot(x2,y2[:,0],label="versicolor")
plt.plot(x3,y3[:,0],label="virginica")
plt.ylabel("sepal length (cm)")
plt.xlabel("x")
plt.legend() #用于生成图例
plt.title("sepal length (cm)")
# sepal width (cm)
plt.subplot(2,2,2)
plt.plot(x1,y1[:,1],label="setosa")
plt.plot(x2,y2[:,1],label="versicolor")
plt.plot(x3,y3[:,1],label="virginica")
plt.ylabel("sepal width (cm)")
plt.xlabel("x")
plt.legend()
plt.title("sepal width (cm)")
# petal length (cm)
plt.subplot(2,2,3)
plt.plot(x1,y1[:,2],label="setosa")
plt.plot(x2,y2[:,2],label="versicolor")
plt.plot(x3,y3[:,2],label="virginica")
plt.ylabel("petal length (cm)")
plt.xlabel("x")
plt.legend()
plt.title("petal length (cm)")
# petal width (cm)
plt.subplot(2,2,4)
plt.plot(x1,y1[:,2],label="setosa")
plt.plot(x2,y2[:,2],label="versicolor")
plt.plot(x3,y3[:,2],label="virginica")
plt.ylabel("petal width (cm)")
plt.xlabel("x")
plt.legend()
plt.title("petal width (cm)")
plt.tight_layout() #自适应分布
plt.show()

在做到后边时，突然想到其实在做完筛选后，其实可以添加一个排序，那么就会使图像非常清晰可见，见右图

image-20221002175529645image-20221002184429334

7.柱状图

1
2
3
4
5
6
7
8
# 提取setosa中的sepal length的范围,绘制柱形图
y1_sl,s=np.histogram(y1[:,0],bins=[4,4.5,5,5.5,6])
fig, ax = plt.subplots()

l = ['4-4.5','4.5-5','5-5.5','5.5-6']
l=np.array(l)
ax.bar(range(len(l)), y1_sl,tick_label=l)
plt.show()

image-20221002183734784

操作时发现，matplotlib 2.x版本中的bar是不支持使用字符串列表作为x，y两轴的，但是3.x版本就允许了，所以在使用包时要注意观察官方文档中给出包版本。

8.散点图

1
2
3
4
5
6
7
8
9
# 绘制散点图
plt.clf()
plt.scatter(y1[:,0],y1[:,1],label="setosa")
plt.scatter(y2[:,0],y2[:,1],label="versicolor")
plt.scatter(y3[:,0],y3[:,1],label="virginica")
plt.legend()
plt.xlabel(iris.feature_names[0])   #x轴名称
plt.ylabel(iris.feature_names[1])   #y轴名称
plt.show()

image-20221002190949615

从图中可以看出不是很好区分，那么可以换一下其他属性继续作图：

image-20221002191244939

9.饼图

1
2
3
4
5
6
# 绘制饼图
plt.figure(figsize=(4,4))
plt.pie(y1_sl,labels=['4-4.5','4.5-5','5-5.5','5.5-6'],autopct='%.2f%%',explode=(0.2, 0, 0, 0),)
zhfont1 = matplotlib.font_manager.FontProperties(fname="SourceHanSerifSC-Light.otf")
plt.title('setosa中的sepal length各个区间占比',fontproperties=zhfont1)
plt.show()

image-20221002192645105

10.3D绘图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D

# 生成数据x
x=np.arange(-4,4,0.125)
# 生成数据y
y=np.arange(-4,4,0.125)
# 对x、y数据执行网格化
x,y= np.meshgrid(x,y)
# 定义z
z1 = np.exp(-x**2 - y**2)
z2 = np.exp(-(x - 1)**2 - (y - 1)**2)
z=(z1+z2)/2
# 绘制3d图像
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(x,y,z,cmap=plt.get_cmap('rainbow'))
plt.xlabel("x")
plt.ylabel("y")
ax.set_zlabel("z")

fig.show()

image-20221002194419179

总结

经过对这两个库的学习发现，这两个库的功能都相当强大。对于numpy来说要理解numpy中ndarray的使用方法，切片、迭代使用频率都比较高。matplotlib主要是理解他的artist对象，理解后其实可以通过官网的example直接套用模板就可以做出很好的效果。理解artist对象后在对官网示例进行编辑、美化，为自己所用。

实验三 数据探索分析

实验目的

  1. 使用矩阵计算和统计量进行EDA
  2. 使用可视化技术进行EDA
  3. 掌握数据相似性和相异性计算

实验场地与设备

实验室4074

实验方式

阅读教程与程序设计

实验设计

要尽可能多的使用可视化方法EDA

实验内容

导入企鹅数据集

1
2
3
4
5
6
7
8
# coding=UTF-8
# This Python file uses the following encoding: utf-8
import matplotlib
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

df = pd.read_csv('../penguins_size.csv')

企鹅数据集的特征及解释：

1
2
3
4
5
6
7
8
9
10
11
'''
变量                 类型       含义
----------------    ------- ---------------------
species             integer	企鹅种类 (Adelie, Gentoo, Chinstrap)
island	            integer	所在岛屿 (Biscoe, Dream, Torgersen)
bill_length_mm	    double	嘴峰长度 (单位毫米)
bill_depth_mm	    double	嘴峰深度 (单位毫米)
flipper_length_mm	integer	鰭肢长度 (单位毫米)
body_mass_g	        integer	体重 (单位克)
sex	                integer	性别
'''

首先验证是否有缺失值：

1
2
# 观察缺失值
print(df.isnull().sum())

结果如下，发现有缺失值，然后直接删除缺失值行，再次观察缺失值。

1
2
3
4
5
6
7
8
species               0
island                0
culmen_length_mm      2
culmen_depth_mm       2
flipper_length_mm     2
body_mass_g           2
sex                  10
dtype: int64
1
2
3
4
# 直接删除缺失值行
df = df.dropna()
# 再次观察缺失值
print(df.isnull().sum())

现在我想要了解各种企鹅的数量，并通过饼状图可视化：

1
2
3
4
5
6
7
8
9
10
# 各种企鹅的数量
species_count=df['species'].value_counts()
print(species_count)
# 绘制饼状图
plt.figure(figsize=(4,4))
labels=species_count.index
plt.pie(species_count,labels=labels,autopct='%.2f%%',explode=(0,0,0.1))
zhfont1 = matplotlib.font_manager.FontProperties(fname="../SourceHanSerifSC-Light.otf")
plt.title('企鹅各个种类占比',fontproperties=zhfont1)
plt.show()

结果：

image-20230204164120934

每个岛屿有多少企鹅？

1
2
3
4
5
6
7
8
9
# 每个岛屿有多少企鹅
island_count=df['island'].value_counts()
print(island_count)
# 绘制饼图
plt.figure(figsize=(4,4))
labels=island_count.index
plt.pie(island_count,labels=labels,autopct='%.2f%%',explode=(0,0,0.1))
plt.title('各个岛屿企鹅数占比',fontproperties=zhfont1)
plt.show()

结果：

image-20230204164243087

每种类型企鹅各种体征属性的均值和分布

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#按种类找出数据
Adelie=df[df['species'].isin(['Adelie'])]
Gentoo=df[df['species'].isin(['Gentoo'])]
Chinstrap=df[df['species'].isin(['Chinstrap'])]

# 每种类型企鹅各种体征属性的均值和分布
Adelie_means=Adelie.mean()
Gentoo_means=Gentoo.mean()
Chinstrap_means=Chinstrap.mean()

#均值可视化
plt.figure(1)
# culmen_length_mm
plt.subplot(2,2,1)
plt.bar(['Adelie','Gentoo','Chinstrap'],
        [Adelie_means['culmen_length_mm'],
         Gentoo_means['culmen_length_mm'],
         Chinstrap_means['culmen_length_mm']]
        )
plt.xlabel('species')
plt.ylabel('culmen_length_mm')
# culmen_depth_mm
plt.subplot(2,2,2)
plt.bar(['Adelie','Gentoo','Chinstrap'],
        [Adelie_means['culmen_depth_mm'],
         Gentoo_means['culmen_depth_mm'],
         Chinstrap_means['culmen_depth_mm']]
        )
plt.xlabel('species')
plt.ylabel('culmen_depth_mm')
# flipper_length_mm
plt.subplot(2,2,3)
plt.bar(['Adelie','Gentoo','Chinstrap'],
        [Adelie_means['flipper_length_mm'],
         Gentoo_means['flipper_length_mm'],
         Chinstrap_means['flipper_length_mm']]
        )
plt.xlabel('species')
plt.ylabel('flipper_length_mm')
# body_mass_g
plt.subplot(2,2,4)
plt.bar(['Adelie','Gentoo','Chinstrap'],
        [Adelie_means['body_mass_g'],
         Gentoo_means['body_mass_g'],
         Chinstrap_means['body_mass_g']]
        )
plt.xlabel('species')
plt.ylabel('body_mass_g ')
plt.tight_layout() #自适应分布
plt.show()

结果

image-20230204164449475

不同种企鹅的分布

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
# culmen_length_mm分布
plt.figure()
plt.subplot(2,2,1)
species=[Adelie,Gentoo,Chinstrap]
labels=['Adelie','Gentoo','Chinstrap']
for i in range(3):
    sns.kdeplot(species[i]['culmen_length_mm'],shade=True,label=labels[i])
plt.xlabel('culmen_length_mm')
plt.ylabel('display')
# culmen_depth_mm
plt.subplot(2,2,2)
for i in range(3):
    sns.kdeplot(species[i]['culmen_depth_mm'],shade=True,label=labels[i])
plt.xlabel('culmen_depth_mm')
plt.ylabel('display')
# flipper_length_mm
plt.subplot(2,2,3)
for i in range(3):
    sns.kdeplot(species[i]['flipper_length_mm'],shade=True,label=labels[i])
plt.xlabel('flipper_length_mm')
plt.ylabel('display')
# body_mass_g
plt.subplot(2,2,4)
for i in range(3):
    sns.kdeplot(species[i]['body_mass_g'],shade=True,label=labels[i])
plt.xlabel('body_mass_g')
plt.ylabel('display')
plt.tight_layout() #自适应分布
plt.show()

结果：

image-20230204164609212

嘴峰长度和深度的关联

1
2
3
4
5
6
7
8
9
# 嘴峰长度和深度的关联
markers=['o','^','>']
plt.figure()
for i in range(3):
    plt.scatter(species[i]['culmen_length_mm'],species[i]['culmen_depth_mm'],label=labels[i],marker=markers[i])
plt.xlabel('culmen_length_mm')
plt.ylabel('culmen_depth_mm')
plt.legend()
plt.show()

image-20230204164705383

相似相异指标

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# 采用数字特征，计算皮尔逊系数
Adelie_corr = Adelie.corr()
Gentoo_corr=Gentoo.corr()
Chinstrap_corr=Chinstrap.corr()

corrs=[Adelie_corr,Gentoo_corr,Chinstrap_corr]
# 并绘制热图

for i in range(3):
    plt.figure()
    sns.heatmap(corrs[i], vmax=1, vmin=-1, center=0)
    plt.xticks(rotation=-15)
    plt.yticks(rotation=75)
    plt.title(labels[i]+'各特征相关性',fontproperties=zhfont1)
    plt.tight_layout() #自适应分布
    plt.show()

image-20230204164830893image-20230204164842868image-20230204164853560

数据散布矩阵

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# 数据散布图矩阵
features=['culmen_length_mm','culmen_depth_mm','flipper_length_mm','body_mass_g']
plt.figure(figsize=(25,25))
n=1
for i in range(4):
    for j in range(4):
        aix=plt.subplot(4,4,n)
        if i!=j:
            for k in range(3):
                plt.scatter(species[k][features[j]],species[k][features[i]],label=labels[k],marker=markers[k])
            plt.xlabel('culmen_length_mm')
            plt.ylabel('culmen_depth_mm')
            n=n+1
        else:
            plt.scatter(0,0,color='white')
            n=n+1
plt.tight_layout() #自适应分布
plt.show()
image-20230204165102779

不同种类企鹅的体重是否有显著性差异

1
2
3
4
5
6
7
8
9
10
11
# 不同种类企鹅的体重是否有显著性差异
body_mass_gs=[Adelie['body_mass_g'],Gentoo['body_mass_g'],Chinstrap['body_mass_g']]
# 绘制箱图
plt.figure()
plt.boxplot(body_mass_gs,labels=labels,showmeans=True)
plt.show()
# 绘制小提琴图
plt.figure()
plt.violinplot(body_mass_gs,showmeans=True)
plt.xticks(ticks = [1, 2,3], labels = labels, fontsize = 11)
plt.show()

image-20230204165141807image-20230204165150819

总结

EDA是为了发现各种问题，全面的了解数据，以直观的方式表达。可以从以下的角度入手：

  1. 特征均值、方差、分布
  2. 样本数据与特征的关系
  3. 特征与特征间的关系

实验四 分类方法-SVM算法

实验目的

  1. 掌握sklearn调用SVM
  2. 了解sklearn中SVM的参数和使用效果
  3. 了解sklearn的正则化方法，以及确认正则化系数的方法
  4. 实现线性SVM算法

实验场地与设备

实验室4074

实验方式

阅读教程与程序设计

实验设计

image-20230204180455883

实验内容

1.基本知识点总结

（1) 基本概念

支持向量机（support vector machines, SVM）是一种二分类模型，它的基本模型是定义在特征空间上的间隔最大的线性分类器；SVM还包括核技巧(kernel trick)，这使它成为实质上的非线性分类器。

对于支持向量机来说，数据点被视为p 维向量，而我们想知道是否可以用 (p-1)维超平面来分开这些点。这就是所谓的线性分类器。可能有许多超平面可以把数据分类。最佳超平面的一个合理选择是以最大间隔把两个类分开的超平面。因此，我们要选择能够让到每边最近的数据点的距离最大化的超平面。如果存在这样的超平面，则称为最大间隔超平面，而其定义的线性分类器被称为最大间隔分类器

（2) 分类

按照决策面类型分为：

  • 核函数为线性的SVM
  • 线性SVM分类器
  • 核函数为径向基核函数的SVM分类器
  • 核函数为多项式的SVM分类器

（3）算法

SVM的学习策略就是间隔最大化，可形式化为一个求解凸二次规划的问题，也等价于正则化的合页损失函数的最小化问题。SVM的的学习算法就是求解凸二次规划的最优化算法。

举出线性支持向量机算法如下：

输入：训练数据集\(T={(x_1,y_1),(x_2,y_2),…,(x_N,y_N)}\)

其中\(x_i∈R^n,y_i∈{+1,-1},i=1,2,...,N\)

输出：分离超平面和分类决策函数

（1）选择惩罚参数 \(C>0\)，构造并求解凸二次规划问题 \[ min 1/2 \sum_{i=1}^N\sum_{j=1}^N{α_i α_j y_i y_j (x_i·x_j)}-\sum_{i=1}^N{α_i} \]

\[ s.t.∑_{i=1}^N{α_i y_i}=0,0≤α_i≤C,i=1,2,...,N \]

得到最优解 \(α^*=(α_1^*,α_2^*,...,α_N^*)\)

（2）计算\(w^*=∑_{i=1}^N{α_i^* y_i x_i }\)

选择 \(\alpha^*\) 的一个分量 \(\alpha^*_j\)满足条件 \(0<\alpha^*_j<C\) ，计算\(b^*=y_i-∑_{i=1}^N{α_i^* y_i }\)

（3）求分离超平面\(w^*·x+b^*=0\)

分类决策函数：\(f(x)=sign(w^*·x+b^*)\)

（4） 模型特点

SVM的潜在缺点包括以下方面：

  • 需要对输入数据进行完全标记

  • 未校准类成员概率

  • SVM仅直接适用于两类任务。因此，必须应用将多类任务减少到几个二元问题的算法

2.线性SVM实现

SMO算法

SMO表示序列最小优化 ( Sequential Minimal Optimization )。Platt的SMO算法是将大优化问题分解为多个小优化问题来求解的。这些小优化问题往往很容易求解, 并且对它们进行顺序求解的结果与将它们作为整体来求解的结果是完全一致的。在结果完全相同的同时, SMO算法的求解时间短很多。

简化版SMO

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
smoSimple.py
input：maxIter,X,y,toler,C
output：b，alphas

 alphas = np.mat(np.zeros((m,1)))#初始化alpha参数，设置为0
 for iterSmo in range(maxIter):
 	for i in range(np.shape(X)[0]):
 		计算Ei
 		if ((y[i]*Ei < -toler) and (alphas[i] < C)) or ((y[i]*Ei > -toler) and (alphas[i] > 0)):
            随机选择一个向量
            j = selectJrand(i,m)
            进行优化
            alphaPairsChanged += 1
            #打印统计信息
            print("第%d次迭代 样本:%d, alpha优化次数:%d" % (iterSmo,i,alphaPairsChanged))
 		if (alphaPairsChanged == 0): iterSmo += 1
        else:
            iterSmo = 0
        print("迭代次数：%d" % iterSmo)
    return b,alphas

Python实现：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
import numpy as np
import random

def selectJrand(i,m):
    j=i
    while(j==i): #选择一个不等于i的j
        j = int(random.uniform(0,m))
    return j

def clipAlpha(aj,H,L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

#简化版SMO算法
def smoSimple(dataMatIn,classLabels,C,toler,maxIter):
    dataMatrix = np.mat(dataMatIn); labelMat = np.mat(classLabels).transpose()
    b = 0; m,n = np.shape(dataMatrix)
    alphas = np.mat(np.zeros((m,1)))#初始化alpha参数，设置为0
    iterSmo = 0 #初始化迭代次数
    while(iterSmo < maxIter):
        alphaPairsChanged = 0#用于记录alpha是否已经进行优化
        #步骤1. 计算误差Ei
        for i in range(m):
            fXi = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
            Ei = fXi - float(labelMat[i])
            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > -toler) and (alphas[i] > 0)):
                j = selectJrand(i,m)
                #步骤1. 计算误差Ej
                fXj = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b
                Ej = fXj - float(labelMat[j])
                #保存更新前的alpha值，使用浅拷贝
                alphaIold = alphas[i].copy()
                alphaJold = alphas[j].copy()
                #步骤2：计算上界H和下界L
                if (labelMat[i] != labelMat[j]):
                    L = max(0,alphas[j] - alphas[i])
                    H = min(C,C + alphas[j] - alphas[i])
                else:
                    L = max(0,alphas[j] + alphas[i] - C)
                    H = min(C,alphas[j] + alphas[i])
                if L==H: print("L==H"); continue
                #步骤3：计算eta
                eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
                if eta >=0 : print("eta>=0");continue
                #步骤4：更新alpha_j
                alphas[j] -= labelMat[j]*(Ei-Ej)/eta
                #步骤5：修剪alpha_j
                alphas[j] = clipAlpha(alphas[j],H,L)
                if (abs(alphas[j] - alphaJold) < 0.00001) : print("j not moving enough") ; continue
                #步骤6：更新alpha_i
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])
                #步骤7：更新b_1和b_2
                b1 = b - Ei - labelMat[i]*(alphas[i] - alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T \
                - labelMat[j]*(alphas[j] - alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T
                b2 = b - Ej - labelMat[i]*(alphas[i] - alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T \
                - labelMat[j]*(alphas[j] - alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T
                #步骤8：根据b_1和b_2更新b
                if (0 < alphas[i]) and (C > alphas[i]) :
                    b = b1
                elif (0 < alphas[j]) and (C > alphas[j]):
                    b = b2
                else:
                    b = (b1 + b2)/2.0
                #统计优化次数
                #如果程序执行到for循环的最后一行都不执行continue语句，那么就已经成功地改变了一对alpha，同时可以增加alphaPairsChanged的值
                alphaPairsChanged += 1
                #打印统计信息
                print("第%d次迭代 样本:%d, alpha优化次数:%d" % (iterSmo,i,alphaPairsChanged))
        #更新迭代次数
        #在for循环之外，需要检查alpha值是否做了更新，如果有更新则将iterSmo设为0后继续运行程序。只有在所有数据集上遍历maxIter次，且不再发生任何alpha修改之后，程序才会停止并退出while循环
        if (alphaPairsChanged == 0): iterSmo += 1
        else:
            iterSmo = 0
        print("迭代次数：%d" % iterSmo)
    return b,alphas

3.支持向量机核方法对比

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
import numpy as np
import matplotlib.pyplot as plt
from sklearn import svm, datasets


def make_meshgrid(x, y, h=.02):
    """创建要绘制的点网格

    参数
    ----------
    x: 创建网格x轴所需要的数据
    y: 创建网格y轴所需要的数据
    h: 网格大小的可选大小，可选填

    返回
    -------
    xx, yy : n维数组
    """
    x_min, x_max = x.min() - 1, x.max() + 1
    y_min, y_max = y.min() - 1, y.max() + 1
    xx, yy = np.meshgrid(np.arange(x_min, x_max, h),
                         np.arange(y_min, y_max, h))
    return xx, yy


def plot_contours(ax, clf, xx, yy, **params):
    """绘制分类器的决策边界。

    参数
    ----------
    ax: matplotlib子图对象
    clf: 一个分类器
    xx: 网状网格meshgrid的n维数组
    yy: 网状网格meshgrid的n维数组
    params: 传递给contourf的参数字典，可选填
    """
    Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])
    Z = Z.reshape(xx.shape)
    out = ax.contourf(xx, yy, Z, **params)
    return out


# 导入数据以便后续使用
iris = datasets.load_iris()
# 采用前两个特征
X = iris.data[:, :2]
y = iris.target

# 创建一个SVM实例并拟合数据。由于要绘制支持向量，因此我们不缩放数据
C = 1.0  # SVM正则化参数
models = (svm.SVC(kernel='linear', C=C),
          svm.LinearSVC(C=C, max_iter=10000),
          svm.SVC(kernel='rbf', gamma=0.7, C=C),
          svm.SVC(kernel='poly', degree=3, gamma='auto', C=C))
models = (clf.fit(X, y) for clf in models)

# 为图像设置标题
titles = ('SVC with linear kernel',
          'LinearSVC (linear kernel)',
          'SVC with RBF kernel',
          'SVC with polynomial (degree 3) kernel')

# 设置一个2x2结构的画布
fig, sub = plt.subplots(2, 2)
plt.subplots_adjust(wspace=0.4, hspace=0.4)

X0, X1 = X[:, 0], X[:, 1]
xx, yy = make_meshgrid(X0, X1)

for clf, title, ax in zip(models, titles, sub.flatten()):
    plot_contours(ax, clf, xx, yy,
                  cmap=plt.cm.coolwarm, alpha=0.8)
    ax.scatter(X0, X1, c=y, cmap=plt.cm.coolwarm, s=20, edgecolors='k')
    ax.set_xlim(xx.min(), xx.max())
    ax.set_ylim(yy.min(), yy.max())
    ax.set_xlabel('Sepal length')
    ax.set_ylabel('Sepal width')
    ax.set_xticks(())
    ax.set_yticks(())
    ax.set_title(title)

plt.show()
image-20230204171813486

4.通过SVM对鸢尾花数据集做预处理、训练、预测并进行准确率估计。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
from sklearn import svm
import numpy as np
from sklearn import model_selection
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib import colors
from sklearn import datasets

#读入数据
iris = datasets.load_iris()
#定义训练集和测试集
X= iris.data
y=iris.target
x=X[:,0:2]
#在 X中取前两列作为特征（为了后期的可视化画图）
x_train,x_test,y_train,y_test=model_selection.train_test_split(x,y,random_state=0,test_size=0.3)
# 用train_test_split将数据随机分为训练集和测试集，测试集占总数据的30%（test_size=0.3),random_state是随机数种子

#搭建模型，训练SVM分类器
classifier=svm.SVC(kernel='rbf',gamma=0.1,decision_function_shape='ovo',C=0.8)
#开始训练
classifier.fit(x_train,y_train.ravel())

def show_accuracy(y_hat,y_train,str):
    pass

#（4）计算svm分类器的准确率
print("SVM-输出训练集的准确率为：",classifier.score(x_train,y_train))
y_hat=classifier.predict(x_train)
show_accuracy(y_hat,y_train,'训练集')
print("SVM-输出测试集的准确率为：",classifier.score(x_test,y_test))
y_hat=classifier.predict(x_test)
show_accuracy(y_hat,y_test,'测试集')


# 查看决策函数，通过decision_function()实现。decision_function中每一列的值代表距离各类别的距离。
print('\npredict:\n', classifier.predict(x_train))


#绘制图像
x1_min, x1_max = x[:, 0].min(), x[:, 0].max()
x2_min, x2_max = x[:, 1].min(), x[:, 1].max()
x1, x2 = np.mgrid[x1_min:x1_max:200j, x2_min:x2_max:200j]

grid_test = np.stack((x1.flat, x2.flat), axis=1)
# 测试点，再通过stack()函数，axis=1，生成测试点

print("grid_test = \n", grid_test)
# 预测分类值
grid_hat = classifier.predict(grid_test)

print("grid_hat = \n", grid_hat)

grid_hat = grid_hat.reshape(x1.shape)
#绘制
cm_light = mpl.colors.ListedColormap(['#A0FFA0', '#FFA0A0', '#A0A0FF'])
cm_dark = mpl.colors.ListedColormap(['g', 'r', 'b'])

alpha=0.5

plt.pcolormesh(x1, x2, grid_hat, cmap=cm_light) # 预测值的显示
# plt.scatter(x[:, 0], x[:, 1], c=y, edgecolors='k', s=50, cmap=cm_dark)  # 样本
plt.plot(x[:, 0], x[:, 1], 'o', alpha=alpha, color='blue', markeredgecolor='k')
plt.scatter(x_test[:, 0], x_test[:, 1], s=120, facecolors='none', zorder=10)  # 圈中测试集样本
plt.xlabel(u'花萼长度', fontsize=13)
plt.ylabel(u'花萼宽度', fontsize=13)
plt.xlim(x1_min, x1_max)
plt.ylim(x2_min, x2_max)
plt.title(u'鸢尾花SVM二特征分类', fontsize=15)
# plt.grid()
plt.show()
image-20230204172358319
image-20230204172435520

5.正则化项选择

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
from sklearn.svm import SVC
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import permutation_test_score
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris

iris = load_iris()
X = iris.data
y = iris.target

clf = SVC(kernel="linear", random_state=7)
cv = StratifiedKFold(2, shuffle=True, random_state=0)

score_iris, perm_scores_iris, pvalue_iris = permutation_test_score(
    clf, X, y, scoring="accuracy", cv=cv, n_permutations=1000
)


print(score_iris)

6.SVR

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
from sklearn.svm import SVR
import numpy as np
import scipy.io as sio
import matplotlib.pyplot as plt

data=sio.loadmat("NIRcorn.mat")
cornwavelength=data["cornwavelength"]
x=[]
# 输出观察发现，是一个ndarray组成的list所以只取第一个元素
for i in range(len(cornwavelength)):
    x.append(cornwavelength[i][0])
x=np.array(x).reshape(-1, 1)
# 获取吸光率
cornspect = data['cornspect']
# 随机取出1组透光率
rd=np.random.randint(0,80,1)
# 存储透光率数据
y = []
for i in range(0, 700):
    y.append(cornspect[rd[0], i])
y=np.array(y).ravel()
#径向基函数拟合模型
svr_rbf = SVR(kernel="rbf", C=100, gamma=0.1, epsilon=0.1)
# 输出结果
fig=plt.figure(1)
plt.plot(x, svr_rbf.fit(x, y).predict(x))
plt.scatter(x,y,s=2,c='green')
plt.show()
#多项式核拟合模型
svr_poly = SVR(kernel="rbf", C=50, gamma='auto', epsilon=0.1)
# 输出结果
fig=plt.figure(2)
plt.plot(x, svr_poly.fit(x, y).predict(x))
plt.scatter(x,y,s=2,c='red')
plt.show()

image-20230204181207801

总结

SVM作为一个经典的模型，他综合了代数、几何、优化、泛函等一系列数学知识，是一个简洁而又漂亮的模型。这种模型对于我继续学习具有很大帮助。

实验五 关联分析算法应用

实验目的

  1. 对给定数据集进行整理，利用apriori， FP-Growth算法对给定数据集进行分类；
  2. 理解不同分类的算法的原理，掌握apriori， FP-Growth 算法的应用；
  3. 掌握apriori算法原理；
  4. 掌握FP-Growth算法原理；
  5. 调用sklearn中的数据集进行相关算法的实验

实验场地与设备

实验室4074

实验方式

阅读教程与程序设计

实验设计

image-20230204175935517

实验内容

1.知识总结

（1）基本概念

支持度：关联规则A->B的支持度support=P(AB)，指的是事件A和事件B同时发生的概率

置信度：置信度confidence=P(B|A)=P(AB)/P(A),指的是发生事件A的基础上发生事件B的概率。

k项集：如果事件A中包含k个元素，那么称这个事件A为k项集，并且事件A满足最小支持度阈值的事件称为频繁k项集

（2）Apriori算法

Apriori算法过程分为两个步骤：

第一步通过迭代，检索出事务数据库中的所有频繁项集，即支持度不低于用户设定的阈值的项集；

第二步利用频繁项集构造出满足用户最小信任度的规则。

（3） FP-Growth算法

FP Tree算法包括三步：

  • 扫描数据，得到所有频繁一项集的的计数。然后删除支持度低于阈值的项，将1项频繁集放入项头表，并按照支持度降序排列。扫描数据，将读到的原始数据剔除非频繁1项集，并按照支持度降序排列。读入排序后的数据集，插入FP树，插入时按照排序后的顺序，插入FP树中，排序靠前的节点是祖先节点，而靠后的是子孙节点。如果有共用的祖先，则对应的公用祖先节点计数加1。插入后，如果有新节点出现，则项头表对应的节点会通过节点链表链接上新节点。直到所有的数据都插入到FP树后，FP树的建立完成。

  • 从项头表的底部项依次向上找到项头表项对应的条件模式基。从条件模式基递归挖掘得到项头表项项的频繁项集。

  • 如果不限制频繁项集的项数，则返回步骤4所有的频繁项集，否则只返回满足项数要求的频繁项集。

  1. 两大算法比较

经典的关联规则挖掘算法包括Apriori算法和FP-growth算法。apriori算法多次扫描交易数据库，每次利用候选频繁集产生频繁集；而FP-growth则利用树形结构，无需产生候选频繁集而是直接得到频繁集，大大减少扫描交易数据库的次数，从而提高了算法的效率，FP Tree算法改进了Apriori算法的I/O瓶颈，巧妙的利用了树结构，这让我们想起了BIRCH聚类，BIRCH聚类也是巧妙的利用了树结构来提高算法运行速度。利用内存数据结构以空间换时间是常用的提高算法运行时间瓶颈的办法。但是apriori的算法扩展性较好，可以用于并行计算等领域。

2.Apriori算法进行关联分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
import sys


def load_data_set():
    data_set = [
        ['beer', 'baby diapers', 'shorts']
        , ['baby diapers', 'shorts']
        , ['baby diapers', 'milk']
        , ['beer', 'baby diapers', 'shorts']
        , ['beer', 'milk']
        , ['baby diapers', 'milk']
        , ['beer', 'milk']
        , ['beer', 'baby diapers', 'milk', 'shorts']
        , ['beer', 'baby diapers', 'milk']
    ]
    return data_set


def create_C1(data_set):
    C1 = set()
    for t in data_set:
        for item in t:
            item_set = frozenset([item])
            C1.add(item_set)
    return C1


def is_apriori(Ck_item, Lksub1):
    for item in Ck_item:
        sub_Ck = Ck_item - frozenset([item])
        if sub_Ck not in Lksub1:
            return False
    return True


def create_Ck(Lksub1, k):
    Ck = set()
    len_Lksub1 = len(Lksub1)
    list_Lksub1 = list(Lksub1)
    for i in range(len_Lksub1):
        for j in range(1, len_Lksub1):
            l1 = list(list_Lksub1[i])
            l2 = list(list_Lksub1[j])
            l1.sort()
            l2.sort()
            if l1[0:k-2] == l2[0:k-2]:
                Ck_item = list_Lksub1[i] | list_Lksub1[j]
                if is_apriori(Ck_item, Lksub1):
                    Ck.add(Ck_item)
    return Ck


def generate_Lk_by_Ck(data_set, Ck, min_support, support_data):
    Lk = set()
    item_count = {}
    for t in data_set:
        for item in Ck:
            if item.issubset(t):
                if item not in item_count:
                    item_count[item] = 1
                else:
                    item_count[item] += 1
    t_num = float(len(data_set))
    for item in item_count:
        if (item_count[item] / t_num) >= min_support:
            Lk.add(item)
            support_data[item] = item_count[item] / t_num
    return Lk


def generate_L(data_set, k, min_support):
    support_data = {}
    C1 = create_C1(data_set)
    L1 = generate_Lk_by_Ck(data_set, C1, min_support, support_data)
    Lksub1 = L1.copy()
    L = []
    L.append(Lksub1)
    for i in range(2, k+1):
        Ci = create_Ck(Lksub1, i)
        Li = generate_Lk_by_Ck(data_set, Ci, min_support, support_data)
        Lksub1 = Li.copy()
        L.append(Lksub1)
    return L, support_data


def generate_big_rules(L, support_data, min_conf):
    big_rule_list = []
    sub_set_list = []
    for i in range(0, len(L)):
        for freq_set in L[i]:
            for sub_set in sub_set_list:
                if sub_set.issubset(freq_set):
                    conf = support_data[freq_set] / support_data[freq_set - sub_set]
                    big_rule = (freq_set - sub_set, sub_set, conf)
                    if conf >= min_conf and big_rule not in big_rule_list:
                        big_rule_list.append(big_rule)
            sub_set_list.append(freq_set)
    return big_rule_list


if __name__ == "__main__":
    """
    Test
    """
    data_set = load_data_set()
    L, support_data = generate_L(data_set, k=3, min_support=0.2)
    big_rules_list = generate_big_rules(L, support_data, min_conf=0.7)
    for Lk in L:
        print("="*50)
        print("frequent " + str(len(list(Lk)[0])) + "-itemsets\t\tsupport")
        print("="*50)
        for freq_set in Lk:
            print(freq_set, support_data[freq_set])
    print
    print("Big Rules")
    for item in big_rules_list:
        print(item[0], "=>", item[1], "conf: ", item[2])

image-20230204173031686

3.FP-Growth算法进行关联分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
class treeNode:
    def __init__(self, nameValue, numOccur, parentNode):
        self.name = nameValue
        self.count = numOccur
        self.nodeLink = None #nodeLink 变量用于链接相似的元素项
        self.parent = parentNode #指向当前节点的父节点
        self.children = {} #空字典，存放节点的子节点

    def inc(self, numOccur):#计数加1
        self.count += numOccur

#将树以文本形式显示
    def disp(self, ind=1):
        print ('  ' * ind, self.name, ' ', self.count)
        for child in self.children.values():
            child.disp(ind + 1)

#构建FP-tree
def createTree(dataSet, minSup=1):
    headerTable = {}
    for trans in dataSet:  #第一次遍历：统计各个数据的频繁度
        for item in trans:
            headerTable[item] = headerTable.get(item, 0) + dataSet[trans]
            #用头指针表统计各个类别的出现的次数，计算频繁量：头指针表[类别]=出现次数
    for k in list(headerTable):  #删除未达到最小频繁度的数据
        if headerTable[k] < minSup:
            del (headerTable[k])
    freqItemSet = set(headerTable.keys())#保存达到要求的数据
    # print ('freqItemSet: ',freqItemSet)
    if len(freqItemSet) == 0:
        return None, None  #若达到要求的数目为0
    for k in headerTable: #遍历头指针表
        headerTable[k] = [headerTable[k], None]  #保存计数值及指向每种类型第一个元素项的指针
    # print ('headerTable: ',headerTable)
    retTree = treeNode('Null Set', 1, None)  #初始化tree
    for tranSet, count in dataSet.items():  # 第二次遍历：
        localD = {}
        for item in tranSet:  # put transaction items in order
            if item in freqItemSet:#只对频繁项集进行排序
                localD[item] = headerTable[item][0]

        #使用排序后的频率项集对树进行填充
        if len(localD) > 0:
            orderedItems = [v[0] for v in sorted(localD.items(), key=lambda p: p[1], reverse=True)]
            updateTree(orderedItems, retTree, headerTable, count)  # populate tree with ordered freq itemset
    return retTree, headerTable  #返回树和头指针表


def updateTree(items, inTree, headerTable, count):
    if items[0] in inTree.children:  # 首先检查是否存在该节点
        inTree.children[items[0]].inc(count)  # 存在则计数增加
    else:  # 不存在则将新建该节点
        inTree.children[items[0]] = treeNode(items[0], count, inTree)#创建一个新节点
        if headerTable[items[0]][1] == None:  # 若原来不存在该类别，更新头指针列表
            headerTable[items[0]][1] = inTree.children[items[0]]#更新指向
        else:#更新指向
            updateHeader(headerTable[items[0]][1], inTree.children[items[0]])
    if len(items) > 1:  #仍有未分配完的树，迭代
        updateTree(items[1::], inTree.children[items[0]], headerTable, count)

#节点链接指向树中该元素项的每一个实例。
# 从头指针表的 nodeLink 开始,一直沿着nodeLink直到到达链表末尾
def updateHeader(nodeToTest, targetNode):
    while (nodeToTest.nodeLink != None):
        nodeToTest = nodeToTest.nodeLink
    nodeToTest.nodeLink = targetNode

# 将读取的文件的一行数据转换成数组
def cast_line_to_array(line):
    lines = line.split()
    line_data = []
    for i in lines:
        line_data.append(int(i))
    return line_data
pass

    # 从文件中读取数据
def load_data():
    data = [
    ['beer', 'baby diapers', 'shorts']
    , ['baby diapers', 'shorts']
    , ['baby diapers', 'milk']
    , ['beer', 'baby diapers', 'shorts']
    , ['beer', 'milk']
    , ['baby diapers', 'milk']
    , ['beer', 'milk']
    , ['beer', 'baby diapers', 'milk', 'shorts']
    , ['beer', 'baby diapers', 'milk']
    ]
#     file = open(filepath)
#     line = file.readline()
#     data.append(cast_line_to_array(line))
#     while line:
#         line = file.readline()
#         data.append(cast_line_to_array(line))
    return data
pass

def loadSimpDat():
    simpDat = [['r', 'z', 'h', 'j', 'p'],
               ['z', 'y', 'x', 'w', 'v', 'u', 't', 's'],
               ['z'],
               ['r', 'x', 'n', 'o', 's'],
               ['y', 'r', 'x', 'z', 'q', 't', 'p'],
               ['y', 'z', 'x', 'e', 'q', 's', 't', 'm']]
    return simpDat
#createInitSet() 用于实现上述从列表到字典的类型转换过程
def createInitSet(dataSet):
    retDict = {}
    for trans in dataSet:
        retDict[frozenset(trans)] = 1
    return retDict

#从FP树中发现频繁项集
def ascendTree(leafNode, prefixPath):  #递归上溯整棵树
    if leafNode.parent != None:
        prefixPath.append(leafNode.name)
        ascendTree(leafNode.parent, prefixPath)


def findPrefixPath(basePat, treeNode):  #参数：指针，节点；
    condPats = {}
    while treeNode != None:
        prefixPath = []
        ascendTree(treeNode, prefixPath)#寻找当前非空节点的前缀
        if len(prefixPath) > 1:
            condPats[frozenset(prefixPath[1:])] = treeNode.count #将条件模式基添加到字典中
        treeNode = treeNode.nodeLink
    return condPats

#递归查找频繁项集
def mineTree(inTree, headerTable, minSup, preFix, freqItemList):
    # 头指针表中的元素项按照频繁度排序,从小到大
    bigL = [v[0] for v in sorted(headerTable.items(), key=lambda p: str(p[1]))]#python3修改
    for basePat in bigL:  #从底层开始
        #加入频繁项列表
        newFreqSet = preFix.copy()
        newFreqSet.add(basePat)
        #print ('finalFrequent Item: ',newFreqSet)
        freqItemList.append(newFreqSet)
        #递归调用函数来创建基
        condPattBases = findPrefixPath(basePat, headerTable[basePat][1])
        #print ('condPattBases :',basePat, condPattBases)

        #2. 构建条件模式Tree
        myCondTree, myHead = createTree(condPattBases, minSup)
        #将创建的条件基作为新的数据集添加到fp-tree
        #print ('head from conditional tree: ', myHead)
        if myHead != None: #3. 递归
            #print('conditional tree for: ',newFreqSet)
            #myCondTree.disp(1)
            mineTree(myCondTree, myHead, minSup, newFreqSet, freqItemList)


if __name__ == "__main__":
    #simpDat = loadSimpDat()
    simpDat = load_data()
    #print(simpDat)
    initSet = createInitSet(simpDat)
    #print(initSet)
    myFPtree, myHeaderTab = createTree(initSet, 3)
    #myFPtree.disp()
    #print(findPrefixPath('t', myHeaderTab['t'][1]))
    freqItems = []
    mineTree(myFPtree, myHeaderTab, 3, set([]), freqItems)
    print('将树以文本形式展示:')
    print(myFPtree.disp())
    print(freqItems)
    print()

image-20230204173146239

总结

这两个算法Aprior比较容易理解，fp-grow的方式虽然有数据结构的基础，在python上实现依旧比较有难度。

实验六 聚类分析

实验目的

  1. 实现k-means算法
  2. 通过sklearn库调用kmeans算法以及改进算法；

实验场地与设备

实验室4074

实验方式

阅读教程与程序设计

实验设计

image-20230204175718032

实验内容

1.Kmeans实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sklearn.datasets import load_wine
plt.style.use('seaborn')
import Select_K


class Kmeans:
    def __init__(self, k=2, tolerance=0.01):
        self.k = k
        self.tol = tolerance
        self.features_count = -1
        self.classifications = None
        self.centroids = None

    def fit(self, data):
        """
        :param data: numpy数组，约定shape为：(数据数量，数据维度)
        :type data: numpy.ndarray
        """
        self.features_count = data.shape[1]
        # 初始化聚类中心（维度：k个 * features种数）
        self.centroids = np.zeros([self.k, data.shape[1]])
        # self.centroids = np.array(gene).reshape(3, 13)
        for i in range(self.k):
            self.centroids[i] = data[i]

        while True:
            # 清空聚类列表
            self.classifications = [[] for i in range(self.k)]
            # 对每个点与聚类中心进行距离计算
            for feature_set in data:
                # 预测分类
                classification = self.predict(feature_set)
                # 加入类别
                self.classifications[classification].append(feature_set)

            # 记录前一次的结果
            prev_centroids = np.ndarray.copy(self.centroids)

            # 更新中心
            for classification in range(self.k):
                self.centroids[classification] = np.average(self.classifications[classification], axis=0)

            # 检测相邻两次中心的变化情况
            for c in range(self.k):
                if np.linalg.norm(prev_centroids[c] - self.centroids[c]) > self.tol:
                    break

            # 如果都满足条件（上面循环没break），则返回
            else:
                return

    def predict(self, data):
        # 距离
        distances = np.linalg.norm(data - self.centroids, axis=1)
        # 最小距离索引
        return distances.argmin()

    def showkmeans(self,labels):
        fig = plt.figure()
        plt.style.use('seaborn')
        colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'w']
        for i in range(self.k):
            color = colors[i % len(colors)]
            for feature_set in self.classifications[i]:
                plt.scatter(feature_set[0], feature_set[1],color=color)
        for centroid in self.centroids:
            plt.scatter(
                        centroid[0],
                        centroid[1],
                        marker="^",
                        color="orange",
                        s=70,
                        edgecolors='k',
                        linewidths=1.5
                        )
        plt.xlabel(labels[0])
        plt.ylabel(labels[1])
        plt.tight_layout()
        plt.show()

    def showkmeansByPCA(self,pca):
        fig = plt.figure()
        colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'w']
        plt.style.use('seaborn')
        for i in range(self.k):
            color = colors[i % len(colors)]
            for feature_set in pca.transform(self.classifications[i]):
                plt.scatter(feature_set[0], feature_set[1],color=color)
        for centroid in pca.transform(self.centroids):
            plt.scatter(
                        centroid[0],
                        centroid[1],
                        marker="^",
                        color="orange",
                        s=70,
                        edgecolors='k',
                        linewidths=1.5
                        )
        plt.xlabel('Principle Componet1')
        plt.ylabel('Principle Componet2')
        plt.tight_layout()
        plt.show()

    def showkmeansByPCA3D(self,pca):
        fig = plt.figure()
        ax = plt.axes(projection='3d')
        colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'w']
        for i in range(self.k):
            color = colors[i % len(colors)]
            for feature_set in pca.transform(self.classifications[i]):
                ax.scatter3D(feature_set[0], feature_set[1]
                            ,feature_set[2],color=color)
        for centroid in pca.transform(self.centroids):
            ax.scatter3D(
                        centroid[0],
                        centroid[1],
                        centroid[2],
                        marker="^",
                        color="orange",
                        s=70,
                        edgecolors='k',
                        linewidths=1.5
                        )
        ax.set_xlabel('Principle Componet1')
        ax.set_ylabel('Principle Componet2')
        ax.set_zlabel('Principle Componet3')
        plt.tight_layout()
        plt.show()

2.Kmeans 的不适用情况

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
import numpy as np
import matplotlib.pyplot as plt

from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs

plt.figure(figsize=(12, 12))

n_samples = 1500
random_state = 170
X, y = make_blobs(n_samples=n_samples, random_state=random_state)

# Incorrect number of clusters
y_pred = KMeans(n_clusters=2, random_state=random_state).fit_predict(X)

plt.subplot(221)
plt.scatter(X[:, 0], X[:, 1], c=y_pred)
plt.title("Incorrect Number of Blobs")

# Anisotropicly distributed data
# 这块本质上就是对坐标进行变换，我们画出的坐标轴是垂直的，
# 但是实际上坐标轴应该是一个平行四边形
# transformation实际上就是对坐标轴进行旋转，通过反三角函数就可以控制
transformation = [[0.2, 0],
                  [0.4, 0.4]]
X_aniso = np.dot(X, transformation)
y_pred = KMeans(n_clusters=3, random_state=random_state).fit_predict(X_aniso)

plt.subplot(222)
plt.scatter(X_aniso[:, 0], X_aniso[:, 1], c=y_pred)
plt.title("Anisotropicly Distributed Blobs")

# Different variance
X_varied, y_varied = make_blobs(
    n_samples=n_samples, cluster_std=[1.0, 2.5, 0.5], random_state=random_state
)
y_pred = KMeans(n_clusters=3, random_state=random_state).fit_predict(X_varied)

plt.subplot(223)
plt.scatter(X_varied[:, 0], X_varied[:, 1], c=y_pred)
plt.title("Unequal Variance")

# Unevenly sized blobs
X_filtered = np.vstack((X[y == 0][:500], X[y == 1][:100], X[y == 2][:10]))
y_pred = KMeans(n_clusters=3, random_state=random_state).fit_predict(X_filtered)

plt.subplot(224)
plt.scatter(X_filtered[:, 0], X_filtered[:, 1], c=y_pred)
plt.title("Unevenly Sized Blobs")

plt.show()
image-20230204175246141

3.kmeans++ 处理图像数据集

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
"""
=========== ========================================================
Shorthand    full name
=========== ========================================================
homo         homogeneity score
compl        completeness score
v-meas       V measure
ARI          adjusted Rand index
AMI          adjusted mutual information
silhouette   silhouette coefficient
=========== ========================================================
"""

# %%
# Load the dataset
# ----------------
#
# We will start by loading the `digits` dataset. This dataset contains
# handwritten digits from 0 to 9. In the context of clustering, one would like
# to group images such that the handwritten digits on the image are the same.

import numpy as np
from sklearn.datasets import load_digits

data, labels = load_digits(return_X_y=True)
(n_samples, n_features), n_digits = data.shape, np.unique(labels).size

print(f"# digits: {n_digits}; # samples: {n_samples}; # features {n_features}")

# %%
# Define our evaluation benchmark
# -------------------------------
#
# We will first our evaluation benchmark. During this benchmark, we intend to
# compare different initialization methods for KMeans. Our benchmark will:
#
# * create a pipeline which will scale the data using a
#   :class:`~sklearn.preprocessing.StandardScaler`;
# * train and time the pipeline fitting;
# * measure the performance of the clustering obtained via different metrics.
from time import time
from sklearn import metrics
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler


def bench_k_means(kmeans, name, data, labels):
    """Benchmark to evaluate the KMeans initialization methods.

    Parameters
    ----------
    kmeans : KMeans instance
        A :class:`~sklearn.cluster.KMeans` instance with the initialization
        already set.
    name : str
        Name given to the strategy. It will be used to show the results in a
        table.
    data : ndarray of shape (n_samples, n_features)
        The data to cluster.
    labels : ndarray of shape (n_samples,)
        The labels used to compute the clustering metrics which requires some
        supervision.
    """
    t0 = time()
    estimator = make_pipeline(StandardScaler(), kmeans).fit(data)
    fit_time = time() - t0
    # 这里的源码是有问题的，estimator作为pipelinem不支持使用index，
    # 我通过debug找到了要的变量
    # 这一点我觉得是可以理解的，从流的角度上看pipeline使用index就相当于截取一段水流
    results = [name, fit_time, estimator.named_steps['kmeans'].inertia_]

    # Define the metrics which require only the true labels and estimator
    # labels
    clustering_metrics = [
        metrics.homogeneity_score,
        metrics.completeness_score,
        metrics.v_measure_score,
        metrics.adjusted_rand_score,
        metrics.adjusted_mutual_info_score,
    ]
    results += [m(labels, estimator.named_steps['kmeans'].labels_) for m in clustering_metrics]

    # The silhouette score requires the full dataset
    results += [
        metrics.silhouette_score(
            data,
            estimator.named_steps['kmeans'].labels_,
            metric="euclidean",
            sample_size=300,
        )
    ]

    # Show the results
    formatter_result = (
        "{:9s}\t{:.3f}s\t{:.0f}\t{:.3f}\t{:.3f}\t{:.3f}\t{:.3f}\t{:.3f}\t{:.3f}"
    )
    print(formatter_result.format(*results))


# %%
# Run the benchmark
# -----------------
#
# We will compare three approaches:
#
# * an initialization using `kmeans++`. This method is stochastic and we will
#   run the initialization 4 times;
# * a random initialization. This method is stochastic as well and we will run
#   the initialization 4 times;
# * an initialization based on a :class:`~sklearn.decomposition.PCA`
#   projection. Indeed, we will use the components of the
#   :class:`~sklearn.decomposition.PCA` to initialize KMeans. This method is
#   deterministic and a single initialization suffice.
from sklearn.cluster import KMeans
from sklearn.decomposition import PCA

print(82 * "_")
print("init\t\ttime\tinertia\thomo\tcompl\tv-meas\tARI\tAMI\tsilhouette")

kmeans = KMeans(init="k-means++", n_clusters=n_digits, n_init=4, random_state=0)
bench_k_means(kmeans=kmeans, name="k-means++", data=data, labels=labels)

kmeans = KMeans(init="random", n_clusters=n_digits, n_init=4, random_state=0)
bench_k_means(kmeans=kmeans, name="random", data=data, labels=labels)

pca = PCA(n_components=n_digits).fit(data)
kmeans = KMeans(init=pca.components_, n_clusters=n_digits, n_init=1)
bench_k_means(kmeans=kmeans, name="PCA-based", data=data, labels=labels)

print(82 * "_")

# %%
# Visualize the results on PCA-reduced data
# -----------------------------------------
#
# :class:`~sklearn.decomposition.PCA` allows to project the data from the
# original 64-dimensional space into a lower dimensional space. Subsequently,
# we can use :class:`~sklearn.decomposition.PCA` to project into a
# 2-dimensional space and plot the data and the clusters in this new space.
import matplotlib.pyplot as plt

reduced_data = PCA(n_components=2).fit_transform(data)
kmeans = KMeans(init="k-means++", n_clusters=n_digits, n_init=4)
kmeans.fit(reduced_data)

# Step size of the mesh. Decrease to increase the quality of the VQ.
h = 0.02  # point in the mesh [x_min, x_max]x[y_min, y_max].

# Plot the decision boundary. For that, we will assign a color to each
x_min, x_max = reduced_data[:, 0].min() - 1, reduced_data[:, 0].max() + 1
y_min, y_max = reduced_data[:, 1].min() - 1, reduced_data[:, 1].max() + 1
xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))

# Obtain labels for each point in mesh. Use last trained model.
Z = kmeans.predict(np.c_[xx.ravel(), yy.ravel()])

# Put the result into a color plot
Z = Z.reshape(xx.shape)
plt.figure(1)
plt.clf()
plt.imshow(
    Z,
    interpolation="nearest",
    extent=(xx.min(), xx.max(), yy.min(), yy.max()),
    cmap=plt.cm.Paired,
    aspect="auto",
    origin="lower",
)

plt.plot(reduced_data[:, 0], reduced_data[:, 1], "k.", markersize=2)
# Plot the centroids as a white X
centroids = kmeans.cluster_centers_
plt.scatter(
    centroids[:, 0],
    centroids[:, 1],
    marker="x",
    s=169,
    linewidths=3,
    color="w",
    zorder=10,
)
plt.title(
    "K-means clustering on the digits dataset (PCA-reduced data)\n"
    "Centroids are marked with white cross"
)
plt.xlim(x_min, x_max)
plt.ylim(y_min, y_max)
plt.xticks(())
plt.yticks(())
plt.show()
image-20230204175455318
image-20230204175505770

总结

kmeans算法具有速度快，易理解的特点，但是同时也有初始点选点的问题。