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://vcdvsql.cn/article46/hehheg.html

成都網站建設公司_創新互聯，為您提供企業網站制作、營銷型網站建設、App開發、建站公司、定制網站、品牌網站建設

廣告

聲明：本網站發布的內容（圖片、視頻和文字）以用戶投稿、用戶轉載內容為主，如果涉及侵權請盡快告知，我們將會在第一時間刪除。文章觀點不代表本網站立場，如需處理請聯系客服。電話：028-86922220；郵箱：631063699@qq.com。內容未經允許不得轉載，或轉載時需注明來源： 創新互聯