Julia统计计算介绍¶
作者：李东风
日期：2020-04-05
Julia比较适合用作数值计算，
编程既有Python、R、Matlab这样的语言的简洁，
又有C++这样的编译语言的运行效率。
统计数据分析、作图需要用到许多复杂的算法，
有些算法耗时很多，
比如MCMC等。
大量数据的分析、计算、测试都需要易用的编程和高效的运行效率，
Julia在这两点都很适合。
Julia用作统计数据分析，
缺点是其问世时间还比较短，
许多统计模型还没有现成的软件包。
但是，
基本的数据管理、作图、统计分布、随机模拟、最优化等统计计算必须的功能已经很完善，
所以能够自己编写的统计计算任务很适合用Julia编写。
Julia目前的版本是1.4.0。
本文承接“Julia语言入门”和“Julia统计数据分析入门”，
从自己编写程序作统计计算的角度简单介绍Julia中与统计计算编程有关的功能，
如向量、矩阵计算，
最优化，
随机模拟，
并行计算等。
向量和矩阵计算¶
Julia的统计功能还不够丰富，
许多统计任务可以通过R语言来完成。
Julia的优点是其运算速度快，
特别是在对向量和矩阵循环时效率很高，
这正是R语言和Python语言的弱项。
所以不需要调用很多现成的统计方法软件包的统计计算编程任务可以用Julia语言完成。
这里讲一些统计计算中常用的编程技巧，
统计计算中常用的一些功能，
提高运行效率的方法等。
首先介绍一些向量和矩阵计算有关知识。
向量和矩阵计算¶
向量是一维数组，
矩阵是二维数组。
数组元素之间的四则运算，
只要使用加点格式，
如x ./ y表示两个数组对应元素之间的除法。
如果函数f(x)可以用于标量，
则f.(x)可以用于数组每个元素。
但是，目前用显式循环而不是向量化写法的程序效率更高。
向量和矩阵运算的许多函数由标准库中的LinearAlgebra包提供。
两个向量x和y的内积用dot(x,y)计算。
两个矩阵A和B的矩阵乘法用A * B表示，
矩阵A的共轭转置用A'表示，
转置用A.'表示。
[A B C]表示三个矩阵横向合并，
[A;B;C]表示三个矩阵纵向合并。
若x是向量，
diagm(x)返回以x为主对角线元素的对角阵，
而Diagonal(x)返回带有特殊标志的对角阵。
inv(A)表示$A^{-1}$。
A \ B表示求解线性方程组A x = B。
当$A$的行数超过列数时，
求最小二乘解。
当$A$的列数超过行数时，
可能有无穷多解，
这时返回其中长度最小的解。
统计计算中用到许多矩阵计算，
高效、精确的矩阵计算常常需要通过矩阵的三角分解、正交-三角分解等实现。
当$A$为方阵时，
lu(A)将矩阵$A$分解成$A = PLU$形式，
其中$P$是置换阵(permutation matrix)，
$L$是下三角阵，
$U$是上三角阵。
设结果为res，
则三个部分分别用res.P, res.L, res.U访问。
置换阵res.P还可以用一维数组res.p表示。
如
In [49]:
A1 = [1 -2 2; 1 -1 2; -1 1 1]
Out[49]:
3×3 Array{Int64,2}:
1 -2 2
1 -1 2
-1 1 1
In [50]:
using LinearAlgebra
res1 = LinearAlgebra.lu(A1)
res1.L
Out[50]:
3×3 Array{Float64,2}:
1.0 0.0 0.0
1.0 1.0 0.0
-1.0 -1.0 1.0
In [51]:
res1.U
Out[51]:
3×3 Array{Float64,2}:
1.0 -2.0 2.0
0.0 1.0 0.0
0.0 0.0 3.0
In [52]:
res1.P
Out[52]:
3×3 Array{Float64,2}:
1.0 0.0 0.0
0.0 1.0 0.0
0.0 0.0 1.0
In [53]:
res1.p
Out[53]:
3-element Array{Int64,1}:
1
2
3
对矩阵$A$，
qr(A)可以做QR分解$A=QR$，
qr(A, true)可以做带枢轴量的QR分解$A=PQR$，
如果结果为res，
则各个部分用res.P, res.Q, res.R访问，
置换阵res.P还可以用res.p，一个一维向量表示。
用eigen(A)可以求方阵A的特征值分解，
如果结果为res，则res.values为特征值，是复数向量；
res.vectors返回各列为特征向量的矩阵，
inv(res), det(res), isposdef(res)可以从分解结果中提取信息。
如
In [54]:
res2 = eigen(A1)
res2.values
Out[54]:
3-element Array{Complex{Float64},1}:
-0.2873715369435107 - 1.3499963980036567im
-0.2873715369435107 + 1.3499963980036567im
1.5747430738870216 + 0.0im
In [55]:
res2.vectors
Out[55]:
3×3 Array{Complex{Float64},2}:
0.783249-0.0im 0.783249+0.0im 0.237883+0.0im
0.493483+0.303862im 0.493483-0.303862im 0.651777+0.0im
-0.0106833-0.22483im -0.0106833+0.22483im 0.720138+0.0im
对称阵(或厄米特阵)的特征值都是实数，
其计算更快。
因为数值误差的关系，
有些对称阵可能保存起来已经不对称，
这时可以用Symmetric(A)将A转换为对称阵，
用Hermitian(A)将A转换为厄米特阵。
向量化与循环¶
Julia中，
显式循环一般比用向量化的计算效率更高。
这是因为向量化运算需要考虑到动态类型问题，
而显示循环则直接访问元素。
考虑下面的例子：
In [56]:
function vecadd1(a,b,c,N)
for i = 1:N
c = a + b
end
end
function vecadd2(a,b,c,N)
for i = 1:N
c .= a .+ b
end
end
function vecadd3(a,b,c,N)
for i = 1:N, j = 1:length(c)
c[j] = a[j] + b[j]
end
end
Out[56]:
vecadd3 (generic function with 1 method)
以上的三个函数都是将向量a与向量b相加赋给向量c，
并重复N次。
第一个版本用了向量化的做法，
第二个版本用了加点格式的向量化做法，
第三个版本是显式循环。
下面比较这三个版本的运行效率：
In [57]:
A = rand(2); B = rand(2); C = zeros(2);
@elapsed vecadd1(A,B,C,10);
@elapsed vecadd2(A,B,C,10);
@elapsed vecadd3(A,B,C,10);
In [58]:
@elapsed vecadd1(A,B,C,10_000_000)
Out[58]:
1.228392501
In [59]:
@elapsed vecadd2(A,B,C,10_000_000)
Out[59]:
0.0907257
In [60]:
@elapsed vecadd3(A,B,C,10_000_000)
Out[60]:
0.035474
在上述演示程序中，
为了避免函数编译的影响，
先用较小的N将三个函数运行了一遍。
可以看出，
第一个版本效率最差，
这是因为计算a + b时需要动态分配结果的内存，
赋值给c也不是存入到c的已有内存而是重新绑定c。
第二个版本用了加点形式，
比第一个版本效率要高，
但是也有动态数据类型的开销。
第三个版本则完全是在访问a, b, c已有内存，
没有做任何的内存重分配或者重绑定，
所以效率最高。
随机数发生器与统计分布¶
随机数发生器¶
在统计建模中经常需要用到各种概率分布，
Julia的Distributions包提供了常见的概率分布的有关的函数，
包括各种分布的随机数发生器。
Julia语言的标准库Random包提供了rand函数和randn等函数。
rand(n)返回n个标准均匀分布U(0,1)随机数，
类型是Float64的一维数组，
无自变量的rand()调用返回一个U(0,1)随机数。
randn(n)和randn()产生标准正态分布的随机数。
Random.seed!(rng, seed)可以设置种子。
为了能够保证模拟的结果可重复，
需要在模拟开始之前设置种子。
但是，随着Julia版本的升级，
随机数种子和随机数生成器可能会改变，
如果特别需要模拟研究结果完全可重复，
可以将生成的模拟样本都备份保存下来。
加载Distributions包后可以生成其它分布的随机数，
为Base的rand函数增加了方法，
还提供了各种与分布的有关函数。
产生随机数的函数调用格式为
rand(distribution, n)或rand(distribution)，
其中distribution是一个分布，用函数调用表示，
比如标准正态分布用Normal()表示，$\text{N}(10, 2^2)$用Normal(10,2)表示。
用Random.seed!(rng, k)设定一个随机数种子序号，
这样可以在模拟时获得可重复的结果。
比如，生成100个Possion(2)分布随机数:
In [61]:
using Random, Distributions
Random.seed!(101)
y = rand(Poisson(2), 100);
结果是Int64型的一维数组。样本的频数统计：
In [62]:
using StatsBase
using DataFrames
tab11 = by(DataFrame(y=y), :y, df -> DataFrame(count=size(df, 1)))
sort!(tab11)
Out[62]:
ycount
Int64Int647 rows × 2 columns
1021
2130
3231
4311
545
651
761
作其离散分布条形图图：
In [63]:
using Plots; Plots.gr()
using StatsPlots
@df tab11 Plots.bar(:y, :count,
title="Poisson(2) Variates",
xticks=0:7,
legend=false)
UndefVarError: select not defined
Stacktrace:
[1] extract_columns_and_names(::DataFrame, ::Symbol, ::Vararg{Symbol,N} where N) at C:\Users\user\.juliapro\JuliaPro_v1.4.0-1\packages\StatsPlots\Y1rGZ\src\df.jl:159
[2] (::var"#7#8")(::DataFrame) at .\none:0
[3] top-level scope at In[63]:3
如下程序生成1000个N($10, 2^2$)随机数，
并作分布直方图：
In [64]:
Random.seed!(101)
y = rand(Normal(10, 2), 1000);
Plots.histogram(y, nbins=30, legend=false)
Out[64]:
3
6
9
12
15
0
25
50
75
100
随机抽样¶
StatsBase.sample(x,n)函数可以从一列值x中随机有放回抽取n个值。
加replace=false作无放回抽样。
如
In [65]:
StatsBase.sample(1:10, 5, replace=false)
Out[65]:
5-element Array{Int64,1}:
4
10
8
7
9
在社会调查抽样中，
经常需要按不等权重抽样。
这时，可以增加一个权重向量。
权重向量的数据类型是AbstractWeights,
不需要总和等于1，
其总和会保存在变量中。
用Weights()函数将普通向量转换成权重向量。
如
In [66]:
StatsBase.sample(1:10, StatsBase.Weights([ones(5); 4 .* ones(5)]), 5; replace=false)
Out[66]:
5-element Array{Int64,1}:
8
5
7
1
9
与分布有关的函数¶
Distributions包不仅包含了密度函数、对数密度函数、分布函数、分位数函数、随机数函数等与分布有关的函数，
还可以输出每个参数分布的各种矩的理论值，
计算矩母函数和特征函数，
给定观测样本后计算参数的最大似然估计，
用共轭先验计算后验分布，对参数作最大后验密度估计等。
比R提供了更多的分布类型和更多的函数类型。
分布类型包括一元随机变量(Univariate)、随机向量(Multivariate)、随机矩阵(Matrixvariate)，
又可分为离散型(Discrete)和连续型(Continuous)。
分布类型的大类名称包括：
DiscreteUnivariateDistribution
ContinuousUnivariateDistribution
DiscreteMultivariateDistribution
ContinuousMultivariateDistribution
DiscreteMatrixDistribution
ContinuousMatrixDistribution
Distributions包中定义了各种概率分布。
例如，
Binomial(n, p)表示二项分布，
Normal(μ, σ)表示正态分布分布，
Multinomial(n, p)表示多项分布，
Wishart(nu, S)表示Wishart分布。
可以从一元分布生成截断分布(truncated)，
如Truncated(Normal(mu, sigma), l, u)。
用params(distribution)返回一个分布的参数，如
In [67]:
params(Binomial(10, 0.2))
Out[67]:
(10, 0.2)
rand(distribution, n)生成n个指定分布的随机数，
对一元分布，结果是一维数组，
连续型分布的元素类型是Float64,
离散型分布的元素类型是Int。
对多元分布，
结果是矩阵，每列为一个抽样。
对随机矩阵分布，
结果是元素为矩阵的数组。
有一些函数用来返回分布的特征量，
对一元分布，有：
maximum() 分布的定义域上界
minimum() 分布的定义域下界
mean() 期望值
var() 方差
std() 标准差,
median() 中位数
modes() 众数(可有多个)
mode() 第一个众数
skewness() 偏度
kurtosis() 峰度，正态为0
entropy() 分布的熵
用mgf(distribution, t)计算矩母函数在t处的值，
用cf(distribution, t)计算特征函数在t处的值。
用pdf(distribution, x)计算分布密度(概率质量)函数在x处的值，
如果x是数组，结果返回数组；
如果x是数组且函数写成pdf!(distribution, x)，则结果保存在x中返回。
logpdf()计算其密度函数的自然对数值。
用cdf(distribution, x)计算分布函数在x处的值，
logcdf计算其对数值，
ccdf计算1减分布函数值，
logccdf计算其对数值。
用quantile(distribution, p)计算分位数函数在p处的值。
类似于pdf，这些函数一般支持加点格式的向量化，
输入x为向量时返回向量，
写成以叹号结尾的函数名时，
输入为向量则结果保存在输入向量的存储中返回。
可这样向量化的函数包括
pdf, logpdf, cdf, logcdf, ccdf, logccdf, quantile, cquantile, invlogcdf, invlogccdf。
比如，计算$\text{N}(10,2^2)$的0.5和0.975分位数：
In [68]:
quantile.(Normal(10, 2), [0.5, 0.975])
Out[68]:
2-element Array{Float64,1}:
10.0
13.919927969080115
In [69]:
10 .+ 2 .* quantile.(Normal(), [0.5, 0.975])
Out[69]:
2-element Array{Float64,1}:
10.0
13.919927969080115
比如，对标准正态分布密度作曲线图：
In [70]:
using Plots; Plots.gr()
Plots.plot(x -> pdf(Normal(), x), -3, 3, legend=false)
Out[70]:
-3
-2
-1
0
1
2
3
0.0
0.1
0.2
0.3
0.4
用rand(distrbution)生成指定分布的一个抽样，
用rand(distrbution, n)生成指定分布的n个抽样，
用rand!(distribution, x)生成指定分布的多个抽样保存在数组x中返回。
用fit(distributionname, x)作参数估计，
一般都是最大似然估计，也可以用fit_mle()规定使用最大似然估计。
比如，
模拟生成$\text{N}(10,2^2)$样本，然后做参数估计：
In [71]:
Random.seed!(101); x=rand(Normal(10, 2), 30)
fit(Normal, x)
Out[71]:
Normal{Float64}(μ=10.349862945450266, σ=2.4124946578030015)
In [72]:
fit_mle(Normal, x)
Out[72]:
Normal{Float64}(μ=10.349862945450266, σ=2.4124946578030015)
支持的分布¶
这里列举出Julia的Distributions包支持的分布。
连续型一元分布¶
Arcsine(a,b)
Beta(α,β)
BetaPrime(α,β)
Biweight(μ, σ)
Chi(ν)
Chisq(ν)
Cosine(μ, σ)
Epanechnikov(μ, σ)
Erlang(α,θ)
Exponential(θ)
FDist(ν1, ν2)
Frechet(α,θ)
Gamma(α,θ)
GeneralizedExtremeValue(μ, σ, ξ)
GeneralizedPareto(μ, σ, ξ)
Gumbel(μ, θ)
InverseGamma(α, θ)
InverseGaussian(μ,λ)
Kolmogorov()
KSDist(n)
KSOneSided(n)
Laplace(μ,θ)
Levy(μ, σ)
Logistic(μ,θ)
LogNormal(μ,σ)
NoncentralBeta(α, β, λ)
NoncentralChisq(ν, λ)
NoncentralF(ν1, ν2, λ)
NoncentralT(ν, λ)
Normal(μ,σ)
NormalCanon(η, λ)
NormalInverseGaussian(μ,α,β,δ)
Pareto(α,θ)
Rayleigh(σ)
SymTriangularDist(μ,σ)
TDist(ν)
TriangularDist(a,b,c)
Triweight(μ, σ)
Uniform(a,b)
VonMises(μ, κ)
Weibull(α,θ)
离散型一元分布¶
Bernoulli(p)
BetaBinomial(n,α,β)
Binomial(n,p)
Categorical(p)
DiscreteUniform(a,b)
Geometric(p)
Hypergeometric(s, f, n)
NegativeBinomial(r,p)
Poisson(λ)
PoissonBinomial(p)
Skellam(μ1, μ2)
截断分布¶
设distribution为一个一元分布，
用
Truncated(distribution, l, u)
可以返回一个截断分布，其中l为下界，
u为上界。
特别地，
TruncatedNormal(mu, sigma, l, u)表示截断正态分布。
多元分布¶
Multinomial
MvNormal
MvNormalCanon
MvLogNormal
Dirichlet
随机矩阵分布¶
Wishart(nu, S)
InverseWishart(nu, P)
混合分布¶
通过MixtureModel()函数构造。
随机模拟软件¶
对于离散事件系统模拟，
比如排队系统的模拟，
SimJulia包可以很容易地建模进行模拟。
一般的独立样本抽样、重要性抽样等可以直接利用Julia的随机数发生器模拟。
贝叶斯推断问题一般使用MCMC，
Julia的Stan包调用Stan程序库进行MCMC计算，
Stan程序会将模型转换为C++然后编译成本地二进制代码运行。
Julia的Mamba包是纯Julia语言进行MCMC模拟的扩展包，
这样模型的描述可以利用Julia程序的所有功能。
最优化与方程求根¶
统计计算的常用工具包括矩阵计算、随机模拟和最优化，
统计计算中许多问题都归结为一个最优化问题，
比如非线性回归、最大似然估计、惩罚最大似然估计等。
Julia提供了若干个优秀的最优化计算扩展包，
比如JuMP、Optim、Convex等扩展包。
Optim包¶
Optim包的好处是完全用Julia实现，
不像其它的最优化包那样需要依赖于其它的可能需要收费的、需要其它编译器的后端引擎。
这部分内容采用了Optim包文档中的例子。
多元函数优化¶
Optim包的optimize()函数类似于R的optim()，
可以用来求多元函数的无约束极值或矩形约束的极值，
一元函数区间内的极值，
复函数极值,
流形上的极值。
可用的优化方法包括:
NelderMead()，是不依赖梯度信息的单纯型方法，只需要目标函数本身；
BFGS(), LBFGS(), 这是著名的拟牛顿法, LBFGS()是内存开销有限版本;
ConjugateGradient(), 共轭梯度法；
GradientDescent(), 最速下降法，还有MomentumGradientDescent(), AcceleratedGradientDescent()；
Newton()，牛顿法，需要梯度函数与海色阵函数，海色阵函数不可省略，
还有NewtonTrustRegion()，此方法克服了当函数不能用二次多项式近似时的困难；
SimulatedAnnealing()，模拟退火法，用于求解全局最值，
只需要目标函数本身，速度慢。
SAMIN(), ParticleSwarm()是类似方法;
Fminbox()，矩形约束的优化；
IPNewton()为非线性约束的内点法加牛顿法；
Brent()和GoldenSection()为一元函数优化；
…………
用户可以提供梯度函数以及海色阵用来加速收敛。
结果是类似R的list类型的一个结构。
举例说明。
考虑二元函数
$f(x, y) = (x - a)^2 + b (y - x^2)^2$的最小值问题。
显然唯一的最小值点为$(a, a^2)$。
用Optim的optimize()函数进行数值求解。
这个二元函数的最小值点很长、很窄、底部平坦的山谷内，
对迭代算法来说有一定困难。
取$a=1, b=100$，
最小值点在$(1,1)$。
取初值为$(0,0)$。
目标函数定义：
In [73]:
function rb(x::Vector{Float64})
return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2;
end
Out[73]:
rb (generic function with 1 method)
使用默认的设置进行优化：
In [74]:
using Optim
ores1 = optimize(rb, [0.0, 0.0])
Out[74]:
* Status: success
* Candidate solution
Minimizer: [1.00e+00, 1.00e+00]
Minimum: 3.525527e-09
* Found with
Algorithm: Nelder-Mead
Initial Point: [0.00e+00, 0.00e+00]
* Convergence measures
√(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 60
f(x) calls: 118
可以看出在没有选用具体优化方法的情况下用了Nelder-Mead方法，
迭代了60次，找到的最小值点的误差在$10^{-4}$以下。
改用BFGS方法，使用数值微分计算梯度：
In [75]:
ores1b = optimize(rb, [0.0, 0.0], BFGS())
Out[75]:
* Status: success
* Candidate solution
Minimizer: [1.00e+00, 1.00e+00]
Minimum: 5.471433e-17
* Found with
Algorithm: BFGS
Initial Point: [0.00e+00, 0.00e+00]
* Convergence measures
|x - x'| = 3.47e-07 ≰ 0.0e+00
|x - x'|/|x'| = 3.47e-07 ≰ 0.0e+00
|f(x) - f(x')| = 6.59e-14 ≰ 0.0e+00
|f(x) - f(x')|/|f(x')| = 1.20e+03 ≰ 0.0e+00
|g(x)| = 2.33e-09 ≤ 1.0e-08
* Work counters
Seconds run: 1 (vs limit Inf)
Iterations: 16
f(x) calls: 53
∇f(x) calls: 53
只迭代了16次，但是计算数值微分也会涉及到比较多的目标函数计算。
可以人为给出梯度函数：
In [76]:
function rb_grad!(G::Vector{Float64}, x::Vector{Float64})
G[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
G[2] = 200.0 * (x[2] - x[1]^2)
end
ores1c = optimize(rb, rb_grad!, [0.0, 0.0], BFGS())
Out[76]:
* Status: success
* Candidate solution
Minimizer: [1.00e+00, 1.00e+00]
Minimum: 7.645563e-21
* Found with
Algorithm: BFGS
Initial Point: [0.00e+00, 0.00e+00]
* Convergence measures
|x - x'| = 3.48e-07 ≰ 0.0e+00
|x - x'|/|x'| = 3.48e-07 ≰ 0.0e+00
|f(x) - f(x')| = 6.91e-14 ≰ 0.0e+00
|f(x) - f(x')|/|f(x')| = 9.03e+06 ≰ 0.0e+00
|g(x)| = 2.32e-09 ≤ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 16
f(x) calls: 53
∇f(x) calls: 53
还可以由用户提供一个海色阵函数，
这时方法为牛顿法，调用如
optimize(rb, rb_grad!, rb_Hess!, [0.0, 0.0])
optimize()不会用数值微分法计算海色阵，
因为两次数值微分造成很大误差。
如果目标函数完全用Julia编写并且不调用C代码和Fortran代码，
则目标函数的梯度函数和海色阵函数很可能可以用符号计算方法获得精确函数。
将上述rb()函数去掉其中的声明才能做符号计算。
如
In [77]:
function rbsym(x)
return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2;
end
optimize(rbsym, [0.0, 0.0], Newton(); autodiff = :forward)
Out[77]:
* Status: success
* Candidate solution
Minimizer: [1.00e+00, 1.00e+00]
Minimum: 3.081488e-31
* Found with
Algorithm: Newton's Method
Initial Point: [0.00e+00, 0.00e+00]
* Convergence measures
|x - x'| = 3.06e-09 ≰ 0.0e+00
|x - x'|/|x'| = 3.06e-09 ≰ 0.0e+00
|f(x) - f(x')| = 9.34e-18 ≰ 0.0e+00
|f(x) - f(x')|/|f(x')| = 3.03e+13 ≰ 0.0e+00
|g(x)| = 1.11e-15 ≤ 1.0e-08
* Work counters
Seconds run: 1 (vs limit Inf)
Iterations: 14
f(x) calls: 44
∇f(x) calls: 44
∇²f(x) calls: 14
optimize()支持比较简单的矩形约束优化，
所谓矩形约束(box constraint)，
是指每个坐标分量属于一个区间，
可以用Inf和-Inf表示无穷区间。
这种约束优化用的是简单的内点法
例如，上面的例子函数本来是无约束的，
但是如果我们要求解在矩形$[1.25, \infty) \times [-2.1, \infty)$中，
这个范围已经不包含全局最小值点。
约束优化程序如下：
In [78]:
lower = [1.25, -2.1]
upper = [Inf, Inf]
initial_x = [2.0, 2.0]
inner_optimizer = GradientDescent()
ores2 = optimize(rb, rb_grad!, lower, upper, initial_x,
Fminbox(inner_optimizer))
Out[78]:
* Status: success
* Candidate solution
Minimizer: [1.25e+00, 1.56e+00]
Minimum: 6.250000e-02
* Found with
Algorithm: Fminbox with Gradient Descent
Initial Point: [2.00e+00, 2.00e+00]
* Convergence measures
|x - x'| = 0.00e+00 ≤ 0.0e+00
|x - x'|/|x'| = 0.00e+00 ≤ 0.0e+00
|f(x) - f(x')| = 0.00e+00 ≤ 0.0e+00
|f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
|g(x)| = 5.00e-01 ≰ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 8
f(x) calls: 84562
∇f(x) calls: 84562
可见结果在约束的左边界处达到。
optimize()的结果如果保存在res变量中，
可以用Optim.minimizer(res)获得最小值点(总是一个向量)，
用Optim.minimum(res)获得最小值。
结果中还有用于检查是否收敛的输出，
详见Optim包的文档。
控制迭代精度等选项，
可以使用optimize()的控制选项，如
In [79]:
optimize(rb, rb_grad!, [0.0, 0.0],
GradientDescent(),
Optim.Options(g_tol = 1e-12,
iterations = 10,
store_trace = true,
show_trace = false))
Out[79]:
* Status: failure (reached maximum number of iterations) (line search failed)
* Candidate solution
Minimizer: [3.73e-01, 1.34e-01]
Minimum: 3.959488e-01
* Found with
Algorithm: Gradient Descent
Initial Point: [0.00e+00, 0.00e+00]
* Convergence measures
|x - x'| = 7.94e-03 ≰ 0.0e+00
|x - x'|/|x'| = 2.13e-02 ≰ 0.0e+00
|f(x) - f(x')| = 7.27e-03 ≰ 0.0e+00
|f(x) - f(x')|/|f(x')| = 1.84e-02 ≰ 0.0e+00
|g(x)| = 1.16e+00 ≰ 1.0e-12
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 10
f(x) calls: 34
∇f(x) calls: 34
此例因为设置了太小的迭代次数限制导致没有收敛。
实际上，
最速下降法迭代一万次也没有收敛。
这是因为最速下降法收敛速度慢，
而且目标函数也比较复杂。
一元函数优化¶
optimize()也可以对一个一元函数在给定的区间内求极值，
可用方法有Brent()和GoldenSection()。
比如函数$f(x) = 2x^2+3x+1 = 2 (x + \frac34)^2 - \frac18$，
显然最小值点为$-\frac34$，
最小值为$-\frac18$。
在区间$[-10, 10]$内求最小值：
In [80]:
funi(x) = 2x^2 + 3x + 1
optimize(funi, -10, 10)
Out[80]:
Results of Optimization Algorithm
* Algorithm: Brent's Method
* Search Interval: [-10.000000, 10.000000]
* Minimizer: -7.500000e-01
* Minimum: -1.250000e-01
* Iterations: 5
* Convergence: max(|x - x_upper|, |x - x_lower|) <= 2*(1.5e-08*|x|+2.2e-16): true
* Objective Function Calls: 6
JuMP包¶
Optim包好处是没有对其它库的依赖，
缺点是能解决的优化问题的类型很少，
主要是无约束优化与最简单的约束优化。
JuMP包是Julia中优化问题建模的一个软件包，
调用一些开源或者商业化的最优化程序包的功能来求解，
所以能解决的优化问题类型很多，
但是其安装设置比较困难，
有些后端软件是收费软件。
涉及的优化问题类型包括线性规划，
(混合)整数规划，
二阶锥规划，
半定规划，
非线性规划。
JuMP可以看成是一个前端工具，
为了使用JuMP包，
需要安装JuMP包，
还要单独安装所需要用的求解包，
这些包有的需要调用二进制代码库，
好在Windows和OS X上有的库已经提供了编译好的包下载，
Linux可以从源程序编译。
比如，GLPK提供了线性规划(LP)和混合整数线性规划(MILP)功能，
用
Pkg.add("GLPK")
安装。
在Windows的Julia中安装需要安装许多辅助库，
因为国内网络问题经常失败，
可以在网络畅通时比如清晨安装。
线性规划¶
这里举例说明用JuMP调用GLPK求解线性规划问题。
例子来自JuMP的文档。
考虑如下的线性规划问题：
$$
\begin{aligned}
& \max 5x + 3y \quad \text{s.t.}\\
& 0 \leq x \leq 2 \\
& 0 \leq y \leq 30 \\
& x + 5y \leq 3.0
\end{aligned}
$$
调用JuMP包和GLPK包：
In [81]:
using JuMP
using GLPK
用Model()函数生成一个优化问题模型，
其中with_optimizer()中选择后端的优化库：
In [82]:
omod21 = Model(with_optimizer(GLPK.Optimizer))
Out[82]:
$$ \begin{alignat*}{1}\text{feasibility}\\
\text{Subject to} \quad\end{alignat*}
$$
下面指定自变量以及自变量取值区间：
In [83]:
@variable(omod21, 0 <= x <= 2 );
@variable(omod21, 0 <= y <= 30 );
指定目标函数以及要求最大化还是最小化：
In [84]:
@objective(omod21, Max, 5x + 3*y )
Out[84]:
$$ 5 x + 3 y $$
注意Julia中5*x可以写成5x但是不要写成5 x。
指定一个线性约束：
In [85]:
@constraint(omod21, con, 1x + 5y <= 3.0 )
Out[85]:
con : $ x + 5 y \leq 3.0 $
显示目前的模型：
In [86]:
print(omod21)
Max 5 x + 3 y
Subject to
x + 5 y <= 3.0
x >= 0.0
y >= 0.0
x <= 2.0
y <= 30.0
以上完成了对一个优化问题的描述，
保存在变量omod21中。
用JuMP.optimize!(model)函数调用后端的优化程序库求解：
In [87]:
JuMP.optimize!(omod21)
求解结束后，结果也保存在模型变量omod21中。
求解结束可能处于如下状态：
求出最优解，或者证明无解；
因数值计算问题不能继续，或者迭代次数超界或者运算时间超界而停止。
用JuMP.termination_status(model)查询求解的结束原因：
In [88]:
println(termination_status(omod21))
OPTIMAL
这个状态的表示找到了最优解。
对线性规划问题，
还可以检查是否找到了primal-dual对，
即可行解集与对偶的零对偶间隙(zero duality gap):
In [89]:
println(primal_status(omod21))
println(dual_status(omod21))
FEASIBLE_POINT
FEASIBLE_POINT
上面两个结果表示primal-dual对都是可行的。
因为求解停止状态是OPTIMAL(最优)，
所以解是最优。
显示最优的目标函数值：
In [90]:
println("最大目标函数值：", objective_value(omod21))
最大目标函数值：10.6
显示最大值点，保存在模型指定的自变量名：
In [91]:
println("最大值点：x=", value(x), " y=", value(y))
最大值点：x=2.0 y=0.2
非线性规划¶
这里给一个利用JuMP调用Ipopt后端的例子。
Ipopt是一个开源的C++非线性最优化库。
考虑二元函数$(1 - x)^2 + 100 * (y - x^2)^2$的无约束最小值问题。
程序如下：
In [92]:
using JuMP
using Ipopt
建立模型：
In [93]:
omod22 = Model(with_optimizer(Ipopt.Optimizer))
Out[93]:
$$ \begin{alignat*}{1}\text{feasibility}\\
\text{Subject to} \quad\end{alignat*}
$$
定义变量：
In [94]:
@variable(omod22, x, start = 0.0);
@variable(omod22, y, start = 0.0);
定义目标函数，选择最小化：
In [95]:
@NLobjective(omod22, Min, (1 - x)^2 + 100 * (y - x^2)^2)
求解：
In [96]:
optimize!(omod22)
This is Ipopt version 3.12.10, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).
Number of nonzeros in equality constraint Jacobian...: 0
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 3
Total number of variables............................: 2
variables with only lower bounds: 0
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 0
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.0000000e+00 0.00e+00 2.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 9.5312500e-01 0.00e+00 1.25e+01 -1.0 1.00e+00 - 1.00e+00 2.50e-01f 3
2 4.8320569e-01 0.00e+00 1.01e+00 -1.0 9.03e-02 - 1.00e+00 1.00e+00f 1
3 4.5708829e-01 0.00e+00 9.53e+00 -1.0 4.29e-01 - 1.00e+00 5.00e-01f 2
4 1.8894205e-01 0.00e+00 4.15e-01 -1.0 9.51e-02 - 1.00e+00 1.00e+00f 1
5 1.3918726e-01 0.00e+00 6.51e+00 -1.7 3.49e-01 - 1.00e+00 5.00e-01f 2
6 5.4940990e-02 0.00e+00 4.51e-01 -1.7 9.29e-02 - 1.00e+00 1.00e+00f 1
7 2.9144630e-02 0.00e+00 2.27e+00 -1.7 2.49e-01 - 1.00e+00 5.00e-01f 2
8 9.8586451e-03 0.00e+00 1.15e+00 -1.7 1.10e-01 - 1.00e+00 1.00e+00f 1
9 2.3237475e-03 0.00e+00 1.00e+00 -1.7 1.00e-01 - 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 2.3797236e-04 0.00e+00 2.19e-01 -1.7 5.09e-02 - 1.00e+00 1.00e+00f 1
11 4.9267371e-06 0.00e+00 5.95e-02 -1.7 2.53e-02 - 1.00e+00 1.00e+00f 1
12 2.8189505e-09 0.00e+00 8.31e-04 -2.5 3.20e-03 - 1.00e+00 1.00e+00f 1
13 1.0095040e-15 0.00e+00 8.68e-07 -5.7 9.78e-05 - 1.00e+00 1.00e+00f 1
14 1.3288608e-28 0.00e+00 2.02e-13 -8.6 4.65e-08 - 1.00e+00 1.00e+00f 1
Number of Iterations....: 14
(scaled) (unscaled)
Objective...............: 1.3288608467480825e-28 1.3288608467480825e-28
Dual infeasibility......: 2.0183854587685121e-13 2.0183854587685121e-13
Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 2.0183854587685121e-13 2.0183854587685121e-13
Number of objective function evaluations = 36
Number of objective gradient evaluations = 15
Number of equality constraint evaluations = 0
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 0
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 14
Total CPU secs in IPOPT (w/o function evaluations) = 0.216
Total CPU secs in NLP function evaluations = 0.125
EXIT: Optimal Solution Found.
上面的显示已经表明找到了最优解。
用termination_status(omod21)查询算法退出状态：
In [97]:
println(termination_status(omod22))
LOCALLY_SOLVED
结果显示得到局部最优，
因为大多数算法只能保证局部最优，
所以没有问题。
显示目标函数值和最小值点：
In [98]:
println("最小化目标函数值=", objective_value(omod22))
println("最小值点：x = ", value(x), " y = ", value(y))
最小化目标函数值=1.3288608467480825e-28
最小值点：x = 0.9999999999999899 y = 0.9999999999999792
增加一个简单的线性约束$x+y=10$，
重新求解：
In [99]:
@constraint(omod22, x + y == 10)
optimize!(omod22);
This is Ipopt version 3.12.10, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).
Number of nonzeros in equality constraint Jacobian...: 2
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 3
Total number of variables............................: 2
variables with only lower bounds: 0
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 1
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.0000000e+00 1.00e+01 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 9.6315968e+05 0.00e+00 3.89e+05 -1.0 9.91e+00 - 1.00e+00 1.00e+00h 1
2 1.6901461e+05 0.00e+00 1.16e+05 -1.0 3.24e+00 - 1.00e+00 1.00e+00f 1
3 2.5433173e+04 1.78e-15 3.18e+04 -1.0 2.05e+00 - 1.00e+00 1.00e+00f 1
4 2.6527756e+03 0.00e+00 7.79e+03 -1.0 1.19e+00 - 1.00e+00 1.00e+00f 1
5 1.1380324e+02 0.00e+00 1.35e+03 -1.0 5.62e-01 - 1.00e+00 1.00e+00f 1
6 3.3745506e+00 0.00e+00 8.45e+01 -1.0 1.50e-01 - 1.00e+00 1.00e+00f 1
7 2.8946196e+00 0.00e+00 4.22e-01 -1.0 1.07e-02 - 1.00e+00 1.00e+00f 1
8 2.8946076e+00 0.00e+00 1.07e-05 -1.7 5.42e-05 - 1.00e+00 1.00e+00f 1
9 2.8946076e+00 0.00e+00 5.91e-13 -8.6 1.38e-09 - 1.00e+00 1.00e+00f 1
Number of Iterations....: 9
(scaled) (unscaled)
Objective...............: 2.8946075504894599e+00 2.8946075504894599e+00
Dual infeasibility......: 5.9130478291535837e-13 5.9130478291535837e-13
Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 5.9130478291535837e-13 5.9130478291535837e-13
Number of objective function evaluations = 10
Number of objective gradient evaluations = 10
Number of equality constraint evaluations = 10
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 1
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 9
Total CPU secs in IPOPT (w/o function evaluations) = 0.005
Total CPU secs in NLP function evaluations = 0.000
EXIT: Optimal Solution Found.
In [100]:
println(termination_status(omod22));
println("最小化目标函数值=", objective_value(omod22))
println("最小值点：x = ", value(x), " y = ", value(y))
LOCALLY_SOLVED
最小化目标函数值=2.89460755048946
最小值点：x = 2.701147124098218 y = 7.2988528759017814
考虑目标函数在简单约束$x \geq 1.25$, $y \geq -2.1$下的优化。
In [101]:
using JuMP
using Ipopt
omod23 = Model(with_optimizer(Ipopt.Optimizer));
@variable(omod23, x, start = 0.0);
@variable(omod23, y, start = 0.0);
@constraint(omod23, x >= 1.25);
@constraint(omod23, y >= -2.1);
@NLobjective(omod23, Min, (1 - x)^2 + 100 * (y - x^2)^2);
optimize!(omod23);
This is Ipopt version 3.12.10, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).
Number of nonzeros in equality constraint Jacobian...: 0
Number of nonzeros in inequality constraint Jacobian.: 2
Number of nonzeros in Lagrangian Hessian.............: 3
Total number of variables............................: 2
variables with only lower bounds: 0
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 0
Total number of inequality constraints...............: 2
inequality constraints with only lower bounds: 2
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.0000000e+00 1.25e+00 1.50e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 3.3580384e+02 0.00e+00 9.92e+02 -1.0 1.35e+00 - 1.00e+00 1.00e+00h 1
2 1.2533591e-01 0.00e+00 2.71e-01 -1.0 1.83e+00 - 1.00e+00 1.00e+00f 1
3 9.6112888e-02 0.00e+00 1.04e+00 -1.7 1.21e-01 - 1.00e+00 1.00e+00f 1
4 8.3379322e-02 0.00e+00 2.37e-01 -1.7 5.21e-02 - 1.00e+00 1.00e+00f 1
5 8.1695749e-02 0.00e+00 4.83e-03 -1.7 7.05e-03 - 1.00e+00 1.00e+00f 1
6 6.4381341e-02 0.00e+00 5.23e-01 -3.8 8.31e-02 - 1.00e+00 1.00e+00f 1
7 6.2674893e-02 0.00e+00 6.35e-03 -3.8 6.91e-03 - 1.00e+00 1.00e+00f 1
8 6.2650298e-02 0.00e+00 1.40e-06 -3.8 1.13e-04 - 1.00e+00 1.00e+00f 1
9 6.2501984e-02 0.00e+00 4.39e-05 -5.7 7.41e-04 - 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 6.2501839e-02 0.00e+00 5.26e-11 -5.7 6.40e-07 - 1.00e+00 1.00e+00f 1
11 6.2499996e-02 0.00e+00 6.78e-09 -8.6 9.21e-06 - 1.00e+00 1.00e+00f 1
Number of Iterations....: 11
(scaled) (unscaled)
Objective...............: 6.2499996278424341e-02 6.2499996278424341e-02
Dual infeasibility......: 6.7847379936480934e-09 6.7847379936480934e-09
Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 2.5284241638252829e-09 2.5284241638252829e-09
Overall NLP error.......: 6.7847379936480934e-09 6.7847379936480934e-09
Number of objective function evaluations = 12
Number of objective gradient evaluations = 12
Number of equality constraint evaluations = 0
Number of inequality constraint evaluations = 12
Number of equality constraint Jacobian evaluations = 0
Number of inequality constraint Jacobian evaluations = 1
Number of Lagrangian Hessian evaluations = 11
Total CPU secs in IPOPT (w/o function evaluations) = 0.005
Total CPU secs in NLP function evaluations = 0.002
EXIT: Optimal Solution Found.
In [102]:
println("约束最小化目标函数值=", objective_value(omod22))
println("约束最小值点：x = ", value(x), " y = ", value(y))
约束最小化目标函数值=2.89460755048946
约束最小值点：x = 1.2499999925568486 y = 1.5624999813819718
方程求根¶
Roots包提供了一元方程求根算法，
主要是二分法之类不依赖于导数的方法。
比如，如下的一元函数
$$
f(x) = \sin(4(x - 0.25)) + x + x^{20} - 1,
\ x \in [0,1]
$$
有唯一一个实根。函数图形如下：
In [103]:
using Plots; Plots.gr()
rf1(x) = sin(4.0*(x - 0.25)) + x + x^20 - 1.0
Plots.plot(rf1, 0.0, 1.0, linewidth=2, legend=false)
Plots.hline!([0.0], color=:black)
Out[103]:
0.00
0.25
0.50
0.75
1.00
-1.5
-1.0
-0.5
0.0
0.5
1.0
Roots提供了若干种一元函数求根算法，比较基本的是二分法：
In [104]:
using Roots
Roots.find_zero(rf1, (0.0, 1.0), Bisection())
Out[104]:
0.40829350427936706
PolynomialRoots包提供了求多项式所有复根的功能。
比如，对多项式$f(x) = x^5 - x - 1$，
用系数升幂表示为[-1, -1, 0, 0, 0, 1]，
用roots()求根：
In [105]:
using PolynomialRoots
PolynomialRoots.roots([-1, -1, 0, 0, 0, 1])
Out[105]:
5-element Array{Complex{Float64},1}:
1.1673039782614187 - 1.1102230246251565e-16im
0.18123244446987524 + 1.083954101317711im
0.18123244446987535 - 1.0839541013177107im
-0.7648844336005848 + 0.3524715460317263im
-0.7648844336005848 - 0.35247154603172626im
并行计算¶
Julia语言是最近几年才发明的语言，
它吸收了现代最新的计算机科学技术成果，
比如，
语言中直接提供了对并行计算的支持。
因为并行计算与硬件密切相关，
所以这些功能在不同运行环境中有些能发挥作用而有些受限。
统计计算中有大量重复计算可以利用并行计算加速，
尤其是随机模拟问题，比较适合用并行计算。
有一些问题可以不依赖于语言支持而并行执行，
比如，
一个模型的随机模拟研究，
往往需要对模型设置许多组真实参数，
然后对每一组真实参数产生许多批次的样本进行模型估计。
这样，
可以运行多个非交互的Julia程序，
每个程序中使用不同的参数模拟，
最后汇总结果。
这样做的好处是不需要了解并行计算，
缺点是需要人为地开始多组程序，
当某些组任务完成时再人为地开始新的任务。
当每个任务耗时都很长，
比如若干个小时的时候，
这样的方法还是可行的。
最简单的并行运算是单一计算机的多个CPU核心(kernels或threads)同时独立地对一个数组的不同部分计算，
进一步还可以将若干台计算机组成集群(cluster)并行运行。
在统计中随机模拟经常需要对多个独立样本并行地执行一些计算，
这种情况下利用Julia的并行能力是有意义的，
即使单个计算机的4个核心同时运行能够将模拟速度提高一倍也是有意义的。
Julia的Distributed包是标准库的一部分，
提供了这样的功能，
实际上是在不同核心、不同CPU、不同计算机上执行单独的进程，
称这样的每个进程为一个工作进程或者节点。
在Linux或者MAC OS X中的ParallelAccelerator包可以将需要并行执行的代码编译为并行运行的本地二进制代码，
从而获得一两个数量级的加速。
Julia提供了coroutine功能，
也是一种并行机制，
离散系统模拟SimJulia就是利用这种机制。
Julia的Base.Threads包试验性地提供了多线程功能，
功能尚未成熟。
并行for循环和独立随机模拟¶
随机模拟程序中，
最简单的是独立样本的平均值法或者重要抽样法。
设随机变量$X$有分布密度$f(x)$，
为了计算积分$I=\int g(x) f(x) \,dx$，
注意到$Eg(X) = I$，
可以产生$X$的独立同分布样本$X_1, X_2, \dots, X_n$，
并估计$I=Eg(X)$为
$$
\hat I = \frac{1}{n} \sum_i g(X_i)
$$
作为示例，
计算积分$\int_0^1 e^x \,dx$,
精确值为$e - 1 \approx 1.718281828459045$。
取$X \sim \text{U}(0,1)$。
不用并行计算，程序为：
function f_para01(n)
s::Float64 = 0.0
for i = 1:n
s += exp(rand())
end
s /= n
return s
end
@time f_para01(1_000)
@time f_para01(1_000_000_000)
结果为
1.718281085841862
21.387429 seconds (7 allocations: 528 bytes)
Julia在标准库的Distributed包提供了@distributed for宏，
可以将循环自动划分范围分配到多个核心或者处理器。
同时，
还可以提供一个汇总操作，
将各个并行部分的结果汇总在一起。
修改为并行版本：
using Distributed
function f_para01b(n)
res = @distributed (+) for i = 1:n
exp(rand())
end
res /= n
return res
end
@time f_para01b(1_000)
@time f_para01b(1_000_000_000)
结果为
1.7182629178900357
7.630154 seconds (829 allocations: 36.438 KiB)
这是用了8个核心的结果，速度提高了大约3倍。
这个结果是在REPL命令行测试，
启动REPL命令行程序时加了选项-p 8表示有8个核心并行运行，
在Windows下可以修改程序的快捷方式在“目标”区域的程序后面加此选项，
如
D:\JuliaPro-1.4.0-1\Julia-1.4.0\bin\julia.exe -p 8
另一种办法是运行单核心的julia.exe命令行程序，
然后用如下命令增加可用的核心个数：
using Distributed
addprocs(7)
注意上述f_para01b()函数中没有利用统一的累加变量s，
这是因为不同的并行进程中无法同时修改同一个公用的变量。
另外，并行的for循环不能保证各步循环的先后次序。
@distributed for的办法是让每步循环都返回一个结果，
在每个并行过程中对其任务的所有循环汇总结果(本例中是求和)，
然后将所有并行过程的汇总结果再汇总(本例中是求和)。
并行运行中每个并行的进程都要有随机数，
所有进程需要使用独立的随机数序列。
经过测试，
上面的串行程序与并行程序得到的估计的标准差基本相同。
并行for循环中对公用的变量或数组写访问一般不可行，
读取则没有问题。
例如，
上面的模拟先将需要用的随机数一次性保存为数组，
然后在并行循环中读取随机数。
在10亿次模拟时内存不够出错，
所以将模拟分批进行，
如果在内层循环使用并行计算，
结果消耗内存较大而且速度比串行版本慢20倍；
下面的版本对各批分配到不同进程并行计算，
先不使用成批的随机数，
与不分批的并行计算速度相仿。
程序：
using Distributed
function f_para01c(n)
nb = 1000
n = Int(ceil(n/nb)*nb)
n1 = Int(n / nb)
res = @distributed (+) for j=1:nb
x = rand(n1)
s = 0.0
for i = 1:n1
s += exp(x[i])
end
s
end
res /= n
return res
end
@time f_para01c(1_000)
@time f_para01c(1_000_000_000)
结果为
1.7182779449240668
7.121955 seconds (841 allocations: 37.031 KiB)
速度与不分批类似，
但是可以在需要占用内存时利用有限的内存完成任务。
所有在有两重循环时应该将外层循环并行化而不是将内层循环并行化。
并行pmap()和互相不影响的计算¶
@distributed for比较适用于循环次数很多但每一步都很简单的情形。
对于每一步都可能比较复杂，
但是各步没有互相影响的计算，
可以用pmap()函数。
下面用不同的方法计算10个大的随机矩阵的奇异值：
using LinearAlgebra
rmatl = Matrix{Float64}[rand(100,100) for i=1:1000]
@time res = map(svdvals, rmatl);
其中调用LinearAlgebra包是Julia v1.0的要求，v0.6并不需要。结果为
3.273278 seconds (227.37 k allocations: 146.472 MiB, 4.24% gc time)
改用pmap()并行计算：
@time res = pmap(svdvals, rmatl);
结果：
1.444302 seconds (1.25 M allocations: 63.989 MiB, 0.38% gc time)
速度提高了1倍；
但是，如果是$1000 \times 1000$的矩阵执行10次重复计算，
速度仅能提高20%左右。
可见pmap()也是在重复次数比较多的时候提高效率明显。
所有节点需要的包和函数¶
并行计算时，
可以将主要的程序看作是“主控”节点，
用using调入的包、自定义的函数对主控节点是可见的，
但是如果其它工作节点也需要用到这些包或者函数，
就需要在所有节点中使其可用。
对某个包例如SharedArrays，
可以用@everywhere宏在所有节点调入，如
@everywhere using SharedArrays
对自定义函数或模块，
可以将源程序文件用@everywhere include("...")在所有节点载入。
比如，设文件f_para03.jl内容为：
f_para03(x) = sin(x) + cos(x)
运行如下代码：
@everywhere include("f_para03.jl")
res = @distributed (+) for i=1:100
f_para03(Float64(i))
end
结果为
-0.6594596218908123
因为用到并行计算的程序往往都是比较耗时的程序，
所以很可能会非交互地运行。
假设所有的并行工作节点都需要fun1.jl和fun2.jl中的程序，
而drv.jl是主控程序，
可以用如下的操作系统命令非交互地并行运行(主控程序中用了@distributed、pmap()等并行控制)：
julia -p 8 -L fun1.jl -L fun2.jl drv.jl
其中-p 8选项表示有8个核心可以并行使用。
-L选项给出每个工作节点需要预先调用的初始化程序。
共享数组¶
并行计算时，
读取共同的数组没有问题；
写入共同的数组则有问题，
因为不同的进程会使用数组的副本。
所以，
下面的程序不会使得公用数组a的第$i$个元素等于$i$:
using Distributed
a = zeros(10)
@distributed for i=1:10
a[i] = i
end
a
结果中a的元素仍然都是零。
如果并行计算的各个部分不会要求写入共同的存储位置，
可以用SharedArray作为公用的可写数组。如
using Distributed
@everywhere using SharedArrays
b = SharedArray{Int,1}(10)
@sync @distributed for i=1:10
b[i] = i
end
b
其中@everywhere语句将SharedArrays包提供给每个工作进程，
@sync关键字要求for循环结束时所有并行执行的循环已经完成，
没有此关键字在查看b结果是可能有些循环步还没有完成。
比例的置信区间模拟研究并行计算示例¶
问题¶
统计计算的一部分工作是用随机模拟比较不同的统计方法的性能。
设总体某属性的真实比例为$p$，
n个样本中具有该属性的个数为x，
Julia的HypothesisTests包提供了$p$的置信区间的若干种计算方法。
如下程序
Using HypothesisTests
testr = BinomialTest(x, n, 0.5)
产生一个二项分布检验结果，
然后用
confint(testr, alpha, method=:wald)
可以用正态近似方法求置信度为1-alpha的Wald置信区间，
这种方法在真实比例$p$接近0或者1时效果很差。
如下程序
confint(testr, alpha, method=:wilson)
计算基于正态近似的Wilson置信区间，
效果优于Wald方法。
希望对不同的样本量$n$，不同的真实比例$p$，
对每一种$(n,p)$组合，
产生$M$批随机样本，
分别用两种方法计算置信度为0.90、0.95和0.99的置信区间，
比较两种方法得到的区间的覆盖率与置信度。
覆盖率是指置信区间包含真实比例的百分比。
设每个组合的重复样本数$M=10^8$，
样本量$n$取100, 500, 1000,
真实比例$p$取0.5, 0.05, 0.01，
共有$3 \times 3 = 9$种参数组合
(实际问题一般有几十种或几百种参数组合)。
一组参数的模拟程序¶
对每一组$(n, p)$值，
可以用如下的程序进行一次模拟：
## 比例的置信区间模拟比较。
## 一个样本量和比例组合的多批样本的模拟。
using Random, Distributions, HypothesisTests
function f_paraci(
n::Int, # 样本量
p::Float64, # 真实比例
M::Int) # 重复模拟次数
# 记录覆盖次数，第一行为Wald方法，第二行为Wilson方法
# 三列分别对应于三个水平alpha值
ncov = [0 0 0;
0 0 0]
alphas = [0.10, 0.05, 0.01]
for i = 1:M
x::Int = rand(Binomial(n, p))
if x==0 || x==n
continue
end
testr = BinomialTest(x, n, 0.5)
for j in eachindex(alphas)
low,hi = confint(testr, alphas[j], method=:wald)
if low < p < hi
ncov[1,j] += 1
end
low,hi = confint(testr, alphas[j], method=:wilson)
if low < p < hi
ncov[2,j] += 1
end
end # for j
end # for i
cover = ncov ./ M
println("n=", n, " p=", p)
println(cover)
cover
end
为了对某一个参数组合计算，
只要在Julia命令行REPL中运行命令如
@time f_paraci(1000, 0.01, 10_000_000)
结果如
n=1000 p=0.01
[0.9074386 0.9269796 0.9706971; 0.9234359 0.9634912 0.9926302]
9.342665 seconds (62 allocations: 4.219 KiB)
2×3 Array{Float64,2}:
0.907439 0.92698 0.970697
0.923436 0.963491 0.99263
取$M=10^7$时每组参数模拟需要10秒左右。
人工并行计算¶
如果有9组参数需要模拟，
假设上面的对一组参数模拟的函数保存在源文件simu-ci.jl中，
可以制作如下的9个源文件，
第一个命名为simu01.jl，内容为
include("simu-ci.jl")
print(f_paraci(100, 0.5, 10_000_000))
其它的8个源文件只是$(n,p)$组合值不同。
在MS Windows操作系统中(Linux类似)，
假设有4个CPU核心，
可以在源程序所在的文件夹打开4个操作系统命令窗口
(MS Windows是CMD或者PowerShell, Linux是Bash Shell)，
分别运行命令如
D:\JuliaPro-1.4.0-1\Julia-1.4.0\bin\julia.exe simu01.jl
其中的julia.exe路径应该替换成自己的Julia可执行程序路径。
等待每个窗口的程序执行结束，
就依次调用下一个模拟源程序。
这种做法好处是不需要利用语言本身的并行计算功能，
在参数组合不太多，
而每个参数组合运行时间很长的时候是很实用的做法。
但是组合太多的时候人为调度任务就太麻烦了。
利用Julia分布式计算¶
设一次模拟的函数f_paraci()定义在simu-ci.jl源文件中。
定义如下的主控文件simu-main.jl：
using Distributed
function simu_all()
M = 10_000_000
pars = [(n, p) for n in [100, 500, 1000], p in [0.5, 0.05, 0.01]]
res = pmap(parone -> f_paraci(parone[1], parone[2], M), pars)
end
@time res = simu_all()
打开一个单独的julia命令行，
用如下命令增加可核心(工作节点)个数，
并在工作节点载入所需的包和模拟一次用的函数：
using Distributed
addprocs(7)
@everywhere using Random
@everywhere using Distributions
@everywhere using StatsBase
@everywhere using HypothesisTests
@everywhere include("simu-ci.jl")
然后
include("simu-main.jl")
可以使用8个节点并行计算。
另一种方法是非交互运行，
在操作系统命令行中运行：
D:\JuliaPro-1.4.0-1\Julia-1.4.0\bin\julia.exe -p 8 -L simu-ci.jl simu-main.jl
但是在我的JuliaPro1.4.0 Window版系统中，
非交互式运行总提示在工作节点处找不到某些要载入的包如Distributions。
经过试验，
取$M=10^7$时在一台Intel 酷睿i7 8550U 8核心笔记本电脑上运行了41秒。
将程序改为如下非并行版本，运行了79秒，
并行版本速度提高了1倍左右。
当参数组合更多，
每次建模耗时更多，
重复数更多时，
并行计算是很有意义的。
include("simu-ci.jl")
function simu_all_single()
M = 10_000_000
pars = [(n, p) for n in [100, 500, 1000], p in [0.5, 0.05, 0.01]]
res = map(parone -> f_paraci(parone[1], parone[2], M), pars)
end
@time res = simu_all_single()
多只股票投资策略实证分析¶
证券投资策略研究中往往需要对多支证券的历史数据进行投资策略的实证分析，
从而找到在历史数据上较优的策略。
对一支股票，
策略研究需要考虑的参数也有计算区间的长度选择，
策略实施的区间长度选择等变量，
所以组合起来有很大的计算量，
有必要采用并行计算加速，
Julia语言对于显式循环的高效率也使得Julia语言比较适用于金融投资策略研究。
这个示例使用Julia v1.4.0运行。
首先模拟生成100支股票1000天的股价数据：
In [106]:
using Distributions
function price_sim(nstocks=100, ndays=1000)
pdata = Matrix{Float64}(undef, (ndays, nstocks))
for j=1:nstocks
# j: 每只股票
# 每只股票有自己的波动率数值
sig = rand(Uniform(0.01, 0.02))
pdata[:,j] .= round.(
rand(Uniform(10.0, 100.0)) .*
exp.(cumsum(randn(ndays) .* sig)),
digits = 2)
end
return pdata
end
pdata = price_sim(100, 1000);
using Plots
Plots.gr()
Plots.plot(pdata[:,1:5])
Out[106]:
0
250
500
750
1000
50
100
150
200
y1
y2
y3
y4
y5
上面的结果显示了模拟的100支股票中前5支的价格序列。
假设对每支股票采用如下的投资策略回测研究：
依次截取$n_1$天的价格统计其20%分位数$\underline{x}$和80%分位数$\bar x$(也可以是更一般的分位数)，
然后在后续的$n_2$天中(经常是一个月)当价格首次低于$\underline{x}$时就开仓买入一个单位，
然后在这$n_2$天中当价格首次高于$\bar x$时就卖出，到$n_2$天的最后一天时不管价位如何都卖出。
如果在$n_2$天中卖出了，卖出的第二天开始后续的$n_1$天的统计，
如果没有开仓，则在$n_2$天结束后开始第二轮的$n_1$天的统计。
如果需要试验不同的$n_1, n_2$的取法组合，
计算量会比较大。
对一支股票的策略回测函数为:
function testone(x::Vector{Float64}, n1=30, n2=30)
ndays = length(x)
t = 1
earn = 0.0
while t < ndays - (n1 + n2)
# 首先进行n1天的统计
lo, hi = quantile(x[t:(t+n1-1)], [0.2, 0.8])
# 然后在后续的n2天寻找开仓时机
t += n1
t1 = t
p_buy = NaN
p_sell = NaN
for t1 in t:(t+n2-2)
if x[t1] < lo
p_buy = x[t1]
break
end
end
if !isnan(p_buy) # 有买入，查找从t1+1到t+n2-1天的卖出时机
for t2=(t1+1):(t+n2-1)
if x[t2] > hi
p_sell = x[t2]
earn += p_sell - p_buy
t = t2 + 1 # 开始下一轮
break
end
end
# 如果没有卖出就在n2天的最后一天平仓
if isnan(p_sell)
p_sell = x[t+n2-1]
earn += p_sell - p_buy
t += n2
end
end
end # while t
return earn
end
function testn(x, n1vec, n2vec)
[testone(x, n1, n2) for n1 in n1vec, n2 in n2vec]
end
设以上函数存入了sim-stock-fun.jl中。
比较串行和并行程序：
using Distributed
addprocs(7)
@everywhere include("sim-stock-fun.jl")
pdata = price_sim(1000, 3000);
串行版本1：
@time for j = 1:size(pdata,2)
x = testn(pdata[:,j], 20:40, 30:50)
end
## 43.397073 seconds (172.90 M allocations: 23.689 GiB, 3.79% gc time)
串行版本2：
@time map(j -> testn(pdata[:,j], 20:40, 30:50),
1:size(pdata,2));
## 49.662468 seconds (173.10 M allocations: 23.699 GiB, 3.44% gc time)
并行版本：
using Distributed
@time pmap(j -> testn(pdata[:,j], 20:40, 30:50),
1:size(pdata,2));
## 19.880122 seconds (1.02 M allocations: 53.250 MiB, 0.03% gc time)
结果发现小规模计算时并行版本更慢，
串行版本的map()函数分配内存次数和总量都很大，
pmap()则相对来说很少。
较大计算量时串行计算CPU占用率仅30%左右，并行计算时CPU占用有80%左右，
3000支股票1000天的数据，建模窗宽与执行窗口共441种组合，
都能在几十秒内完成，
并行算法快一倍多。
更复杂的建模可能会速度慢得多，
这时并行的速度提升就很有意义；
即使仅使用串行算法，
使用Julia会比同样方法的R、Python等程序快10倍以上。
895
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
