Python最小二乘法拟合与作图

在函数拟合中，如果用p表示函数中需要确定的参数，那么目标就是找到一组p，使得下面函数S的值最小：

为维西等地区用户提供了全套网页设计制作服务，及维西网站建设行业解决方案。主营业务为网站设计制作、网站建设、维西网站设计，以传统方式定制建设网站，并提供域名空间备案等一条龙服务，秉承以专业、用心的态度为用户提供真诚的服务。我们深信只要达到每一位用户的要求，就会得到认可，从而选择与我们长期合作。这样，我们也可以走得更远！

这种算法称为最小二乘法拟合。Python的Scipy数值计算库中的optimize模块提供了 leastsq() 函数，可以对数据进行最小二乘拟合计算。

此处利用该函数对一段弧线使用圆方程进行了拟合，并通过Matplotlib模块进行了作图，程序内容如下：

Python的使用中需要导入相应的模块，此处首先用 import 语句

分别导入了numpy, leastsq与pylab模块，其中numpy模块常用用与数组类型的建立，读入等过程。leastsq则为最小二乘法拟合函数。pylab是绘图模块。

接下来我们需要读入需要进行拟合的数据，这里使用了 numpy.loadtxt() 函数：

其参数有：

进行拟合时，首先我们需要定义一个目标函数。对于圆的方程，我们需要圆心坐标(a,b)以及半径r三个参数，方便起见用p来存储：

紧接着就可以进行拟合了， leastsq() 函数需要至少提供拟合的函数名与参数的初始值：

返回的结果为一数组，分别为拟合得到的参数与其误差值等，这里只取拟合参数值。

leastsq() 的参数具体有：

输出选项有：

最后我们可以将原数据与拟合结果一同做成线状图，可采用 pylab.plot() 函数：

pylab.plot() 函数需提供两列数组作为输入，其他参数可调控线条颜色，形状，粗细以及对应名称等性质。视需求而定，此处不做详解。

pylab.legend() 函数可以调控图像标签的位置，有无边框等性质。

pylab.annotate() 函数设置注释，需至少提供注释内容与放置位置坐标的参数。

pylab.show() 函数用于显示图像。

最终结果如下图所示：

用Python作科学计算

numpy.loadtxt

scipy.optimize.leastsq

python 常用的系统函数有哪些

1.常用内置函数：(不用import就可以直接使用)

help(obj) 在线帮助, obj可是任何类型

callable(obj) 查看一个obj是不是可以像函数一样调用

repr(obj) 得到obj的表示字符串，可以利用这个字符串eval重建该对象的一个拷贝

eval_r(str) 表示合法的python表达式，返回这个表达式

dir(obj) 查看obj的name space中可见的name

hasattr(obj,name) 查看一个obj的name space中是否有name

getattr(obj,name) 得到一个obj的name space中的一个name

setattr(obj,name,value) 为一个obj的name space中的一个name指向vale这个object

delattr(obj,name) 从obj的name space中删除一个name

vars(obj) 返回一个object的name space。用dictionary表示

locals() 返回一个局部name space,用dictionary表示

globals() 返回一个全局name space,用dictionary表示

type(obj) 查看一个obj的类型

isinstance(obj,cls) 查看obj是不是cls的instance

issubclass(subcls,supcls) 查看subcls是不是supcls的子类

类型转换函数

chr(i) 把一个ASCII数值,变成字符

ord(i) 把一个字符或者unicode字符,变成ASCII数值

oct(x) 把整数x变成八进制表示的字符串

hex(x) 把整数x变成十六进制表示的字符串

str(obj) 得到obj的字符串描述

list(seq) 把一个sequence转换成一个list

tuple(seq) 把一个sequence转换成一个tuple

dict(),dict(list) 转换成一个dictionary

int(x) 转换成一个integer

long(x) 转换成一个long interger

float(x) 转换成一个浮点数

complex(x) 转换成复数

max(...) 求最大值

min(...) 求最小值

用于执行程序的内置函数

complie 如果一段代码经常要使用,那么先编译,再运行会更快。

2.和操作系统相关的调用

系统相关的信息模块 import sys

sys.argv是一个list,包含所有的命令行参数.

sys.stdout sys.stdin sys.stderr 分别表示标准输入输出,错误输出的文件对象.

sys.stdin.readline() 从标准输入读一行 sys.stdout.write("a") 屏幕输出a

sys.exit(exit_code) 退出程序

sys.modules 是一个dictionary，表示系统中所有可用的module

sys.platform 得到运行的操作系统环境

sys.path 是一个list,指明所有查找module，package的路径.

操作系统相关的调用和操作 import os

os.environ 一个dictionary 包含环境变量的映射关系 os.environ["HOME"] 可以得到环境变量HOME的值

os.chdir(dir) 改变当前目录 os.chdir('d:\\outlook') 注意windows下用到转义

os.getcwd() 得到当前目录

os.getegid() 得到有效组id os.getgid() 得到组id

os.getuid() 得到用户id os.geteuid() 得到有效用户id

os.setegid os.setegid() os.seteuid() os.setuid()

os.getgruops() 得到用户组名称列表

os.getlogin() 得到用户登录名称

os.getenv 得到环境变量

os.putenv 设置环境变量

os.umask 设置umask

os.system(cmd) 利用系统调用，运行cmd命令

操作举例：

os.mkdir('/tmp/xx') os.system("echo 'hello' /tmp/xx/a.txt") os.listdir('/tmp/xx')

os.rename('/tmp/xx/a.txt','/tmp/xx/b.txt') os.remove('/tmp/xx/b.txt') os.rmdir('/tmp/xx')

用python编写一个简单的shell

#!/usr/bin/python

import os, sys

cmd = sys.stdin.readline()

while cmd:

os.system(cmd)

cmd = sys.stdin.readline()

用os.path编写平台无关的程序

os.path.abspath("1.txt") == os.path.join(os.getcwd(), "1.txt")

os.path.split(os.getcwd()) 用于分开一个目录名称中的目录部分和文件名称部分。

os.path.join(os.getcwd(), os.pardir, 'a', 'a.doc') 全成路径名称.

os.pardir 表示当前平台下上一级目录的字符 ..

os.path.getctime("/root/1.txt") 返回1.txt的ctime(创建时间)时间戳

os.path.exists(os.getcwd()) 判断文件是否存在

os.path.expanduser('~/dir') 把~扩展成用户根目录

os.path.expandvars('$PATH') 扩展环境变量PATH

os.path.isfile(os.getcwd()) 判断是否是文件名，1是0否

os.path.isdir('c:\Python26\temp') 判断是否是目录,1是0否

os.path.islink('/home/huaying/111.sql') 是否是符号连接 windows下不可用

os.path.ismout(os.getcwd()) 是否是文件系统安装点 windows下不可用

os.path.samefile(os.getcwd(), '/home/huaying') 看看两个文件名是不是指的是同一个文件

os.path.walk('/home/huaying', test_fun, "a.c")

遍历/home/huaying下所有子目录包括本目录,对于每个目录都会调用函数test_fun.

例：在某个目录中，和他所有的子目录中查找名称是a.c的文件或目录。

def test_fun(filename, dirname, names): //filename即是walk中的a.c dirname是访问的目录名称

if filename in names: //names是一个list,包含dirname目录下的所有内容

print os.path.join(dirname, filename)

os.path.walk('/home/huaying', test_fun, "a.c")

文件操作

打开文件

f = open("filename", "r") r只读 w写 rw读写 rb读二进制 wb写二进制 w+写追加

读写文件

f.write("a") f.write(str) 写一字符串 f.writeline() f.readlines() 与下read类同

f.read() 全读出来 f.read(size) 表示从文件中读取size个字符

f.readline() 读一行,到文件结尾,返回空串. f.readlines() 读取全部，返回一个list. list每个元素表示一行，包含"\n"\

f.tell() 返回当前文件读取位置

f.seek(off, where) 定位文件读写位置. off表示偏移量，正数向文件尾移动，负数表示向开头移动。

where为0表示从开始算起,1表示从当前位置算,2表示从结尾算.

f.flush() 刷新缓存

关闭文件

f.close()

regular expression 正则表达式 import re

简单的regexp

p = re.compile("abc") if p.match("abc") : print "match"

上例中首先生成一个pattern(模式),如果和某个字符串匹配，就返回一个match object

除某些特殊字符metacharacter元字符，大多数字符都和自身匹配。

这些特殊字符是 。^ $ * + ? { [ ] \ | ( )

字符集合(用[]表示)

列出字符,如[abc]表示匹配a或b或c,大多数metacharacter在[]中只表示和本身匹配。例：

a = ".^$*+?{\\|()" 大多数metachar在[]中都和本身匹配，但"^[]\"不同

p = re.compile("["+a+"]")

for i in a:

if p.match(i):

print "[%s] is match" %i

else:

print "[%s] is not match" %i

在[]中包含[]本身，表示"["或者"]"匹配.用

和

表示.

^出现在[]的开头,表示取反.[^abc]表示除了a,b,c之外的所有字符。^没有出现在开头，即于身身匹配。

-可表示范围.[a-zA-Z]匹配任何一个英文字母。[0-9]匹配任何数字。

\在[]中的妙用。

\d [0-9]

\D [^0-9]

\s [ \t\n\r\f\v]

\S [^ \t\n\r\f\v]

\w [a-zA-Z0-9_]

\W [^a-zA-Z0-9_]

\t 表示和tab匹配, 其他的都和字符串的表示法一致

\x20 表示和十六进制ascii 0x20匹配

有了\，可以在[]中表示任何字符。注：单独的一个"."如果没有出现[]中，表示出了换行\n以外的匹配任何字符,类似[^\n].

regexp的重复

{m,n}表示出现m个以上(含m个),n个以下(含n个). 如ab{1,3}c和abc,abbc,abbbc匹配，不会与ac,abbbc匹配。

m是下界，n是上界。m省略表下界是0,n省略，表上界无限大。

*表示{,} +表示{1,} ?表示{0,1}

最大匹配和最小匹配 python都是最大匹配，如果要最小匹配，在*,+,?,{m,n}后面加一个?.

match object的end可以得到匹配的最后一个字符的位置。

re.compile("a*").match('aaaa').end() 4 最大匹配

re.compile("a*?").match('aaaa').end() 0 最小匹配

使用原始字符串

字符串表示方法中用\\表示字符\.大量使用影响可读性。

解决方法：在字符串前面加一个r表示raw格式。

a = r"\a" print a 结果是\a

a = r"\"a" print a 结果是\"a

使用re模块

先用re.compile得到一个RegexObject 表示一个regexp

后用pattern的match,search的方法,得到MatchObject

再用match object得到匹配的位置,匹配的字符串等信息

RegxObject常用函数:

re.compile("a").match("abab") 如果abab的开头和re.compile("a")匹配，得到MatchObject

_sre.SRE_Match object at 0x81d43c8

print re.compile("a").match("bbab")

None 注：从str的开头开始匹配

re.compile("a").search("abab") 在abab中搜索第一个和re_obj匹配的部分

_sre.SRE_Match object at 0x81d43c8

print re.compile("a").search("bbab")

_sre.SRE_Match object at 0x8184e18 和match()不同,不必从开头匹配

re_obj.findall(str) 返回str中搜索所有和re_obj匹配的部分.

返回一个tuple,其中元素是匹配的字符串.

MatchObject的常用函数

m.start() 返回起始位置,m.end()返回结束位置(不包含该位置的字符).

m.span() 返回一个tuple表示(m.start(), m.end())

m.pos(), m.endpos(), m.re(), m.string()

m.re().search(m.string(), m.pos(), m.endpos()) 会得到m本身

m.finditer()可以返回一个iterator,用来遍历所有找到的MatchObject.

for m in re.compile("[ab]").finditer("tatbxaxb"):

print m.span()

高级regexp

| 表示联合多个regexp. A B两个regexp，A|B表示和A匹配或者跟B匹配.

^ 表示只匹配一行的开始行首,^只有在开头才有此特殊意义。

$ 表示只匹配一行的结尾

\A 表示只匹配第一行字符串的开头 ^匹配每一行的行首

\Z 表示只匹配行一行字符串的结尾 $匹配第一行的行尾

\b 只匹配词的边界 例：\binfo\b 只会匹配"info" 不会匹配information

\B 表示匹配非单词边界

示例如下：

print re.compile(r"\binfo\b").match("info ") #使用raw格式 \b表示单词边界

_sre.SRE_Match object at 0x817aa98

print re.compile("\binfo\b").match("info ") #没有使用raw \b表示退格符号

None

print re.compile("\binfo\b").match("\binfo\b ")

_sre.SRE_Match object at 0x8174948

分组(Group) 示例：re.compile("(a(b)c)d").match("abcd").groups() ('abc', 'b')

#!/usr/local/bin/python

import re

x = """

name: Charles

Address: BUPT

name: Ann

Address: BUPT

"""

#p = re.compile(r"^name:(.*)\n^Address:(.*)\n", re.M)

p = re.compile(r"^name:(?P.*)\n^Address:(?P.*)\n", re.M)

for m in p.finditer(x):

print m.span()

print "here is your friends list"

print "%s, %s"%m.groups()

Compile Flag

用re.compile得到RegxObject时，可以有一些flag用来调整RegxObject的详细特征.

DOTALL, S 让.匹配任意字符,包括换行符\n

IGNORECASE, I 忽略大小写

LOCALES, L 让\w \W \b \B和当前的locale一致

MULTILINE, M 多行模式，只影响^和$(参见上例)

VERBOSE, X verbose模式

Python 函数 - 返回生成器

如果函数要返回一系列结果，我们常见的方法就是将结果放到一份列表中，然后返回给调用者。比如下面的函数，返回字符串中每个单词的首字母在真个字符串中的索引：

运行结果：

上述的结果完全符合我们的预期，但 get_word_index 函数不够简洁。下面我们尝试使用生成器来实现：

运行结果：

改写之后，不仅运行结果符合要求，由于不需要和 result 列表交互，函数也变得非常简洁。下面我们就来详细学习下生成器吧~

生成器是指使用 yield 表达式的函数，调用生成器函数时，它并不会真的运行，而是会返回迭代器。每次在这个迭代器上面调用内置的 next 函数时，迭代器就会把生成器推进到下一个 yield 表达式那里。生成器传给 yield 的值均会由迭代器返回给调用者。

此外，如果输入量非常大，使用列表作为返回值，那么程序就有可能耗尽内存并崩溃。相反，使用生成器之后，则可以应对任意长度的输入数据。

例如，下面这个生成器函数可以获取文件中单词的索引，而不管文件内容多大，该函数执行时消耗的内存，只由单行的文本长度决定：

其中 test_generator.txt 中的内容如下：

运行结果：

下面这句话特别重要： 生成器函数返回的迭代器，是由状态的，及调用者不应该反复使用它 。我们那 word_index_iter 来说明：

如果想重复调用，请将其封装成容器：

运行结果：

关于上述自定义容器的实现原理，我的另外一篇文章做了详细介绍，链接奉上：

Python高级数据结构——堆

在一个 最小堆 (min heap) 中，如果 P 是 C 的一个父级节点，那么 P 的 key（或 value) 应小于或等于 C 的对应值。 正因为此，堆顶元素一定是最小的，我们会利用这个特点求最小值或者第 k 小的值。

在一个 最大堆 (max heap) 中，P 的 key（或 value) 大于或等于 C 的对应值。

以python为例，说明堆的几个常见操作，这里需要用到一个内置的包：heapq

python中使用堆是通过传入一个数组，然后调用一个函数，在原地让传入的数据具备堆的特性

需要注意的是，heapify默认构造的是小顶堆（min heap），如果要构造大顶堆，思路是把所有的数值倒转，既* -1，例如：

使用heapq提供的函数： heappop 来实现

具体使用方式参考 初始化Heapify

使用heapq提供的函数： heappush 来实现

同时heapq还提供另外一个函数： heappushpop ，能够在一个函数实现pushpop两个操作；顺序是：先push再pop

根据官方文档的描述，这个函数会比先在外围先调用heappush，再调用heappop，效率更高

先pop数据再push数据，和heappushpop的顺序是反着的； 同样的，这样调用的性能也会比先调用heappop再调用heappush更好

如果pop的时候队列是空的，会抛出一个异常

可以通过 heapq.merge 将多个 已排序 的输入合并为一个已排序的输出，这个本质上不是堆；其实就是用两个指针迭代

对于这个问题，有一个算法题可以实现相同的功能

从 iterable 所定义的数据集中返回前 n 个最大/小元素组成的列表。

函数为： heapq.nlargest() | heapq.nsmallest()

heapq - Heap queue algorithm - Python 3.10.4 documentation

Python气象数据处理与绘图(2)：常用数据计算方法

对于气象绘图来讲，第一步是对数据的处理，通过各类公式，或者统计方法将原始数据处理为目标数据。

按照气象统计课程的内容，我给出了一些常用到的统计方法的对应函数：

在计算气候态，区域平均时均要使用到求均值函数，对应NCL中的dim_average函数，在python中通常使用np.mean()函数

numpy.mean(a, axis, dtype)

假设a为[time,lat,lon]的数据，那么

需要特别注意的是，气象数据中常有缺测，在NCL中，使用求均值函数会自动略过，而在python中，当任意一数与缺测(np.nan)计算的结果均为np.nan，比如求[1,2,3,4，np.nan]的平均值，结果为np.nan

因此，当数据存在缺测数据时，通常使用np.nanmean()函数，用法同上，此时[1,2,3,4，np.nan]的平均值为(1+2+3+4)/4 = 2.5

同样的，求某数组最大最小值时也有np.nanmax(), np.nanmin()函数来补充np.max(), np.min()的不足。

其他很多np的计算函数也可以通过在前边加‘nan’来使用。

另外，

也可以直接将a中缺失值全部填充为0。

np.std(a, axis, dtype)

用法同np.mean()

在NCL中有直接求数据标准化的函数dim_standardize()

其实也就是一行的事，根据需要指定维度即可。

皮尔逊相关系数：

相关可以说是气象科研中最常用的方法之一了，numpy函数中的np.corrcoef(x, y)就可以实现相关计算。但是在这里我推荐scipy.stats中的函数来计算相关系数：

这个函数缺点和有点都很明显，优点是可以直接返回相关系数R及其P值，这避免了我们进一步计算置信度。而缺点则是该函数只支持两个一维数组的计算，也就是说当我们需要计算一个场和一个序列的相关时，我们需要循环来实现。

其中a[time,lat,lon]，b[time]

(NCL中为regcoef()函数)

同样推荐Scipy库中的stats.linregress(x,y)函数：

slop: 回归斜率

intercept：回归截距

r_value： 相关系数

p_value： P值

std_err： 估计标准误差

直接可以输出P值，同样省去了做置信度检验的过程，遗憾的是仍需同相关系数一样循环计算。

分享标题：python的p32函数 python2与3
新闻来源：http://hbruida.cn/article/hehheg.html