julia 入门学习

1.打印 println

2.字符串 $ 借鉴perl

3.特点

极其卓越的性能。Julia在很多数据分析任务以及其他编程实践中都表现出了令人难以置信的性能。它的表现可以和C语言媲美，C语言经常被用来作为衡量运算速度的标准。

强大的基础库。Julia有一个强大的基础库，它不需要其他平台，就可以进行所有的线性代数运算，这些运算是数据分析模块的必备组件。

支持多分派。Julia实现了多分派机制，这使它可以使用同一种函数实现不同的过程，使函数更容易扩展，并可以对不同类型的输入重复使用。

容易上手。特别是对于那些从Python、R、Matlab或Octave迁移过来的使用者，学习Julia特别容易。

用户友好的界面。不论是在本地还是云上，Julia的用户界面都非常友好，在所有的流程中，用户与Julia的交流都非常顺畅。Julia还对所有的功能和数据类型提供了方便易用的帮助文件。

与其他语言无缝对接。这些语言包括（但不限于）R、Python和C。这使你不需要进行完整的迁移，就可以使用现有的代码库。

开源。Julia以及它的所有文档与教程都是开源的，非常易于获取，详尽而又全面。

开发者承诺。Julia的开发者承诺会一直加强这门语言的性能，并对使用者提供尽可能的帮助。他们提供了大量的讨论，组织年度会议，并提供咨询服务。

自定义函数。Julia的自定义函数可以和内置在基础代码中的函数一样快速而简洁。

并行能力。Julia具有强大的并行能力，这使得在多核计算机和集群上的部署非常容易。

极大的灵活性。Julia在开发新程序方面极其灵活，不论是编程新手，还是专家级用户，Julia适合各种编程水平的使用者，这个特性在其他语言中是很难得的。

• 变量

变量名必须以字母（a-z 或 A-Z），下划线，或一个 Unicode 编码指针中指向比 00A0 更大的指针子集开始；特别是 Unicode 字符 Lu/Ll/Lt/Lm/Lo/Nl（字母），Sc/So （货币和其他符号），和其他一些可以看做字符的一些输入（例如 Sm 数学符号的子集）是允许的。首位之后的字符也包括 ！和数字（0-9 和其他字符 Nd/No ），以及其他 Unicode 编码指针：变音符号和其他修改标记（字母 Mn/Mc/Me/Sk），一些标点连接器（字母 PC），素数，和其他的一些字符。

运算符类似 + 也是有效的标识符，但需要特别解析。在某些情况下，运算符可以像变量一样使用；例如 (+) 是指增加功能，和 (+) = f 将重新定义这个运算。大多数的 Unicode 中缀操作符（在 Sm 中），如 ⊕ ，会被解析为中缀操作符，同时可以自定义方法（例如，你可以使用 ⊗ = kron 定义 ⊕ 成为一个中缀 Kronecker 积）。

尽管 Julia 对命名本身只有很少的限制, 但尽量遵循一定的命名规范吧：

变量名使用小写字母

单词间使用下划线 (‘_’) 分隔，但不鼓励

类型名首字母大写, 单词间使用驼峰式分隔.

函数名和宏名使用小写字母, 不使用下划线分隔单词.

修改参数的函数结尾使用 ! . 这样的函数被称为 mutating functions 或 in-place functions

• 字符串

关于 Julia 字符串，有一些值得注意的高级特性：

String 是个抽象类型，不是具体类型
Julia 的 Char 类型代表单字符，是由 32 位整数表示的 Unicode 码位
与 Java 中一样，字符串不可更改：String 对象的值不能改变。要得到不同的字符串，需要构造新的字符串
概念上，字符串是从索引值映射到字符的部分函数，对某些索引值，如果不是字符，会抛出异常
Julia 支持全部 Unicode 字符: 文本字符通常都是 ASCII 或 UTF-8 的，但也支持其它编码

Julia 函数的参数遵循 “pass-by-sharing” 的惯例，即不传递值，而是传递引用。函数参数本身，有点儿像新变量绑定（引用值的新位置），但它们引用的值与传递的值完全相同。对可变值（如数组）的修改，会影响其它函数。

Julia 提供一系列控制流：

复合表达式 ： begin 和 (;)
条件求值 ： if-elseif-else 和 ?: (ternary operator)
短路求值 ： &&, || 和 chained comparisons
重复求值: 循环 ： while 和 for
异常处理 ： try-catch ， error 和 throw
任务（也称为协程） ： yieldto

任务（也称为协程）

任务是一种允许计算灵活地挂起和恢复的控制流，有时也被称为对称协程、轻量级线程、协同多任务等。

如果一个计算（比如运行一个函数）被设计为 Task，有可能因为切换到其它 Task 而被中断。原先的 Task 在以后恢复时，会从原先中断的地方继续工作。切换任务不需要任何空间，同时可以有任意数量的任务切换，而不需要考虑堆栈问题。任务切换与函数调用不同，可以按照任何顺序来进行。

任务比较适合生产者-消费者模式，一个过程用来生产值，另一个用来消费值。消费者不能简单的调用生产者来得到值，因为两者的执行时间不一定协同。在任务中，两者则可以正常运行。

Julia 提供了 produce 和 consume 函数来解决这个问题。生产者调用 produce 函数来生产值：

整理部分代码

python常见代码的归纳总结

判断文件或者文件夹是否存在

if(os.path.exists(rootdir) == False)

创建文件夹

os.mkdir(rootdir)

调用系统命令

os.system(cmd)

字典循环

for key,value in dict.items()

打开文件并读取内容进行处理

fd = open('xxxx.txt', encoding='utf-8')
for line in fd:
    print line
fd.close()

创建文件并写入内容

fd = open('xxxx.txt', 'a+', encoding='utf-8')
fd.write('aaaaa' + '\n')
fd.close()

使用xlrd读取EXCEL

导入

import xlrd

打开excel

data = xlrd.open_workbook('demo.xls')     # 注意这里的workbook首字母是小写

查看文件中包含sheet的名称

data.sheet_names()

得到第一个工作表，或者通过索引顺序 或 工作表名称

table = data.sheets()[0]
table = data.sheet_by_index(0)
table = data.sheet_by_name(u'Sheet1')

获取行数和列数

nrows = table.nrows
ncols = table.ncols

获取整行和整列的值（数组）

table.row_values(i)
table.col_values(i)

循环行,得到索引的列表

for rownum in range(table.nrows):
    print table.row_values(rownum)

单元格

cell_A1 = table.cell(0,0).value
cell_C4 = table.cell(2,3).value

分别使用行列索引

cell_A1 = table.row(0)[0].value
cell_A2 = table.col(1)[0].value

简单的写入

row = 0
col = 0
ctype = 1 # 类型 0 empty,1 string, 2 number, 3 date, 4 boolean, 5 error
value = 'lixiaoluo'
xf = 0 # 扩展的格式化 (默认是0)
table.put_cell(row, col, ctype, value, xf)
table.cell(0,0) # 文本:u'lixiaoluo'
table.cell(0,0).value # 'lixiaoluo'

使用xlwt写入EXCEL

导入xlwt

import xlwt

新建一个excel文件

file = xlwt.Workbook() #注意这里的Workbook首字母是大写，无语吧

新建一个sheet

table = file.add_sheet('sheet name')

写入数据table.write(行,列,value)

table.write(0,0,'test')

如果对一个单元格重复操作，会引发

returns error:
# Exception: Attempt to overwrite cell:
# sheetname=u'sheet 1' rowx=0 colx=0

所以在打开时加cell_overwrite_ok=True解决

table = file.add_sheet('sheet name',cell_overwrite_ok=True)

保存文件

file.save('demo.xls')

另外，使用style

style = xlwt.XFStyle() #初始化样式

font = xlwt.Font() #为样式创建字体

font.name = 'Times New Roman'

font.bold = True

style.font = font #为样式设置字体

table.write(0, 0, 'some bold Times text', style) # 使用样式

命令行 getopt

try:
     options,args = getopt.getopt(sys.argv[1:],"hp:i:",["help","ip=","port="])
except getopt.GetoptError:
     sys.exit()
for name,value in options:
     if name in ("-h","--help"):
          usage()
     if name in ("-i","--ip"):
          print(value)
     if name in ("-p","--port"):
          print(value)

简单爬虫

import requests
AGENT = 'Mozilla/5.0 (Macintosh; Intel Mac OS X 10_10_3) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/52.0.2743.116 Safari/537.36'
HEADERS = {
        'User-Agent': AGENT,
        'Content-Type':'application/x-www-form-urlencoded; charset=UTF-8',
        'X-Requested-With':'XMLHttpRequest',
        'Accept':'*/*'
session = requests.session()

模拟登录

postdata = {
    'defaults':'xxx',
    'fromLogin':'xxx',
    'userName':'xxx',
    'password':'xxxx'
}
url = 'xxxxxxxx'
login_info = session.post(url, headers = HEADERS, data = postdata,verify = False)
if(login_info.status_code == requests.codes.ok):
    print('login success')
    return True
else:
    print('login err')
    return False
}

下载html页面

def downloadUrl(rootdir, url, orgid, page):
    html = session.get(url, headers=global_config.HEADERS, verify=False)
    if(html.text[1:7] == 'script'):
        print(html.text)
        return "err"
    if(len(html.text) < 60):
        return "err"
    sample = open(rootdir + "/" + str(orgid) + '_' + str(page) + ".html", "w", encoding='utf-8')
    sample.write(html.text)
    sample.close()
    return 'ok'

解析JOSN文件内容

def scrapy_by_file(json_file_name):
    #读取JSON文件的内容
    text = open(json_file_name, encoding='utf-8').read()
    #特殊处理,去除从WINDOWS系统带过来的BOM特殊字符
    if text.startswith(u'\ufeff'):
        text = text.encode('utf8')[3:].decode('utf8')
    #将文本内容的JSON数据转换成自定义的JSON对象
    try:
        json_data = json.loads(text)
    except:
        print(json_file_name)
        return
    for row in json_data['rows']:

def scrapy_by_row(row):
    try:
        orgid = row['organization']['id']
        familyid = row['censusRegisterFamily']['id']
    except:
        print('errrr')
        return
        scrapy_by_row(row)

遍历文件夹

遍历目录(rootdir) 遍历到的每个文件都执行dirFunc

def waklThroughDir(rootdir, dirFunc):
    for parent, dirnames, filenames in os.walk(rootdir):
        for filename in filenames:
            print(filename)
            #获取后缀为txt的文件
            if(filename.split('.')[-1] == 'html'):
                dirFunc(os.path.join(parent, filename))

##采集温州房产网基本信息

# -*- coding: utf-8 -*-
import re
import requests
import time

#-----------------------------用于解析的正则表达式常量------------------------------------------------------------------
#解析页数
PAGE_NUM = '共找到 (.*?) 符合条件的记录'
#解析小区名称
NAME = 'texttext_title"><ahref(.*?)</a></div><divclass="texttext_moreinfo">'
#解析小区价格
PRICE = 'class="hot_price">(.*?)</span>'
#解析小区地址
ADDRESS = 'text_moreinfo">(.*?)</div><divclass="texttext_moreinfo"><span>'
#文件生成路径
ROOTDIR = 'F:\\test\\'

#-----------------------------模拟请求的头部信息，否则将被识别出是程序抓包而被拦截--------------------------------------
HEADERS = {
    'Accept': 'text/html,application/xhtml+xml,application/xml;q=0.9,image/webp,*/*;q=0.8',
    'Accept-Encoding': 'gzip, deflate, sdch',
    'User-Agent': 'Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/51.0.2704.106 Safari/537.36',
    'Host': 'www.0577home.net',
    'Upgrade-Insecure-Requests': '1'
}

#-----------------------------抓取某一页的房产信息，pageNo为页号--------------------------------------------------------
def getHouseListByPageno(pageNo):
    #建立一个连接用于后续发起请求
    session = requests.session()
    url = 'http://www.0577home.net/xiaoqu/list_0_0_0_0_0_0_0_' + str(pageNo) + '.html'
    houseList = session.get(url, headers = HEADERS, verify = False)
    #以写入模式打开文件
    fh = open(ROOTDIR + "houseList_pageNo" + str(pageNo) + ".txt",  'w' ,encoding='utf-8')
    #将movieList写入文件
    fh.write(houseList.text)
    #关闭文件
    fh.close()

#-------------------------------获取需要抓取的页面总数------------------------------------------------------------------
def getPageNum():
    #打开已经下载好的第一页房产内容
    f = open(ROOTDIR + 'houseList_pageNo1.txt', encoding='utf-8')
    #获取文件内容
    rawContent = f.read()
    #用正则表达式解析页面内容
    pageNum = re.findall(PAGE_NUM, rawContent)
    #返回页面号
    return int(pageNum[0]) / 20 + 1

def parseHouseListToFile(srcFile, dstFile):
    #打开待解析的文件
    f = open(srcFile, encoding='utf-8')
    #读取文件内容以备解析
    rawContent = f.read()
    p = re.compile('\s+')
    content = re.sub(p, '', rawContent)
    dnames = re.findall(NAME, content)
    names = []
    for dname in dnames:
        idx = dname.rfind('>')
        names.append(dname[idx + 1:])
    prices = re.findall(PRICE, content)
    daddress = re.findall(ADDRESS, content)
    address = []
    for daddr in daddress:
        id = daddr.rfind('>')
        address.append(daddr[id + 1:])
    i = 0
    for x in names:
        #写入时用'$'做分割，结尾加上回车符
        dstFile.write(names[i] + '$' + prices[i] + '$' + address[i] + '\n')
        i = i + 1

# -------------------------------主函数，下载并解析房产信息--------------------------------------------------------------
if __name__ == '__main__':
    #---------------------抓取页面-----------------------------
    #抓取第一页房产信息
    getHouseListByPageno(1)
    #通过第一页房产信息获取总共要抓取的页面数量
    pageNum = getPageNum()
    #抓取剩余的页面
    for i in range(2, int(pageNum) + 1):
        getHouseListByPageno(str(i))
    #---------------------解析页面-----------------------------
    #获取当前年月日
    localtime = time.strftime('%Y%m%d', time.localtime(time.time()))
    #创建一个文件，文件名前面带上年月日
    f = open(ROOTDIR + localtime + '_houseList.txt', 'a+', encoding='utf-8')
    #解析所有的页面
    #for k in range(1, int(pageNum) + 1):
    for k in range(1, 115):
        parseHouseListToFile(ROOTDIR + "houseList_pageNo" + str(k) + ".txt", f)
    #关闭文件
    f.close()
    f = open(ROOTDIR + localtime + '_houseList.txt', encoding='utf-8')
    fd = open(ROOTDIR + localtime + '_houseInfo.txt', 'w', encoding='utf-8')
    k = 0
    for line in f:
        data = line.strip('\n')
        data = data.split('$')
        idx = data[3]
        getHouseInfoByPageno(idx, k)
        houseInfo = parseHouseInfo(ROOTDIR + "houseInfo_pageNo" + str(idx) + ".html")
        print(str(k) + "$".join(data) + '$' + "$".join(houseInfo))
        fd.write("$".join(data) + '$' + "$".join(houseInfo) + '\n')
        k += 1
    f.close()
    fd.close()

调用java

import sys
import jpype

name = sys.argv[1]
jarpath = '/home/dsadm/why/python'
jpype.startJVM(jpype.getDefaultJVMPath(), "-Djava.ext.dirs=%s" % jarpath)
DECRYPT = jpype.JClass('why.fmrt.decrypt.DECRYPT')
upperName =DECRYPT.decrypt(name)
print(upperName)
jpype.shutdownJVM()

构建 web 页面

import os

import tornado.httpserver
import tornado.ioloop
import tornado.options
import tornado.web

from view import *

from tornado.options import define, options
define("port", default=8000, help="run on the given port", type=int)

class Application(tornado.web.Application):
    def __init__(self):
        handlers = [
            (r"/", Indexhandler),
        ]
        settings = dict(
            template_path=os.path.join(os.path.dirname(__file__), 'templates'),
            autoescape=None,
            debug=False,
        )
        tornado.web.Application.__init__(self, handlers, **settings)

if __name__ == "__main__":
    tornado.options.parse_command_line()
    http_server = tornado.httpserver.HTTPServer(Application(), xheaders=True)
    http_server.listen(options.port)
    tornado.ioloop.IOLoop.instance().start()

多进程

import multiprocessing
for process_id in range(PROCESS_NUM):
    p = multiprocessing.Process(target=worker, args=(process_id,))
    jobs.append(p)
    p.start()

多线程

# -*- coding: utf-8 -*-

import threading

class Threadconfig():
    def __init__(self, thread_size):
        self.thread_size = thread_size

    def topen(self):
        self.thread_tasks = []

    def build(self, func, **kwargs):
        self.thread_task = threading.Thread(target=func, kwargs=(kwargs))
        self.thread_tasks.append(self.thread_task)

    def run(self):
        for thread_task in self.thread_tasks:
            thread_task.setDaemon(True)
            thread_task.start()
        while 1:
            alive = False
            for thread_num in range(0, self.thread_size):
                alive = alive or self.thread_tasks[thread_num].isAlive()
            if not alive:
                break

    def __del__(self):
        self.thread_tasks = []

移除中文分隔符号

  cmd = "sed ':a;N;$ s/\\r\\n//g;ba' " + oldfile + " > " + newfile
os.system(cmd)  

线性回归是所有模型的基石

手工玩偶店的小何童鞋，想分析下玩偶的数量和成本之间的关系，他得到了如下表格的数据。将这些数据放在了图像之上，似乎是一个线性的关系，但是感觉并不严格，像是在一条直线上下随机波动。其实数据是由“自然之力”按照下面的公式来产生的。 其中b是一个随机变量，服从期望为0，方差为1的正态分布。

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

# 10,7.7
# 10,9.87
# 11,11.18
# 12,10.43
# 13,12.36
# 14,14.15
# 15,15.73
# 16,16.4
# 17,18.86
# 18,16.13
# 19,18.21
# 20,18.37
# 21,22.61
# 22,19.83
# 23,22.67
# 24,22.7
# 25,25.16
# 26,25.55
# 27,28.21
# 28,28.12
# ...

$$y{i} = x{i} + b_{i}$$

机器学习的解法

步骤如下：
1.确定场景的类型
2.定义损失函数，使得模型预测的成本和实际的成本相近
3.提取特征（可能需要除去记错或者是特别异常的数据）
4.确定模型并估计参数（直接上线性模型）
5.评估模型效果（均方差要达到最小）

1

%matplotlib inline

1
2

import sys
sys.version

'3.6.4 (default, Jan 21 2018, 16:48:17) \n[GCC 4.2.1 Compatible Apple LLVM 9.0.0 (clang-900.0.39.2)]'

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

# _*_ coding=utf8 _*_

"""
机器学习方法实现
"""
import os
import sys

import matplotlib.pyplot as plt
import pandas as pd 
import numpy as np
from sklearn import linear_model


def data_load(path):

    # data = pd.read_csv(path)
    # data = pd.read_csv(path, delimiter='', header=0)
    # # 选择列
    # data['x']
    # data[['x', 'y']]
    # data.x
    #
    #
    # # 选择行
    # data[:10]
    # data[10:]
    # data[2:3]
    pass


def read_data(path):
#     data = pd.read_csv(path, delimiter=',', header=0)
    data = pd.read_csv(path)
    return data


def linearModel(data):
    features = ["x"]
    labels = ["y"]

    train_data = data[:15]
    test_data = data[15:]
    # 产生并且训练模型
    model = train_model(train_data, features, labels)
    # 评价模型效果
    error, score = evaluate_model(model, test_data, features, labels)
    # 图形化模型结果
    visual_model(model, data, features, labels, error, score)


def train_model(train_data, features, labels):
    model = linear_model.LinearRegression()
    model.fit(train_data[features], train_data[labels])
    return model


def evaluate_model(model, test_data, features, labels):
    error = np.mean(
        (model.predict(test_data[features]) - test_data[labels])**2)
    score = model.score(test_data[features], test_data[labels])
    return error, score


def visual_model(model, data, features, labels, error, score):
    """
    模型可视化
    """
    # 为在Matplotlib中显示中文，设置特殊字体
    plt.rcParams['font.family'] = ['Heiti']
    # 创建一个图形框
    fig = plt.figure(figsize=(6, 6), dpi=80)
    # 在图形框里只画一幅图
    ax = fig.add_subplot(111)
    # 在Matplotlib中显示中文，需要使用unicode
    # 在Python3中，str不需要decode
    if sys.version_info[0] == 3:
        ax.set_title(u'%s' % "线性回归示例")
    else:
        ax.set_title(u'%s' % "线性回归示例".decode("utf-8"))
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')
    # 画点图，用蓝色圆点表示原始数据
    # 在Python3中，str不需要decode
    if sys.version_info[0] == 3:
        ax.scatter(data[features], data[labels], color='b',
            label=u'%s: $y = x + \epsilon$' % "真实值")
    else:
        ax.scatter(data[features], data[labels], color='b',
            label=u'%s: $y = x + \epsilon$' % "真实值".decode("utf-8"))
    # 根据截距的正负，打印不同的标签
    if model.intercept_ > 0:
        # 画线图，用红色线条表示模型结果
        # 在Python3中，str不需要decode
        if sys.version_info[0] == 3:
            ax.plot(data[features], model.predict(data[features]), color='r',
                label=u'%s: $y = %.3fx$ + %.3f'\
                % ("预测值", model.coef_, model.intercept_))
        else:
            ax.plot(data[features], model.predict(data[features]), color='r',
                label=u'%s: $y = %.3fx$ + %.3f'\
                % ("预测值".decode("utf-8"), model.coef_, model.intercept_))
    else:
        # 在Python3中，str不需要decode
        if sys.version_info[0] == 3:
            ax.plot(data[features], model.predict(data[features]), color='r',
                label=u'%s: $y = %.3fx$ - %.3f'\
                % ("预测值", model.coef_, abs(model.intercept_)))
        else:
            ax.plot(data[features], model.predict(data[features]), color='r',
                label=u'%s: $y = %.3fx$ - %.3f'\
                % ("预测值".decode("utf-8"), model.coef_, abs(model.intercept_)))
    legend = plt.legend(shadow=True)
    legend.get_frame().set_facecolor('#6F93AE')
    # 显示均方差和决定系数
    # 在Python3中，str不需要decode
    if sys.version_info[0] == 3:
        ax.text(0.99, 0.01, 
            u'%s%.3f\n%s%.3f'\
            % ("均方差：", error, "决定系数：", score),
            style='italic', verticalalignment='bottom', horizontalalignment='right',
            transform=ax.transAxes, color='m', fontsize=13)
    else:
         ax.text(0.99, 0.01, 
            u'%s%.3f\n%s%.3f'\
            % ("均方差：".decode("utf-8"), error, "决定系数：".decode("utf-8"), score),
            style='italic', verticalalignment='bottom', horizontalalignment='right',
            transform=ax.transAxes, color='m', fontsize=13)
    # 展示上面所画的图片。图片将阻断程序的运行，直至所有的图片被关闭
    # 在Python shell里面，可以设置参数"block=False"，使阻断失效。
    plt.show()


# 为什么搭建模型 有两个目的 使用模型对未知数据做预测；或者利用模型分析数据，找出数据的内在规律，为决策提供支持。
# 在利用模型做预测时，我们希望模型的准确度越高越好，但是在利用模型做分析时， 我们则更关心模型是否是可靠的

1
2
3
4
5
6
7
8

if __name__ == '__main__':
#     home_path = os.path.dirname(os.path.abspath(__file__))
    if os.name == "posix":
        data_path = "./data/simple_example.csv"

    data = read_data(data_path)
    print(data.shape)
#     linearModel(data)

(20, 2)

1

linearModel(data)

/usr/local/Cellar/python3/3.6.4_2/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/scipy/linalg/basic.py:1226: RuntimeWarning: internal gelsd driver lwork query error, required iwork dimension not returned. This is likely the result of LAPACK bug 0038, fixed in LAPACK 3.2.2 (released July 21, 2010). Falling back to 'gelss' driver.
  warnings.warn(mesg, RuntimeWarning)

统计学

1.假设条件概率
2.估计参数
3.推导参数的分布
4.假设检验和置信区间

使用第三方库 statsmodels 来训练线性回归模型

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

import statsmodels.api as sm


def linear_model2(data):
    features = ["x"]
    labels = ["y"]
    Y = data[labels]
    X = sm.add_constant(data[features])
    result = train_model2(X, Y)
    model_summary(result)
    
    
def train_model2(X, Y):
    model = sm.OLS(Y, X)
    result = model.fit()
    return result


def model_summary(result):
    # 整体统计分析结果
    print(result.summary())
    print(result.f_test("x=0"))
    print(result.f_test("const=0"))
    print(result.f_test(["x=0", "const=0"]))

1
2
3

datapath = "./data/simple_example.csv"
data = read_data(datapath)
linear_model2(data)

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.962
Model:                            OLS   Adj. R-squared:                  0.960
Method:                 Least Squares   F-statistic:                     460.5
Date:                Tue, 07 Aug 2018   Prob (F-statistic):           2.85e-14
Time:                        11:40:29   Log-Likelihood:                -31.374
No. Observations:                  20   AIC:                             66.75
Df Residuals:                      18   BIC:                             68.74
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.9495      0.934     -1.017      0.323      -2.912       1.013
x              1.0330      0.048     21.458      0.000       0.932       1.134
==============================================================================
Omnibus:                        0.745   Durbin-Watson:                   2.345
Prob(Omnibus):                  0.689   Jarque-Bera (JB):                0.673
Skew:                           0.074   Prob(JB):                        0.714
Kurtosis:                       2.113   Cond. No.                         66.3
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
<F test: F=array([[ 460.4584822]]), p=2.848465414495684e-14, df_denom=18, df_num=1>
<F test: F=array([[ 1.03355794]]), p=0.32279564008314576, df_denom=18, df_num=1>
<F test: F=array([[ 2442.62159921]]), p=1.2108814742372977e-22, df_denom=18, df_num=2>

得到参数b的估计值为-0.9495, 但是这个值在b=0这个假设下的P-value高达32.3%，统计学上认为这种参数是不显著的，应该舍弃此参数。同理a的估计值是1.033，P-value小于0.01，因此a是显著的，应该被纳入模型。因此，需要调整模型。

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

import os
import sys

import numpy as np
import statsmodels.api as sm
from statsmodels.sandbox.regression.predstd import wls_prediction_std
import matplotlib.pyplot as plt
import pandas as pd


def modelSummary(re):
    """
    分析线性回归模型的统计性质
    """
    # 整体统计分析结果
    print(re.summary())
    # 在Windows下运行此脚本需确保Windows下的命令提示符(cmd)能显示中文
    # 用f test检测x对应的系数a是否显著
    print("检验假设x的系数等于0：")
    print(re.f_test("x=0"))
    # 用f test检测常量b是否显著
    print("检测假设const的系数等于0：")
    print(re.f_test("const=0"))
    # 用f test检测a=1, b=0同时成立的显著性
    print("检测假设x的系数等于1和const的系数等于0同时成立：")
    print(re.f_test(["x=1", "const=0"]))


def visualizeModel(re, data, features, labels):
    """
    模型可视化
    """
    # 计算预测结果的标准差，预测下界，预测上界
    prstd, preLow, preUp = wls_prediction_std(re, alpha=0.05)
    # 为在Matplotlib中显示中文，设置特殊字体
    plt.rcParams['font.family'] = ['Heiti']
    # 创建一个图形框
    fig = plt.figure(figsize=(6, 6), dpi=80)
    # 在图形框里只画一幅图
    ax = fig.add_subplot(111)
    # 在Matplotlib中显示中文，需要使用unicode
    # 在Python3中，str不需要decode
    if sys.version_info[0] == 3:
        ax.set_title(u'%s' % "线性回归统计分析示例")
    else:
        ax.set_title(u'%s' % "线性回归统计分析示例".decode("utf-8"))
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')
    # 画点图，用蓝色圆点表示原始数据
    # 在Python3中，str不需要decode
    if sys.version_info[0] == 3:
        ax.scatter(data[features], data[labels], color='b',
            label=u'%s: $y = x + \epsilon$' % "真实值")
    else:
        ax.scatter(data[features], data[labels], color='b',
            label=u'%s: $y = x + \epsilon$' % "真实值".decode("utf-8"))
    # 画线图，用红色虚线表示95%置信区间
    # 在Python3中，str不需要decode
    if sys.version_info[0] == 3:
        ax.plot(data[features], preUp, "r--", label=u'%s' % "95%置信区间")
        ax.plot(data[features], re.predict(data[features]), color='r',
            label=u'%s: $y = %.3fx$'\
            % ("预测值", re.params[features]))
    else:
        ax.plot(data[features], preUp, "r--", label=u'%s' % "95%置信区间".decode("utf-8"))
        ax.plot(data[features], re.predict(data[features]), color='r',
            label=u'%s: $y = %.3fx$'\
            % ("预测值".decode("utf-8"), re.params[features]))
    ax.plot(data[features], preLow, "r--")
    legend = plt.legend(shadow=True)
    legend.get_frame().set_facecolor('#6F93AE')
    plt.show()


def trainModel(X, Y):
    """
    训练模型
    """
    model = sm.OLS(Y, X)
    re = model.fit()
    return re


def linearModel2(data):
    """
    线性回归统计性质分析步骤展示

    参数
    ----
    data : DataFrame，建模数据
    """
    features = ["x"]
    labels = ["y"]
    Y = data[labels]
    # 加入常量变量
    X = sm.add_constant(data[features])
    # 构建模型
    re = trainModel(X, Y)
    # 分析模型效果
    modelSummary(re)
    # const并不显著，去掉这个常量变量
    resNew = trainModel(data[features], Y)
    # 输出新模型的分析结果
    print(resNew.summary())
    # 将模型结果可视化
    visualizeModel(resNew, data, features, labels)

1
2

linearModel2(data)

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.962
Model:                            OLS   Adj. R-squared:                  0.960
Method:                 Least Squares   F-statistic:                     460.5
Date:                Tue, 07 Aug 2018   Prob (F-statistic):           2.85e-14
Time:                        11:42:08   Log-Likelihood:                -31.374
No. Observations:                  20   AIC:                             66.75
Df Residuals:                      18   BIC:                             68.74
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.9495      0.934     -1.017      0.323      -2.912       1.013
x              1.0330      0.048     21.458      0.000       0.932       1.134
==============================================================================
Omnibus:                        0.745   Durbin-Watson:                   2.345
Prob(Omnibus):                  0.689   Jarque-Bera (JB):                0.673
Skew:                           0.074   Prob(JB):                        0.714
Kurtosis:                       2.113   Cond. No.                         66.3
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
检验假设x的系数等于0：
<F test: F=array([[ 460.4584822]]), p=2.848465414495684e-14, df_denom=18, df_num=1>
检测假设const的系数等于0：
<F test: F=array([[ 1.03355794]]), p=0.32279564008314576, df_denom=18, df_num=1>
检测假设x的系数等于1和const的系数等于0同时成立：
<F test: F=array([[ 0.99654631]]), p=0.3886267976063851, df_denom=18, df_num=2>
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.996
Model:                            OLS   Adj. R-squared:                  0.996
Method:                 Least Squares   F-statistic:                     4876.
Date:                Tue, 07 Aug 2018   Prob (F-statistic):           2.26e-24
Time:                        11:42:08   Log-Likelihood:                -31.933
No. Observations:                  20   AIC:                             65.87
Df Residuals:                      19   BIC:                             66.86
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
x              0.9862      0.014     69.825      0.000       0.957       1.016
==============================================================================
Omnibus:                        0.489   Durbin-Watson:                   2.218
Prob(Omnibus):                  0.783   Jarque-Bera (JB):                0.561
Skew:                           0.033   Prob(JB):                        0.755
Kurtosis:                       2.182   Cond. No.                         1.00
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

来对比下两种解法，数据科学家在搭建模型时，通常会定义一些技术指标来衡量模型预测的准确度。需要模型参数的估计值是可靠的，但是事实并非如此，机器学习的模型结果显示，玩偶生产的固定成本是负数，这违背了事实。模型没有抓住数据真正的内在关系。为了提高预测的准确度，常常提取更多的特征，并以此搭建复杂的模型。大家都热衷于复杂度更高的模型。一旦发生过拟合，模型越复杂，错得越多。

1.数据整合困难

半导体制造涉及的工序和设备非常多，各种英文名字的系统很多到现在都还没有完全搞清楚，数据也分散在各个系统中，有结构化的，也有非结构化的（各种wafer map），很多数据存在于设备的日志中，基本没有被使用过。即使要使用这些数据，也面临一些困难，比如有些设备数据我们没有权限访问，需要跟Vendor申请，有些数据Vendor做过加密，需要解密（如PDF的FDC raw trace数据，与厂商沟通很久，才能得到你想要的），有些设备数据服务器挂掉了，都没人知道，甚至找不到服务器维护人员（如Tel公司的WIS数据，应该很有价值，但是数据服务器直接挂掉好久了）。另外，这些数据一般缺少元数据，没有数据字典这种东西，你拿到这些数据，也要花费很长时间去理解数据，更不用说将所有这些系统的数据整合打通了。还好目前已经整清楚一部分数据了，包括传统的YMS数据（WAT,CP,Metrology,WIP）和tool数据（iEMS，FDC），后续有机会把半导体行业数据一点一点理清楚，设计一套公共数据模型（类似互联网数据中台），可以满足各种数据应用

2.组织架构不够完善

随着大数据兴起，公司高层还是非常重视大数据的，但在组织架构上并没有提升到很高的位置，毕竟对于半导体制造，大数据仍然属于边缘的Support部门。在IT部门下面有传统的EA（Engineering Analysis）部门，很多资源在做数据接口、数据解析、报表等工作，在YE，YAE等偏业务部门也有数据分析团队，他们更接近Fab的用户，有实际的use case，能积累一定的行业经验。现在公司又冒出了大数据团队，既有IT的EA下面一个小团队（主要是搞个大数据平台hadoop，在上面做一些数据分析和应用），也有YE下面的一个小团队（面向业务需求，提供大数据应用），这几波人经常没有形成合力，事情推动起来比较困难，IT里面搞大数据的，接触用户的机会比较少，经常会搞一些研究型的项目（如做各种correlation），并没有跟实际业务问题对接起来。业务部门的大数据团队面临数据整合的困难，面临巧妇难为无米之炊的困境，在数据没有整合好之前，很难做出大的成果，IT在数据方面是有优势的，它负责各种数据库系统，对数据更清楚一些。这几个团队如果能够整合成一个大数据部门，并在组织架构上上升一层，价值会体现的更加明显。另外，团队内部角色定位和分工也不是很清晰，按照EA1,EA2这种划分有意义吗？团队内部也没有形成很好的工作流程，各个小组工作互不相干，在一个部门内部甚至很少会有接触到。智能制造不仅体现在IT自动化和信息化水平上，未来更多要体现在数据应用上，要让数据提供智能决策和生产。

3、依赖外部大数据供应商

IT部门比较忙的事情，就是跟各种大数据分析厂商搞POC，采购各种分析工具，大数据分析技能没有提升多少，各种厂商的工具倒接触了不少（PDF的extensio，bistel eDatalyzer，nanometrics，spotfire），主要是部门没有资源做自主分析，很多大数据人才不会选择制造业，还有部门领导也不相信自己的团队能做出来，借着POC名义让大家学习人家是怎么做的。目前大数据分析厂商大多数是国外的，主要是国外半导体行业发展好，有大量的应用机会，这其实是国内大数据的一个机会，制造业的问题相对比较通用一些，是比较容易做出解决方案和软件工具的，随着国内半导体行业的发展，是有制造业大数据创业机会的。

4、传统方法使用大数据

半导体制造分析使用传统统计分析方法较多（DOE，ANOVA，boxplot，假设检验），这跟行业特点是有关系的，不像互联网更关注相关关系，制造业更关注因果关系，发现了问题，要找到问题的根本原因，才能解决问题，所以DOE的思想是比较好理解的。但是这些方法的缺陷也比较明显，半导体制程太复杂，影响因素太多，每个因素都做DOE，是非常低效的，这时候机器学习、深度学习的方法是可以大大提高解决问题效率的，用户对这块的接受度还是比较低的，他们对可解释性要求较高，其实机器学习有些基于规则的模型还是比较有价值的，要让用户接受这些将会是一个漫长的过程，要逐渐改变用户的观念。另外，在与用户沟通时，要主动引导用户，不要问用户需要什么（用户想要的就是各种correlation，boxplot，regression），要搞清楚用户面临的痛点，我们提供更好的解决方案帮他们解决就好了，而不是为用户提供他们所要的那些功能。

劝退文

先看一下作者的背景。作者叫 Alex Irpan，现为谷歌大脑机器人团队的软件工程师。他从伯克利拿到的计算机科学本科学位，本科的时候曾经在伯克利人工智能实验室（Berkeley AI Research (BAIR) Lab）进行本科科研，导师是 DRL 大牛 Pieter Abbeel，他还和 John Schulman 工作过。

这篇文章一上来就指出深度强化学习是个大坑。它的成功案例其实很少，但每个都太有名了，例如用 Deep Q Network（DQN）在 Atari games 上用原始像素图片作为状态达到甚至超越人类专家的表现、通过左右互搏（self-play）等方式在围棋上碾压人类、大大降低了谷歌能源中心的能耗等等。造成的结果就是没有从事过深度强化学习的研究人员对它产生了很大的错觉，高估了它的能力，低估了它的难度。

强化学习本身是一个非常通用的人工智能范式，在直觉上让人觉得非常适合用来模拟各种时序决策任务，如语音、文本类任务。当它和深度神经网络这种只要给我足够层和足够多的神经元，可以逼近任何函数的非线性函数近似模型结合在一起感觉要上天啊，无怪乎 DeepMind 经常号称人工智能=深度学习+强化学习。

然而 Alex 告诉我们别急，让我们先来审视一些问题：

1.它的样本利用率非常低。换言之为了让模型的表现达到一定高度需要极为大量的训练样本。

2.最终表现很多时候不够好。在很多任务上用非强化学习甚至非学习的其它方法，如基于模型的控制（model based control），线性二次型调节器（Linear Quadratic Regulator）等等可以获得好得多的表现。最气人的是这些模型很多时候样本利用率还高。当然这些模型有的时候会有一些假设比如有训练好的模型可以模仿，比如可以进行蒙特卡洛树搜索等等。

3.DRL 成功的关键离不开一个好的奖励函数（reward function），然而这种奖励函数往往很难设计。在 Deep Reinforcement Learning That Matters 作者提到有时候把奖励乘以一个常数模型表现就会有天和地的区别。但奖励函数的坑爹之处还不止如此。奖励函数的设计需要保证：

加入了合适的先验，良好的定义了问题和在一切可能状态下的对应动作。坑爹的是模型很多时候会找到作弊的手段。Alex 举的一个例子是有一个任务需要把红色的乐高积木放到蓝色的乐高积木上面，奖励函数的值基于红色乐高积木底部的高度而定。结果一个模型直接把红色乐高积木翻了一个底朝天。仔啊，你咋学坏了，阿爸对你很失望啊。

奖励函数的值太过稀疏。换言之大部分情况下奖励函数在一个状态返回的值都是 0。这就和我们人学习也需要鼓励，学太久都没什么回报就容易气馁。都说 21 世纪是生物的世纪，怎么我还没感觉到呢？21 世纪才刚开始呢。我等不到了啊啊啊啊啊。

有的时候在奖励函数上下太多功夫会引入新的偏见（bias）。

要找到一个大家都使用而又具有好的性质的奖励函数。这里Alex没很深入地讨论，但链接了一篇陶神（Terence Tao）的博客，大家有兴趣可以去看下。

4.局部最优/探索和剥削（exploration vs. exploitation）的不当应用。Alex举的一个例子是有一个连续控制的环境里，一个类似马的四足机器人在跑步，结果模型不小心多看到了马四脚朝天一顿乱踹后结果较好的情况，于是你只能看到四脚朝天的马了。

5.对环境的过拟合。DRL 少有在多个环境上玩得转的。你训练好的 DQN 在一个 Atari game上work 了，换一个可能就完全不 work。即便你想要做迁移学习，也没有任何保障你能成功。

6.不稳定性。

读 DRL 论文的时候会发现有时候作者们会给出一个模型表现随着尝试 random seed 数量下降的图，几乎所有图里模型表现最终都会降到 0。相比之下在监督学习里不同的超参数或多或少都会表现出训练带来的变化，而 DRL 里运气不好可能很长时间你模型表现的曲线都没有任何变化，因为完全不 work。

即便知道了超参数和随机种子，你的实现只要稍有差别，模型的表现就可以千差万别。这可能就是 Deep Reinforcement Learning That Matters 一文里 John Schulman 两篇不同文章里同一个算法在同一个任务上表现截然不同的原因。

即便一切都很顺利，从我个人的经验和之前同某 DRL 研究人员的交流来看只要时间一长你的模型表现就可能突然从很好变成完全不 work。原因我不是完全确定，可能和过拟合和 variance 过大有关。

特别是上述第六点，几乎是灾难性的。作者提到自己实习的时候一开始实现 Normalized Advantage Function (NAF)，为了找出 Theano 本身的 bugs 花了六周，这还是在 NAF 作者就在他旁边可以供他骚扰的情况下的结果。原因就是DRL的算法很多时候在没找好超参数的情况下就是不 work 的，所以你很难判断自己的代码到底有没有 bug 还是运气不好。

作者也回顾了 DRL 成功的案例，他认为 DRL 成功的案例其实非常少，大体包括：

各类游戏：Atari Games, Alpha Go/Alpha Zero/Dota2 1v1/超级马里奥/日本将棋，其实还应该有 DRL 最早的成功案例，93年的西洋双陆棋（backgammon）。

DeepMind 的跑酷机器人。

为 Google 的能源中心节能。

Google 的 AutoML。

作者认为从这些案例里获得的经验教训是 DRL 可能在有以下条件的情况下更可能有好的表现，条件越多越好：

• 数据获取非常容易，非常 cheap。

• 不要急着一上来就攻坚克难，可以从简化的问题入手。

• 可以进行左右互搏。

• 奖励函数容易定义。

• 奖励信号非常多，反馈及时。

他也指出了一些未来潜在的发展方向和可能性：

• 局部最优或许已经足够好。未来某些研究可能会指出我们不必过于担心大部分情况下的局部最优。因为他们比起全局最优并没有差很多。

• 硬件为王。在硬件足够强的情况下我们或许就不用那么在乎样本利用率了，凡事硬刚就可以有足够好的表现。各种遗传算法玩起来。

• 人为添加一些监督信号。在环境奖励出现频次太低的情况下可以引入自我激励（intrinsic reward）或者添加一些辅助任务，比如DeepMind就很喜欢这套，之前还写了一篇 Reinforcement Learning with Unsupervised Auxiliary Tasks（https://arxiv.org/abs/1611.05397） 。LeCun 不是嫌蛋糕上的樱桃太少吗，让我们多给他点樱桃吧！

• 更多融合基于模型的学习从而提高样本使用率。这方面的尝试其实已经有很多了，具体可以去看 Alex 提到的那些工作。但还远不够成熟。

• 仅仅把 DRL 用于 fine-tuning。比如最初 Alpha Go 就是以监督学习为主，以强化学习为辅。

• 自动学习奖励函数。这涉及到 inverse reinforcement learning 和 imitation learning。

• 迁移学习和强化学习的进一步结合。

• 好的先验。

• 有的时候复杂的任务反而更容易学习。Alex 提到的例子是 DeepMind 经常喜欢让模型学习很多同一环境的变种来减小对环境的过拟合。我觉得这也涉及 curriculum learning，即从简单的任务开始逐步加深难度。可以说是层层递进的迁移学习。另外一个可能的解释是很多时候人觉得困难的任务和机器觉得困难的任务是相反的。比如人觉得倒水很简单，你让机器人用学习的路子去学倒水就可以很难。但反过来人觉得下围棋很简单而机器学习模型却在下围棋上把人击败了。

最后 Alex 总体还是非常乐观的。他说尽管现在有很多困难，使得 DRL 或许还不是一个强壮（robust）到所有人都可以轻易加入的研究领域并且很多时候一些问题用DRL远没有监督学习简单和表现好，但或许过几年你再回来 DRL 就 work 了也未知啊。这还是很振奋人心的。田渊栋老师也表达过类似的想法，觉得正因为这个领域还不够成熟所以还有很多机会。他们都是了不起的研究人员。