看机器学习看到了回归函数，看了一半看不下去了，看到能用方差进行函数回归，又手痒痒了，自己推公式写代码验证：

常见的最小二乘法是一阶函数回归

回归方法就是寻找方差的最小值

y = kx + b

xi, yi

y-yi = kxi+b-yi

方差为

∑(kxi + b - yi )^2

f = k^2∑xi^2 + b^2 + ∑yi^2 +2kb∑xi - 2k∑xi*yi - 2yib

求极值需要对其求微分，因为是二元函数，因此使用全微分公式，其极值点应该在两个元的偏微分都为0处

δf/δk = 2k∑(xi^2)+2(∑xib-∑xi*yi) = 0

δf/δb = 2b +2(k∑xi -∑yi) = 0

b = ∑yi - k∑xi

2k∑(xi^2) + 2(∑xi∑yi - k(∑xi)^2-∑xi*yi)

k = (∑xi*yi - ∑xi∑yi)/(∑(xi^2)-(∑xi)^2)

以上就是最小二乘法的推导

那么扩展到二阶函数拟合回归

y = c0+c1*x+c2*x^2

xi, yi

y - yi = c0 + c1*xi + c2*xi^2 - yi

c0^2 + c1^2∑xi^2 + c2^2∑xi^4 +yi^2...

其三元偏微分表示如下

0 = c0 + c1∑xi + c2∑xi^2 - ∑yi

0 = c1∑xi^2 + c0∑xi + c2∑xi^3 - ∑xi*yi

0 = c2∑xi^4 + c0∑xi^2 + c1∑xi^3 - ∑xi^2 * yi

解方程可以求出三个参数

更高阶的回归可就不能手解方程了，需要用线性代数的知识

y = c0 + c1 * x + c2 * x^2 + c3 * x^3

0 = c0 + c1∑xi +c2∑xi^2 + c3∑xi^3 - ∑yi

0 = c0∑xi + c1∑xi^2 + c2∑xi^3 +c3∑xi^4 - ∑xi*yi

0 = c0∑xi^2 + c1∑xi^3 + c2∑xi^4 + c3∑xi^5 - ∑xi^2*yi

0 = c0∑xi^3 + c1∑xi^4 +c2∑xi^5 +c3∑xi^6 - ∑xi^3*yi

M * [c0, c1, c2, c3]T = [∑yi, ∑xi*yi, ∑xi^2*yi, ∑xi^3*yi]T

[c0, c1, c2, c3]T = M^-1 * [∑yi, ∑xi*yi, ∑xi^2*yi, ∑xi^3*yi]T

实际上这个最终转为了求矩阵逆的过程

因此多阶函数的拟合回归本质是矩阵求逆，这对于应用数学工具而言是比较容易的

更普遍的，对于y = ∑cj*gj(x)的形式，有

Yj = ∑(i)(gj(xi) * yi)

Mj,k = ∑(i)(gj(xi) * gk(xi))

M * C = Y; C = M^(-1) * Y (M 为矩阵，C和Y为向量)

求解Ck可以得到线性参数集

数学最优美的地方就是可以化繁为简，而同样优美的是代码，使用python写出以上过程的算法验证

以下程序基于numpy和matplotlib，其实简单应用用python确实足以替代matlab了

为了防止高次计算溢出，对x, y进行归一，之前没做这步，算到四次回归就溢出了

x' = x/xmax

y' = y/ymax

y = ∑(Cj * gj(x / xmax))*ymax

对：

x = 1, 2, 4, 6, 7, 9, 12, 15, 18, 20

y = 1.5, 5, 14, 22, 25, 29, 32, 37, 41, 45

进行3阶函数回归计算，完整代码如下，为了强制符合PEP8规范，英文注释写的很蹩脚...

终于上代码了：

"""

name: regression

author: robocky

create time: 2017-1-1

description: Use polynomial function to fitting data

"""

import numpy as np

from numpy.linalg import inv

import matplotlib.pyplot as plt

def polynomial(c: np.array, x: np.array, gs: list):

"""y = c0 * g0(x) + c1 * g1(x) + c2 * g2(x)...+ cn * gn(x)"""

return c.dot(np.array([g(x) for g in gs]))

def regression(x: np.array, y: np.array, gs: list):

"""Use inv of matrix to calculate regression

y = ∑ci * gi(x)

"""

y_res = np.array([sum(g(x) * y) for g in gs])

x_matrix = np.array([[sum(g(x) * h(x)) for g in gs] for h in gs])

return inv(x_matrix).dot(y_res)

def func_gen(n):

"""Generate power functions"""

return lambda x: x ** n

def func_gen_sin(n):

"""Generate sin functions"""

return (lambda x: x ** 0) if n == 0 else lambda x: np.sin(n * x)

if __name__ == '__main__':

# Test

x_list = np.array([1, 2, 4, 6, 7, 9, 12, 15, 18, 20])

y_list = np.array([1.5, 5, 14, 22, 25, 29, 32, 37, 41, 45])

g_list = [func_gen(i) for i in range(4)]

g_list2 = [func_gen_sin(i) for i in range(4)]

# Set numbers to unit, in order to avoid overflow

x_max, y_max = max(x_list), max(y_list)

c_list = regression(x_list / x_max, y_list / y_max, g_list)

c_list2 = regression(x_list / x_max, y_list / y_max, g_list2)

# plot

plt.figure()

plt.plot(x_list, y_list, 'bo')

x_line = np.arange(0, 1.02, 0.01)

plt.plot(x_line * x_max, polynomial(c_list, x_line, g_list) * y_max, color='red', label='x^n')

plt.plot(x_line * x_max, polynomial(c_list2, x_line, g_list2) * y_max, color='green', label='sin(x)')

plt.legend()

plt.show()

回归算法只有三行，当然计算过程都省略了，多项式也很简单

做了个函数生成器，理论上可以用任何函数进行回归运算，不过结果可能会发散，比如纯正弦函数，但加个常数项结果就收敛了

最后上个结果图

微博传播数量和传播深度的预测--基于pyspark和某个回归算法

8-28决定参加一下这个千万条的数据处理任务,因为场景和自己做过的一个回归分析预测差不多,第一天开始在小规模的数据上做准备工作. 第二次大修改版本 date 20160829 星期一¶ 原始数据处理, ...

从0开始学正则表达式－基于python

关于正则表达式,当我们了解它就不难,不了解就很难,其实任何事情都是这样,没有人一生下来就啥都会,说白了,每个人都是一个学习了解进步的过程.学习和掌握正则表达式可能并不是太简单,因为它确实是有点像“外星 ...

基于python的七种经典排序算法

参考书目: 一.排序的基本概念和分类 所谓排序,就是使一串记录,按照其中的某个或某些关键字的大小,递增或递减的排列起来的操作.排序算法,就是如何使得记录按照要求排列的方法. ...

Python‘ssl.match_hostname()’函数SSL证书验证安全绕过漏洞

漏洞名称: Python‘ssl.match_hostname()’函数SSL证书验证安全绕过漏洞 CNNVD编号: CNNVD-201312-033 发布时间: 2013-12-04 更新时间: 2 ...

基于python的七种经典排序算法(转)

一.排序的基本概念和分类 所谓排序,就是使一串记录,按照其中的某个或某些关键字的大小,递增或递减的排列起来的操作.排序算法,就是如何使得记录按照要求排列的方法. 排序的稳定性:经过某种排序后,如果两个 ...

基于Python协同过滤算法的认识

Contents    1. 协同过滤的简介    2. 协同过滤的核心    3. 协同过滤的实现    4. 协同过滤的应用 1. 协同过滤的简介 关于协同过滤的一个最经典的例子就是看电影,有时候 ...

基于面绘制的MC算法以及基于体绘制的 Ray-casting 实现Dicom图像的三维重建(python实现)

加入实验室后,经过张老师的介绍,有幸与某公司合共共同完成某个项目,在此项目中我主要负责的是三维 pdf 报告生成.Dicom图像上亮度.对比度调整以及 Dicom图像三维重建.今天主要介绍一下完成Di ...

机器学习_线性回归和逻辑回归_案例实战：Python实现逻辑回归与梯度下降策略_项目实战：使用逻辑回归判断信用卡欺诈检测

线性回归: 注:为偏置项,这一项的x的值假设为[1,1,1,1,1....] 注:为使似然函数越大,则需要最小二乘法函数越小越好 线性回归中为什么选用平方和作为误差函数?假设模型结果与测量值 误差满足 ...

从Theano到Lasagne：基于Python的深度学习的框架和库

从Theano到Lasagne:基于Python的深度学习的框架和库 摘要:最近,深度神经网络以“Deep Dreams”形式在网站中如雨后春笋般出现,或是像谷歌研究原创论文中描述的那样:Incept ...

随机推荐

KVM 介绍(3)：I/O 全虚拟化和准虚拟化 [KVM I/O QEMU Full-Virtualizaiton Para-virtualization]

学习 KVM 的系列文章: (1)介绍和安装 (2)CPU 和 内存虚拟化 (3)I/O QEMU 全虚拟化和准虚拟化(Para-virtulizaiton) (4)I/O PCI/PCIe设备直接分 ...

python chr() unichr() ord()

了解一下python chr(),unichr(),ord()函数的用法. 参考链接: http://crazier9527.iteye.com/blog/411001 chr() 输入参数(取值范围 ...

zw版_Halcon-delphi系列教程_卫星航拍精确打击目标自动识别

卫星航拍精确打击目标自动识别 这几天,俄罗斯空军在叙利亚大显神威,美军有意见了,说俄罗斯是狂轰滥炸,不是精确打击.      不过,战斗 ...

UVa 10161 Ant on a Chessboard

一道数学水题,找找规律. 首先要判断给的数在第几层,比如说在第n层.然后判断(n * n - n + 1)(其坐标也就是(n,n)) 之间的关系. 还要注意n的奇偶.  Problem A.Ant o ...

《Genesis-3D开源游戏引擎--横版格斗游戏制作教程02：关键帧动画导入与切割》

2. 关键帧动画导入与切割 动画的分割与导入概述: 在游戏当中,游戏角色在不同状态下会有不同的动作,这些动作在引擎里相当于一段段的动画片段.当导入模型资源的时候,连同模型动画都会一并导入到引擎中.开发 ...

div+css位置绝对定位和相对定位

position:absolute: 当div中被隔着些元素的话那么用此方法将把元素重叠在一起,所以元素可以不在容器中加也能重叠在一起

2017-2018-2 20155314《网络对抗技术》Exp6 信息搜集与漏洞扫描

2017-2018-2 20155314Exp6 信息搜集与漏洞扫描 目录 实验目标 实验内容 实验环境 基础问题回答 预备知识 实验步骤 1 信息搜集 1.1 外围信息搜 ...

Scala--包和引入

一.包 同一个包可以定义在多个文件中,你也可以在同一个文件中定义多个包 文件名为:pack.scala package com{ package scala{ package test{ class ...

Markdown 入门指南

导语: Markdown是一种轻量级的标记语言,语法简单,学习成本不算太高,但确实可以让你专注于文字,不用太分心与排版等等. Markdown 官方文档 这里可以看到官方的Markdown语法规则: ...

java web 验证码-字符变形(推荐)

该文章转载自:http://www.cnblogs.com/jianlun/articles/5553452.html 因为在我做的这个系统中发现验证码有点偏上,整体效果看起来不太好,就做了一些修改. ...