Python科学计算
统计-stats & 用Python做科学计算-基础篇
统计-stats
SciPy的stats模块包含了多种概率分布的随机变量，随机变量分为连续和离散两种。所有的连续随机变量都是rv_continuous的派生类的对象，而所有的离散随机变量都是rv_discrete的派生类的对象。
连续和离散概率分布
可以使用下面的语句获得stats模块中所有的连续随机变量：
&&& from scipy import stats
&&& [k for k,v in stats.__dict__.items() if isinstance(v, stats.rv_continuous)]
['genhalflogistic','triang','rayleigh','betaprime', ...]
连续随机变量对象都有如下方法：
rvs：对随机变量进行随机取值，可以通过size参数指定输出的数组的大小。
pdf：随机变量的概率密度函数。
cdf：随机变量的累积分布函数，它是概率密度函数的积分。
sf：随机变量的生存函数，它的值是1-cdf(t)。
ppf：累积分布函数的反函数。
stat：计算随机变量的期望值和方差。
fit：对一组随机取样进行拟合，找出最适合取样数据的概率密度函数的系数。
概率密度函数、直方图统计和累积分布函数
下面以正规分布为例，简单地介绍随机变量的用法。下面的语句获得缺省正规分布的随机变量的期望值和方差，我们看到缺省情况下它是一个均值为0，方差为1的随机变量：
&&& stats.norm.stats()
(array(0.0), array(1.0))
norm可以像函数一样调用，通过loc和scale参数可以指定随机变量的偏移和缩放参数。对于正态分布的随机变量来说，这两个参数相当于指定其期望值和标准差：
&&& X = stats.norm(loc=1.0, scale=2.0)
&&& x.stats()
(array(1.0), array(4.0))
下面调用随机变量X的rvs()方法，得到包含一万次随机取样值的数组x，然后调用NumPy的mean()和var()计算此数组的均值和方差，其结果符合随机变量X的特性：
&&& x = X.rvs(size=10000) # 对随机变量取10000个值
&&& np.mean(x) # 期望值
&&& np.var(x)
我们也可以使用fit()方法对随机取样序列x进行拟合，它返回的是与随机取样值最吻合的随机变量的参数：
&&& stats.norm.fit(x) # 得到随机序列期望值和标准差
array([ 1.,
接下来比较随机变量X的概率密度函数和对数组x进行直方图统计的结果：
&&& t = np.arange(-10, 10, 0.01)
&&& pl.plot(t, X.pdf(t)) # 绘制概率密度函数的理论值
&&& p, t2 = np.histogram(x, bins=100, normed=True)
&&& t2 = (t2[:-1] + t2[1:])/2
&&& pl.plot(t2, p) # 绘制统计所得到的概率密度
其中histogram()对数组x进行直方图统计，它将数组x的取值范围分为100个区间，并统计x中的每个值落入各个区间中的次数。histogram()返回两个数组p和t2，其中p表示各个区间的取样值出现的频数，由于normed参数为True，因此p的值是正规化之后的结果。t2表示区间，由于其中包括区间起点和终点，因此t2的长度为101。(左)显示了概率密度函数和直方图统计的结果，可以看出二者是一致的。
下面的程序绘制随机变量X的累积分布函数和数组p的累加结果，其结果如(右)所示。
&&& pl.plot(t, X.cdf(t))
&&& pl.plot(t2, np.add.accumulate(p)*(t2[1]-t2[0]))
正规分布的概率密度函数(左)和累积分布函数(右)
有些随机分布除了loc和scale参数之外，还需要额外的形状参数。例如伽玛分布可用于描述等待个独立的随机事件发生所需的时间，就是伽玛分布的形状参数。下面计算形状参数为1和2时的伽玛分布的期望值和方差：
&&& stats.gamma.stats(1.0)
(array(1.0), array(1.0))
&&& stats.gamma.stats(2.0)
(array(2.0), array(2.0))
伽玛分布的尺度参数和随机事件发生的频率相关，由scale参数指定：
&&& stats.gamma.stats(2.0, scale=2)
(array(4.0), array(8.0))
根据伽玛分布的数学定义可知其期望值为，方差为。上面的程序验证了这两个公式。
当随机分布有额外的形状参数时，它所对应的rvs()、pdf()等方法都会增加额外的参数接收形状参数。例如下面的程序调用rvs()对的伽玛分布取4个随机值：
&&& x = stats.gamma.rvs(2, scale=2, size=4)
array([ 2.,
接下来调用pdf()查看上面4个随机值所对应的概率密度：
&&& stats.gamma.pdf(x, 2, scale=2)
array([ 0.,
0.1498734 ,
也可以先创建将形状参数和尺度参数固定的随机变量，然后再调用其pdf()计算概率密度：
&&& X = stats.gamma(2, scale=2)
&&& X.pdf(x)
array([ 0.,
0.1498734 ,
当分布函数的值域为离散时我们称之为离散概率分布。例如投掷有六个面的骰子时，只能获得1到6的整数，因此所得到的概率分布为离散的。对于离散随机分布，通常使用概率质量函数(PMF)描述其分布情况。
在stats库中所有描述离散分布的随机变量都从rv_discrete类继承。也可以直接用rv_discrete类自定义离散概率分布。例如假设有一个不均匀的骰子，它的各点出现的概率不相等。我们可以用下面的数组x保存骰子的所有可能值，数组p保存每个值出现的概率：
&&& x = range(1,7)
&&& p = (0.4, 0.2, 0.1, 0.1, 0.1, 0.1)
于是可以用下面的语句定义表示这个特殊骰子的随机变量，并调用其rvs()方法投掷此骰子20次，获得符合概率p的随机数：
&&& dice = stats.rv_discrete(values=(x,p))
&&& dice.rvs(size=20)
array([2, 5, 1, 2, 1, 1, 2, 4, 1, 3, 1, 1, 4, 3, 1, 1, 1, 2, 6, 4])
二项、泊松、伽玛分布
本节用几个实例程序对概率论中的二项分布、泊松分布以及伽玛分布进行一些实验和讨论。
二项分布是最重要的离散概率分布之一。假设有一种只有两个结果的试验，其成功概率为，那么二项分布描述了进行次这样的独立试验，成功次的概率。二项分布的概率质量函数公式如下：
例如可以通过二项分布的概率质量公式计算投掷5次骰子出现3次6点的概率。投掷一次骰子，点数为6的概率(即试验成功的概率)为，试验次数为。使用二项分布的概率质量函数pmf()可以很容易计算出现次6点的概率。和概率密度函数pdf()类似，pmf()的第一个参数为随机变量的取值，后面的参数为描述随机分布所需的参数。对于二项分布来说，参数分别为和，而取值范围则为0到之间的整数。下面的程序计算为0到6所对应的概率：
&&& stats.binom.pmf(range(6), 5, 1/6.0)
array([0......000129])
由结果可知：出现0或1次6点的概率为40.2%，而出现3次6点的概率为3.215%。
在二项分布中，如果试验次数很大，而每次试验成功的概率很小，其乘积比较适中时，那么试验成功的次数的概率可以用泊松分布近似描述。
在泊松分布中使用描述单位时间(或单位面积)中随机事件的平均发生率。如果我们将二项分布中的试验次数看作单位时间中所做的试验次数，那么它和事件出现的概率的乘积就是事件的平均发生率，即。泊松分布的概率质量函数公式如下：
下面的程序分别计算二项分布和泊松分布的概率质量函数，其结果如所示，可以看出当足够大时，二者是十分接近的。
比较二项分布和泊松分布的概率质量函数
程序中事件平均发生率恒等于10。根据二项分布的试验次数，计算每次事件出现的概率。程序中运算部分大致如下：
&&& _lambda = 10.0
&&& k = np.arange(20)
&&& possion = stats.poisson.pmf(k, _lambda)
# 泊松分布
&&& binom100 = stats.binom.pmf(k, 100, _lambda/100)
# 二项分布 n=100
&&& binom1000 = stats.binom.pmf(k, 1000, _lambda/1000) # 二项分布 n=1000
&&& np.max(np.abs(binom100-possion)) # 计算最大误差
&&& np.max(np.abs(binom1000-possion)) # n为1000时，误差较小
当n足够大时二项分布和泊松分布近似相等
泊松分布适合描述单位时间内随机事件发生的次数的分布情况。例如某一设施在一定时间内的使用次数，机器出现故障的次数，自然灾害发生的次数等等。
为了加深读者对泊松分布概念的理解，下面我们使用随机数模拟泊松分布，并与其概率质量函数进行比较，其结果如所示。图中，每秒内事件的平均发生次数为10，即。其中左图的观察时间为1000秒，而右图的观察时间为50000秒。可以看出观察时间越长，每秒内事件发生的次数越符合泊松分布。
模拟泊松分布
模拟泊松分布
由于上面的程序中包含了许多绘图方面的代码，下面我们直接在IPython中介绍泊松分布的模拟过程。首先定义事件发生率和观察时间:
&&& _lambda = 10
&&& time = 10000
可以用NumPy的随机数生成函数rand()，产生平均分布于0到time之间的_lambda*time个事件所发生的时刻。由于rand()产生的是0到1之间的平均分布的随机数，因此需要对其结果扩大time倍：
&&& t = np.random.rand(_lambda*time)*time
用histogram()可以统计数组t中每秒之内的事件发生的次数count：
&&& count, time_edges = np.histogram(t, bins=time, range=(0,time))
array([10,
8, ..., 11, 10, 18])
根据泊松分布的定义，count数组中的数值的分布情况应该符合泊松分布。下面统计事件次数在0到20区间内的概率分布。当histogram()的normed参数为True并且每个统计区间的长度为1时，其结果和概率质量函数相等。
&&& dist, count_edges = np.histogram(count, bins=20, range=(0,20), normed=True)
&&& poisson = stats.poisson.pmf(x, _lambda)
&&& np.max(np.abs(dist-possion)) # 最大误差很小，符合泊松分布
还可以换一个角度看随机事件的分布问题。我们可以观察相邻两个事件之间的时间间隔的分布情况，或者隔k个事件的时间间隔的分布情况。根据概率论，事件之间的时间间隔应符合伽玛分布，由于时间间隔可以是任意数值，因此伽玛分布是一种连续概率分布。伽玛分布的概率密度函数公式如下，它描述第个事件发生所需的等待时间的概率分布。是伽玛函数，当为整数时，它的值和的阶乘相等。
下面的程序模拟了事件的时间间隔的伽玛分布，其结果如所示。图中的观察时间为1000秒，平均每秒产生10个事件。左图中“k=1”，它表示相邻两个事件间的间隔的分布，而“k=2”则表示相隔一个事件的两个事件间的间隔的分布，可以看出它们都符合伽玛分布。
模拟伽玛分布
模拟伽玛分布
下面我们直接在IPython中模拟伽玛分布。首先在10000秒之内产生100000个随机事件发生的时刻。因此事件的平均发生次数为每秒10次：
&&& _lambda = 10
&&& time = 10000
&&& t = np.random.rand(_lambda*time)*time
为了计算事件前后的时间间隔，需要先对随机时刻进行排序：
&&& t.sort()
然后分别计算“k=1”和“k=2”时的时间间隔：
&&& s1 = t[1:] - t[:-1]
#相邻两事件的时间间隔
&&& s2 = t[2:] - t[:-2]
#相隔一个事件的两个事件的时间间隔
对s1和s2分别调用histogram()进行概率统计，设置normed为True可以直接统计概率密度：
&&& dist1, x1 = np.histogram(s1, bins=100, normed=True)
&&& dist2, x2 = np.histogram(s2, bins=100, normed=True)
histogram()返回的第二个值为统计区间的边界，下面gamma.pdf()计算伽玛分布的概率密度时，使用各个区间的中值进行计算。pdf()的第二个参数为k值，scale参数为：
&&& gamma1 = stats.gamma.pdf((x1[:-1]+x1[1:])/2, 1, scale=1.0/_lambda)
&&& gamma2 = stats.gamma.pdf((x2[:-1]+x2[1:])/2, 2, scale=1.0/_lambda)
&&& np.max(np.abs(gamma1 - dist1))
&&& np.max(np.abs(gamma2 - dist2))
由于概率密度函数的值本身比较大，因此上面的误差已经很小了：
&&& np.max(gamma1), np.max(gamma2)
本站点由提供动力linux中安装python科学计算环境-numpy、scipy、matplotlib、OpenCV...
http://blog.csdn.net/pipisorry/article/details/
在Ubuntu中安装numpy、scipy、matplotlib、OpenCV等
和Python(x,y)不一样，在Ubuntu中需要手工安装科学计算的各个模块，
如何安装IPython, NumPy, SciPy, matplotlib, PyQt4, Spyder, Cython, SWIG, ETS, OpenCV：
在Ubuntu下安装Python模块通常可以使用apt-get和pip命令。apt-get命令是Ubuntu自带的包管理命令，而pip则是Python安装扩展模块的工具，通常pip会扩展模块的源代码并编译安装。
Ubuntu 12.04中缺省安装了Python2.7.3，首先通过下面的命令安装pip，pip是Python的一个安装和管理扩展库的工具。
sudo apt-get install python-pip
安装Python开发环境【亦可参见linux中搭建python开发环境安装pycharm IDE】，方便今后编译其他扩展库，占用空间92.8M：
sudo apt-get install python-dev
为了安装最新版的IPython 0.13beta，需要下载IPython源代码，并执行安装命令。在IPython 0.13beta中提供了改进版本的IPython notebook。下面的命令首先安装版本管理软件git，然后通过git命令从IPython的开发代码库中下载最新版本的IPython源代码，并执行安装命令：
sudo apt-get install git
git clone /ipython/ipython.git
cd ipython
sudo python setup.py install
如果安装目前的最新稳定版本，可以输入:
sudo apt-get install ipython
安装完毕之后，请输入ipython命令测试是否能正常启动。
为了让IPython notebook工作，还还需要安装tornado和pyzmq：
sudo pip install tornado
sudo apt-get install libzmq-dev
sudo pip install pyzmq
sudo pip install pygments
下面测试IPython：
mkdir notebook
cd notebook
ipython notebook
为了在IPython中离线使用LaTeX数学公式，需要安装mathjax，首先输入下面的命令启动ipython notebook：
sudo ipython notebook
在IPython notebook界面中输入：
from IPython.external.mathjax import install_mathjax
install_mathjax()
安装NumPy，SciPy和matplotlib
通过apt-get命令可以快速安装这三个库：
sudo apt-get install python-numpy
sudo apt-get install python-scipy
sudo apt-get install python-matplotlib
如果需要通过pip编译安装，可以先用apt-get命令安装所有编译所需的库：
sudo apt-get build-dep python-numpy
sudo apt-get build-dep python-scipy
然后通过pip命令安装：
sudo pip install numpy
sudo pip install scipy
通过build-dep会安装很多库，包括Python 3.2。
PyQt4和Spyder
下面的命令安装PyQt4，Qt界面设计器，PyQt4的开发工具以及文档：
sudo apt-get install python-qt4
sudo apt-get install qt4-designer
sudo apt-get install pyqt4-dev-tools
sudo apt-get install python-qt4-doc
安装完毕之后，文档位于：
/usr/share/doc/python-qt4-doc
安装好PyQt4之后通过下面的命令安装Spyder：
sudo apt-get install spyder
由于Spyder经常更新，通过下面的命令可以安装最新版：
sudo pip install spyder --upgrade
cython和SWIG
Cython和SWIG是编写Python扩展模块的工具：
sudo pip install cython
sudo apt-get install swig
输入 cython --version 和 swig -version 查看版本。
ETS是enthought公司开发的一套科学计算软件包，其中的Mayavi通过VTK实现数据的三维可视化。
首先通过下面的命令安装编译ETS所需的库：
sudo apt-get install python-dev libxtst-dev scons python-vtk pyqt4-dev-tools python2.7-wxgtk2.8 python-configobj
sudo apt-get install libgl1-mesa-dev libglu1-mesa-dev
创建ets目录，并在此目录下下载ets.py，运行ets.py可以复制最新版的ETS源程序，并安装：
wget /enthought/ets/raw/master/ets.py
python ets.py clone
sudo python ets.py develop
#sudo python ets.py install 或者运行install安装
如果一切正常，那么输入 mayavi2 命令则会启动mayavi。
为了编译OpenCV需要下载cmake编译工具，和一些依赖库：
sudo apt-get install build-essential
sudo apt-get install cmake
sudo apt-get install cmake-gui
sudo apt-get install libavcodec-dev libavformat-dev libswscale-dev
sudo apt-get install libjpeg-dev libpng-dev libtiff-dev libjer-dev
然后从 http://sourceforge.net/projects/opencvlibrary/ 下载最新版的OpenCV源代码，并解压。然后创建编译用的目录release，并启动cmake-gui:
mkdir release
在界面中选择OpenCV源代码的目录，和编译输出目录release，然后按Configure按钮，并根据需要设置各个编译选项，最后点Generate按钮，退出cmake-gui界面。进入编译路径，执行下面的命令：
cd release
sudo make install
安装完毕之后，启动I，并输入 import cv2 测试OpenCV是否能正常载入。
from:http://blog.csdn.net/pipisorry/article/details/
ref:Ubuntu中安装Python科学计算环境
Ubuntu-Python2.7安装 scipy,numpy,matplotlib
您对本文章有什么意见或着疑问吗？请到您的关注和建议是我们前行的参考和动力&&
您的浏览器不支持嵌入式框架，或者当前配置为不显示嵌入式框架。科学计算：最优化问题-1（Python）
scipy中的optimize子包中提供了常用的最优化算法函数实现。我们可以直接调用这些函数完成我们的优化问题。optimize中函数最典型的特点就是能够从函数名称上看出是使用了什么算法。下面optimize包中函数的概览：
1.非线性最优化
--& 简单Nelder-Mead算法
fmin_powell --& 改进型Powell法
fmin_bfgs& --&
拟Newton法
fmin_cg -- 非线性共轭梯度法
fmin_ncg& --&
线性搜索Newton共轭梯度法
leastsq& -- 最小二乘
2.有约束的多元函数问题
fmin_l_bfgs_b
---使用L-BFGS-B算法
---梯度信息
fmin_cobyla ---线性逼近
fmin_slsqp ---序列最小二乘法
nnls ---解|| Ax - b
||_2 for x&=0
3.全局优化
---模拟退火算法
4.标量函数
curve_fit--
使用非线性最小二乘法拟合
6.标量函数求根
brentq ---classic Brent
brenth ---A variation on the
classic Brent（1980）
---Ridder是提出这个算法的人名
bisect ---二分法
newton ---牛顿法
fixed_point
7.多维函数求根
fsolve ---通用
broyden1 ---Broyden’s first
Jacobian approximation.
broyden2 ---Broyden’s second
Jacobian approximation
newton_krylov ---Krylov
approximation for inverse Jacobian
anderson ---extended Anderson
excitingmixing ---tuned diagonal
Jacobian approximation
linearmixing ---scalar Jacobian
approximation
diagbroyden ---diagonal Broyden
Jacobian approximation
8.实用函数
line_search
---找到满足强Wolfe的alpha值
check_grad
---通过和前向有限差分逼近比较检查梯度函数的正确性
二、实战非线性最优化
fmin完整的调用形式是：
fmin(func, x0,
args=(), xtol=0.0001, ftol=0.0001,
maxiter=None, maxfun=None,
full_output=0, disp=1, retall=0,
callback=None)
不过我们最常使用的就是前两个参数。一个描述优化问题的函数以及初值。后面的那些参数我们也很容易理解。如果您能用到，请自己研究。下面研究一个最简单的问题，来感受这个函数的使用方法：f(x)=x**2-4*x+8，我们知道，这个函数的最小值是4，在x=2的时候取到。
from scipy.optimize import
#引入优化包
def myfunc(x):
x**2-4*x+8&&&
x0 = [1.3]&&&
#猜一个初值&
xopt = fmin(myfunc,
x0)&&& #求解
print xopt&&&
运行之后，给出的结果是：
Optimization terminated
successfully.
Current function value: 4.000000
Iterations: 16
Function evaluations: 32
程序准确的计算得出了最小值，不过最小值点并不是严格的2，这应该是由二进制机器编码误差造成的。
除了fmin_ncg必须提供梯度信息外，其他几个函数的调用大同小异，完全类似。我们不妨做一个对比：
from scipy.optimize import
fmin,fmin_powell,fmin_bfgs,fmin_cg
def myfunc(x):
return x**2-4*x+8
x0 = [1.3]&
xopt1 = fmin(myfunc, x0)
print xopt1
xopt2 = fmin_powell(myfunc, x0)
print xopt2
xopt3 = fmin_bfgs(myfunc, x0)
print xopt3
xopt4 = fmin_cg(myfunc,x0)
print xopt4
给出的结果是：
Optimization terminated
successfully.
Current function value: 4.000000
Iterations: 16
Function evaluations: 32
Optimization terminated
successfully.
Current function value: 4.000000
Iterations: 2
Function evaluations: 53
Optimization terminated
successfully.
Current function value: 4.000000
Iterations: 2
Function evaluations: 12
Gradient evaluations: 4
Optimization terminated
successfully.
Current function value: 4.000000
Iterations: 2
Function evaluations: 15
Gradient evaluations: 5
我们可以根据给出的消息直观的判断算法的执行情况。每一种算法数学上的问题，请自己看书学习。个人感觉，如果不是纯研究数学的工作，没必要搞清楚那些推导以及定理云云。不过，必须了解每一种算法的优劣以及能力所及。在使用的时候，不妨多种算法都使用一下，看看效果分别如何，同时，还可以互相印证算法失效的问题。
在from scipy.optimize import
fmin之后，就可以使用help(fmin)来查看fmin的帮助信息了。帮助信息中没有例子，但是给出了每一个参数的含义说明，这是调用函数时候的最有价值参考。
有源码研究癖好的，或者当你需要改进这些已经实现的算法的时候，可能需要查看optimize中的每种算法的源代码。在这里：
聪明的你肯定发现了，顺着这个链接往上一级、再往上一级，你会找到scipy的几乎所有源码！
我的更多文章：
( 21:26:05)( 19:46:30)( 13:03:50)( 13:41:19)( 16:45:22)
已投稿到：
以上网友发言只代表其个人观点，不代表新浪网的观点或立场。}