-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path执行效率对比 matlab vs python.Rmd
152 lines (110 loc) · 4.3 KB
/
执行效率对比 matlab vs python.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
Matlab 在数值计算方面有很大的优势，用 Matlab 和 python 分别解一个带参数的
积分方程结果非常出乎意料。
F(p; b, k) 是 f(x, p) 从 0 到 1 的积分，从而是一个 p 的函数。b, k 为参数。
这里 f 是一个和威布尔分布相关的函数，而 F 则是两个威布尔随机变量和的
cdf。其中 b 是两个变量共同的形状参数，k 是它们尺度参数的比值。由于位置参数
不影响 F 的取值，因此无需纳入考虑。
问题的背景及 F 的推导不作具体讨论。现在只需知道：
$$
F(p; b, k) = \int_0^1 f(x, p; b, k) dx
$$
而我们关心的是，函数 F 的不动点，即：
$F(p) = p$
的根 rp。
因为 0 和 1 是函数 f 的瑕点，因此我们只关心 rp 在 (0,1) 之间的情况。b, k 是 F
的参数，因此最终 rp 是 b, k 的函数，记为 rp(b, k) 下面编程对 rp 的值进行数值
求解。
python:
```{python}
import pandas as pd
pd.set_option("display.precision", 5)
import numpy as np
from scipy.stats import weibull_min
from scipy.integrate import quad,quadrature
from scipy.optimize import brentq
rv = weibull_min
def f(x, p, b=1, k=1):
return rv.cdf((1+k)*rv.ppf(p, b) - k*rv.ppf(x, b), b)
def F(p, b=1, k=1, eps=1e-8):
return quadrature(f, 0, 1, args=(p, b, k), tol=eps, rtol=eps,maxiter=200)[0]
# return quad(f, 0, 1, args=(p, b, k), epsabs=eps, epsrel=eps)[0]
def root_p(b=1, k=1, eps=1e-8, xtol=1e-6):
return brentq(lambda p: F(p, b, k, eps=eps)-p, 0.01, 0.99, xtol=xtol, rtol=xtol)
v = np.vectorize(root_p)
eps = 1e-8
xtol = 1e-6
if __name__ == '__main__':
from time import clock
start = clock()
k = np.linspace(0.01, 0.99, 6)
b = [0.5, 1, 2.5, 4, 10]
R = pd.DataFrame(index=k, columns=b)
for ib in b:
R.loc[:, ib] = v(ib, k, eps=eps, xtol=xtol)
print(R)
print("elapsed: {:.2f} s".format(clock() - start))
'''
quad()
Result
0.5 1.0 2.5 4.0 10.0
0.010 0.76246 0.63397 0.5239 0.49067 0.45448
0.206 0.84415 0.6718 0.52915 0.4875 0.44145
0.402 0.86338 0.69725 0.53342 0.48473 0.43143
0.598 0.87042 0.70916 0.5364 0.48276 0.4254
0.794 0.87312 0.71406 0.538 0.48171 0.42256
0.990 0.8738 0.71533 0.53847 0.48141 0.42179
elapsed: 76.15 s
quadrature(maxiter=200)
0.5 1.0 2.5 4.0 10.0
0.010 0.76232 0.63393 0.52388 0.49065 0.4545
0.206 0.84378 0.67177 0.52915 0.48749 0.44145
0.402 0.86381 0.69726 0.53342 0.48473 0.43143
0.598 0.87026 0.70917 0.53641 0.48276 0.42539
0.794 0.87271 0.71404 0.538 0.48171 0.42256
0.990 0.87379 0.71534 0.53847 0.48141 0.42179
elapsed: 19.03 s
'''
```
matlab：
```{r,engine='matlab',engine.opts='--no-display -r'}
start_tic = tic;
rv = prob.WeibullDistribution();
k = linspace(0.01, 0.99, 6);
b = [0.5 1 2.5 4 10];
R = zeros(length(k), length(b));
eps = 1e-8;
opt = optimset('TolX', 1e-6);
% matlab can't defines functions in script, use func handles instead
% note, the parameter b is passed as a property of rv.
f = @(x, p, k, rv) rv.cdf((1+k).*rv.icdf(p) - k.*rv.icdf(x));
F = @(p, k, rv)integral(@(x)f(x, p, k, rv), 0, 1, 'abstol', eps, 'reltol', eps);
root_p = @(k, rv)fzero(@(p)F(p, k, rv)-p, [0.01, 0.99], opt);
for j = 1:length(b)
rv.B = b(j);
for i = 1:length(k)
R(i,j) = root_p(k(i), rv);
end
end
R
fprintf('%s elapsed: %f s\n', mfilename, toc(start_tic));
exit
```
```{r}
system("matlab -nodisplay -r \"run('matpy/test.m'); exit\"",intern = T)
# "R ="
# [13] ""
# [14] " 0.7625 0.6340 0.5239 0.4907 0.4545"
# [15] " 0.8441 0.6718 0.5291 0.4875 0.4415"
# [16] " 0.8634 0.6973 0.5334 0.4847 0.4314"
# [17] " 0.8704 0.7092 0.5364 0.4828 0.4254"
# [18] " 0.8731 0.7141 0.5380 0.4817 0.4226"
# [19] " 0.8738 0.7153 0.5385 0.4814 0.4218"
# [20] ""
# [20] ""
# [21] "test elapsed: 2.016592 s"
```
执行耗时如下：
> matlab: 0.941343 s
>
> python: 111.63 s
matlab 的计算效率大约是 python 的 120 倍！！