teratail
回答率
85.47%
質問するログイン新規登録
トップに関する質問optimize.basinhoppingの[setting an array element with a sequence]エラー

編集履歴

質問編集履歴

1

大幅に質問内容を変更しました。

2019/12/14 16:23

投稿

Jagali
Jagali

スコア4

test CHANGED
@@ -1 +1 @@
1
- least_squaresの使い方
1
+ optimize.basinhoppingの[setting an array element with a sequence]エラー
test CHANGED
@@ -4,17 +4,23 @@
4
4
 
5
5
  逆問題の求解をしています。
6
6
 
7
- 連立微分方程式を数値積分し、実測で得られているデータと最もFitするように、
7
+ 未知のパラメーターを含む連立微分方程式が、実測で得られているデータと最もFitするように、
8
-
9
- ｚを最小二乗法で最適化したいです。
8
+
10
-
11
-
12
-
13
- 最後の行の最適化のleastsqがどうしてもうまく動作させられません。
9
+ scipy.optimize.basinhoppingで最適パラメーターを見つけたいです。
10
+
14
-
11
+ リファレンス：
12
+
13
+ https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.basinhopping.html#scipy.optimize.basinhopping
14
+
15
+
16
+
17
+ プログラムを実行すると、for ループが2回終了したところで、
18
+
15
- 実行するとエラーは出ないのですが、最適化結果が初期値のまま返されてしまいます。
19
+ 「setting an array element with a sequence.」と表示されてしまいます。
20
+
21
+
22
+
16
-
23
+ 途中からエラーを起こす理由がよくわかりません。
17
-
18
24
 
19
25
  修正点を教えていただけないでしょうか？
20
26
 
@@ -22,382 +28,284 @@
22
28
 
23
29
 
24
30
 
25
- ### 発生している問題・エラーメッセージ
31
+ ### コード
26
-
32
+
27
- ```Python
33
+ ```Python3
28
34
 
29
35
  import numpy as np
30
36
 
31
- import matplotlib as mpl
32
-
33
- import matplotlib.pyplot as plt
34
-
35
37
  from scipy.integrate import odeint
36
38
 
37
39
  import pandas as pd
38
40
 
39
41
  from scipy import optimize
40
42
 
43
+ from scipy.optimize import basinhopping
44
+
45
+
46
+
47
+ n = 3 #データ数
48
+
49
+ date = "1812" #シート名
50
+
51
+
52
+
53
+ D = pd.read_excel('前データ.xlsx', sheet_name=date,dtype='float64')
54
+
41
- import scipy.stats
55
+ DD = np.array(D)
56
+
42
-
57
+ c = DD[5:12,0]
43
-
44
-
58
+
59
+
60
+
45
- def func(y, t, z):
61
+ def f(y, t, z):
46
62
 
47
63
  #単位式
48
64
 
49
- R1 = z[3] * y[0] * y[1]
65
+ R1 = c[3] * y[0] * y[1]
50
-
66
+
51
- R2 = z[4] * z[0] * y[1]
67
+ R2 = c[4] * c[0] * y[1]
52
-
68
+
53
- R3 = z[5] * z[1] * y[1]
69
+ R3 = c[5] * c[1] * y[1]
54
-
70
+
55
- R4 = z[6] * y[4] * y[1]
71
+ R4 = c[6] * y[4] * y[1]
56
-
72
+
57
- R5 = z[7] * y[0] * y[2]
73
+ R5 = z[0] * y[0] * y[2]
58
-
74
+
59
- R6 = z[8] * y[0] * y[3]
75
+ R6 = z[1] * y[0] * y[3]
60
-
76
+
61
- R7 = z[9] * y[1] * y[2]
77
+ R7 = z[2] * y[1] * y[2]
62
-
78
+
63
- R8 = z[10] * y[1] * y[3]
79
+ R8 = z[3] * y[1] * y[3]
64
80
 
65
81
  #連立微分方程式
66
82
 
67
- y0 = -R1 -z[2]*R5 -z[2]*R6
68
-
69
- y1 = -R1 -R2 -R3 -R4 +z[11]*z[2]*R5 +z[12]*z[2]*R6 -z[2]*R7 -z[2]*R8
70
-
71
- y2 = -z[13]*R5 -z[14]*R7
72
-
73
- y3 = -z[15]*R6 -z[16]*R8
74
-
75
- y4 = -R4
76
-
77
- return [y0, y1, y2, y3, y4]
78
-
79
-
80
-
81
- #微分方程式の初期値
82
-
83
- Y = pd.read_excel('TEST.xlsx', skiprows = [6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22])
84
-
85
- y = Y["初期値"]
86
-
87
- Z = pd.read_excel('TEST.xlsx', skiprows = [1,2,3,4,5])
88
-
89
- z = Z["初期値"]
90
-
91
-
92
-
93
- #実測データ
94
-
95
- D = pd.read_excel('DATA.xlsx', skiprows = [10,11,12,13,14,15,16,17,18,19,20,21,22])
96
-
97
- data_time = D["Time"]
98
-
99
- data_O = np.array(D["O"])
100
-
101
- data_P = np.array(D["P"])
102
-
103
- data_OP = scipy.stats.zscore(np.concatenate([data_O, data_P])) #結合して標準化
104
-
105
-
106
-
107
- #積分の実行と誤差を計算
108
-
109
- t = np.arange(0, 900, 0.01)
110
-
111
- def ode_era(z):
112
-
113
- ode = odeint(func, y, t, args=(z,))
114
-
115
- ode_O = ode[data_time,0]
116
-
117
- ode_P = ode[data_time,4]
118
-
119
- ode_OP = scipy.stats.zscore(np.concatenate([ode_O, ode_P])) #結合して標準化
120
-
121
- return ode_OP - data_OP
122
-
123
-
124
-
125
- #最適化の実行
126
-
127
- ANS = optimize.least_squares(ode_era,z)
83
+ e0 = -R1 -c[2]*R5 -c[2]*R6
84
+
85
+ e1 = -R1 -R2 -R3 -R4 +z[4]*c[2]*R5 +z[5]*c[2]*R6 -c[2]*R7 -c[2]*R8
86
+
87
+ e2 = -z[6]*R5 -z[7]*R7
88
+
89
+ e3 = -z[8]*R6 -z[9]*R8
90
+
91
+ e4 = -R4
92
+
93
+ return [e0, e1, e2, e3, e4]
94
+
95
+
96
+
97
+ def ODE(x, teta, i): 　　　　　　　　#数値積分の実行
98
+
99
+ f2 = lambda y,t: f(y, t, teta)
100
+
101
+ y0 = DD[0:5, i]
102
+
103
+ r = odeint(f2,y0,x)
104
+
105
+ r_O = scipy.stats.zscore(r[:,0])
106
+
107
+ r_P = scipy.stats.zscore(r[:,4])
108
+
109
+ return np.concatenate([r_O, r_P])
110
+
111
+
112
+
113
+
114
+
115
+ def ODE_resid(p):　　　　　　　　　　#実測値と計算値の比較
116
+
117
+ for i in range(n):
118
+
119
+ x_data = np.array(DD[22:31, i])
120
+
121
+ data_Oa = scipy.stats.zscore(DD[33:42,i])
122
+
123
+ data_Pa = scipy.stats.zscore(DD[44:53,i])
124
+
125
+ y_data = np.concatenate([data_Oa, data_Pa])
126
+
127
+ if i == 0:
128
+
129
+ print(i) #何回ぐらい計算されているのか把握するために書いています。
130
+
131
+ Res0 = y_data - ODE(x_data,p,i))
132
+
133
+ if i == 1:
134
+
135
+ print(i)
136
+
137
+ Res1 = y_data - ODE(x_data,p,i))
138
+
139
+ if i == 2:
140
+
141
+ Res2 = y_data - ODE(x_data,p,i))
142
+
143
+ print(i)
144
+
145
+ Q = np.concatenate([Res0,Res1,Res2])
146
+
147
+ return Q
148
+
149
+
150
+
151
+ g = DD[12:22,0] #推定するパラメーター
152
+
153
+ j = basinhopping(ODE_resid,g)
128
154
 
129
155
  ```
130
156
 
131
- ###実行結果
132
-
133
- active_mask: array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
134
-
135
- cost: 2.1925926509284848
136
-
137
- fun: array([-1.50504475, -0.21617593, 0.04835059, 0.18241003, 0.3565212 ,
138
-
139
- 0.39604178, 0.53540293, 0.71837111, 0.86037334, -0.16863658,
140
-
141
- -0.15435401, -0.15355698, -0.1520566 , -0.15161362, -0.1505089 ,
142
-
143
- -0.14947951, -0.14875215, -0.14729194])
144
-
145
- grad: array([-9.73025931e+01, -7.58921428e+02, -4.54999004e-01, -8.92974021e-11,
146
-
147
- -9.75097653e-10, -9.75373824e-13, -2.14089495e-13, -1.37899984e+01,
148
-
149
- -2.51843576e+01, 3.10083970e-06, -2.06437191e-06, 1.42880394e-01,
150
-
151
- 1.16386908e-02, 3.38799083e-07, 1.21195919e-07, 3.93826473e-09,
152
-
153
- 1.55404385e-08])
154
-
155
- jac: array([[ 1.95728561e+01, 1.52660147e+02, 1.24124534e-01,
156
-
157
- 2.23517418e-11, 1.96146071e-10, 1.96198622e-13,
158
-
159
- 3.66836786e-14, 4.09854011e+00, 6.34926283e+00,
160
-
161
- -6.25430600e-07, 4.19032802e-07, -2.92027146e-02,
162
-
163
- -2.29141116e-03, -6.75410848e-08, -2.47595839e-08,
164
-
165
- -7.09038485e-10, -2.83547004e-09],
166
-
167
- [ 1.83482961e+01, 1.43109124e+02, 8.77128472e-02,
168
-
169
- 1.71197785e-11, 1.83873177e-10, 1.83819196e-13,
170
-
171
- 3.86714935e-14, 2.67085925e+00, 4.78946735e+00,
172
-
173
- -5.86305459e-07, 3.91952268e-07, -2.72613913e-02,
174
-
175
- -2.15850770e-03, -6.34588676e-08, -2.31257390e-08,
176
-
177
- -6.81272890e-10, -2.72354359e-09],
178
-
179
- [ 1.57356596e+01, 1.22731712e+02, 5.59137823e-02,
180
-
181
- 1.20864974e-11, 1.57690644e-10, 1.57646644e-13,
182
-
183
- 3.58372927e-14, 1.50582054e+00, 3.34205659e+00,
184
-
185
- -5.01829900e-07, 3.34177260e-07, -2.31257975e-02,
186
-
187
- -1.87598169e-03, -5.46998450e-08, -1.96204631e-08,
188
-
189
- -6.27096119e-10, -2.48563097e-09],
190
-
191
- [ 1.21088136e+01, 9.44438857e+01, 2.78799429e-02,
192
-
193
- 7.25189845e-12, 1.21344030e-10, 1.21310735e-13,
194
-
195
- 2.95013189e-14, 5.46571702e-01, 1.98685271e+00,
196
-
197
- -3.84589334e-07, 2.54502104e-07, -1.74638182e-02,
198
-
199
- -1.47801638e-03, -4.24756260e-08, -1.48147614e-08,
200
-
201
- -5.43122125e-10, -2.11557529e-09],
202
-
203
- [ 7.74420722e+00, 6.04018116e+01, 2.95658473e-03,
204
-
205
- 2.59942479e-12, 7.76058435e-11, 7.76006625e-14,
206
-
207
- 2.06261873e-14, -2.49508262e-01, 7.08309278e-01,
208
-
209
- -2.43832929e-07, 1.59352461e-07, -1.07533783e-02,
210
-
211
- -9.90629196e-04, -2.76807394e-08, -9.11719171e-09,
212
-
213
- -4.23256020e-10, -1.61549013e-09],
214
-
215
- [ 2.84809741e+00, 2.22141979e+01, -1.93690994e-02,
216
-
217
- -1.80469619e-12, 2.85416842e-11, 2.86942873e-14,
218
-
219
- 9.92715359e-15, -9.14803527e-01, -5.05721875e-01,
220
-
221
- -8.64181073e-08, 5.34606476e-08, -3.33912671e-03,
222
-
223
- -4.33608890e-04, -1.09799695e-08, -2.82214791e-09,
224
-
225
- -2.70545248e-10, -9.90889080e-10],
226
-
227
- [-7.95559648e+00, -6.20502277e+01, -5.77755334e-02,
228
-
229
- -1.00831191e-11, -7.97238946e-11, -7.95492759e-14,
230
-
231
- -1.49369240e-14, -1.94722860e+00, -2.77682558e+00,
232
-
233
- 2.59154170e-07, -1.77395636e-07, 1.26663744e-02,
234
-
235
- 8.31454992e-04, 2.62584159e-08, 1.07628339e-08,
236
-
237
- 1.23252153e-10, 6.00708700e-10],
238
-
239
- [-2.53094646e+01, -1.97403638e+02, -1.03826956e-01,
240
-
241
- -2.13252174e-11, -2.53632814e-10, -2.53587197e-13,
242
-
243
- -5.62906265e-14, -2.98512446e+00, -5.88672283e+00,
244
-
245
- 8.09406471e-07, -5.41013136e-07, 3.75363827e-02,
246
-
247
- 2.96240300e-03, 8.72297927e-08, 3.18533940e-08,
248
-
249
- 9.10846955e-10, 3.71502998e-09],
250
-
251
- [-5.41890573e+01, -4.22653836e+02, -1.60715689e-01,
252
-
253
- -3.72032324e-11, -5.43043762e-10, -5.43452990e-13,
254
-
255
- -1.25680864e-13, -3.94150492e+00, -1.05371802e+01,
256
-
257
- 1.71269955e-06, -1.12816190e-06, 7.70774856e-02,
258
-
259
- 6.77932799e-03, 1.92192504e-07, 6.53288732e-08,
260
-
261
- 2.69393992e-09, 1.03251690e-08],
262
-
263
- [ 4.10598353e-01, 3.20267879e+00, 5.18598632e-03,
264
-
265
- 8.27842289e-13, 4.11465764e-12, 4.14557946e-15,
266
-
267
- 2.89461017e-13, 1.80118084e-01, 3.10704015e-01,
268
-
269
- -6.06131548e-08, 7.99999150e-09, -6.00621104e-04,
270
-
271
- -4.88609076e-05, -4.81845800e-09, -1.70236933e-09,
272
-
273
- -4.06325779e-11, -1.56448916e-10],
274
-
275
- [ 6.22727722e-01, 4.85718916e+00, 5.07100766e-03,
276
-
277
- 8.69234403e-13, 6.24045730e-12, 6.26612932e-15,
278
-
279
- 2.17191875e-13, 1.66466482e-01, 3.03462721e-01,
280
-
281
- -5.55091318e-08, 1.28484495e-08, -9.43087041e-04,
282
-
283
- -7.12499022e-05, -4.79148855e-09, -1.69171703e-09,
284
-
285
- -4.09711827e-11, -1.58148948e-10],
286
-
287
- [ 8.16178910e-01, 6.36601689e+00, 4.97088482e-03,
288
-
289
- 9.02348095e-13, 8.17909837e-12, 8.21474271e-15,
290
-
291
- 1.50910020e-13, 1.54745534e-01, 2.96802878e-01,
292
-
293
- -5.07464678e-08, 1.72061830e-08, -1.24298036e-03,
294
-
295
- -9.25511122e-05, -4.73654145e-09, -1.66991468e-09,
296
-
297
- -4.19869972e-11, -1.63432832e-10],
298
-
299
- [ 9.93139587e-01, 7.74622534e+00, 4.88314937e-03,
300
-
301
- 9.43740209e-13, 9.95248556e-12, 9.97231557e-15,
302
-
303
- 8.98763537e-14, 1.44635670e-01, 2.90637828e-01,
304
-
305
- -4.62837633e-08, 2.11348463e-08, -1.50686502e-03,
306
-
307
- -1.12913549e-04, -4.65903578e-09, -1.63950090e-09,
308
-
309
- -4.40186261e-11, -1.71657312e-10],
310
-
311
- [ 1.15549623e+00, 9.01252853e+00, 4.80581234e-03,
312
-
313
- 9.76853900e-13, 1.15795434e-11, 1.15961601e-14,
314
-
315
- 3.34680080e-14, 1.35878466e-01, 2.84896560e-01,
316
-
317
- -4.20850729e-08, 2.46871321e-08, -1.74012780e-03,
318
-
319
- -1.32463872e-04, -4.56250027e-09, -1.60244809e-09,
320
-
321
- -4.63888598e-11, -1.82225081e-10],
322
-
323
- [ 1.30487667e+00, 1.01776232e+01, 4.73729292e-03,
324
-
325
- 9.85132323e-13, 1.30763650e-11, 1.30862762e-14,
326
-
327
- -1.88454986e-14, 1.28263615e-01, 2.79520512e-01,
328
-
329
- -3.81228795e-08, 2.79103566e-08, -1.94719434e-03,
330
-
331
- -1.51269138e-04, -4.45071567e-09, -1.56034292e-09,
332
-
333
- -4.94363031e-11, -1.94630722e-10],
334
-
335
- [ 1.57016771e+00, 1.22467628e+01, 4.62174413e-03,
336
-
337
- 1.05135971e-12, 1.57351792e-11, 1.57608436e-14,
338
-
339
- -1.12906098e-13, 1.15801804e-01, 2.69677393e-01,
340
-
341
- -3.08125423e-08, 3.35201723e-08, -2.29682773e-03,
342
-
343
- -1.86972320e-04, -4.19059909e-09, -1.46567317e-09,
344
-
345
- -5.65470043e-11, -2.22979908e-10],
346
-
347
- [ 1.90061721e+00, 1.48241058e+01, 4.48848240e-03,
348
-
349
- 1.10930867e-12, 1.90466642e-11, 1.90658447e-14,
350
-
351
- -2.32651830e-13, 1.02232523e-01, 2.56664746e-01,
352
-
353
- -2.11213688e-08, 4.02820613e-08, -2.69923359e-03,
354
-
355
- -2.36697495e-04, -3.73854571e-09, -1.30980252e-09,
356
-
357
- -6.90753825e-11, -2.71040280e-10],
358
-
359
- [ 2.32238586e+00, 1.81136940e+01, 4.33520776e-03,
360
-
361
- 1.18381447e-12, 2.32733786e-11, 2.33069444e-14,
362
-
363
- -3.90923023e-13, 8.82357359e-02, 2.38134868e-01,
364
-
365
- -7.56195624e-09, 4.84992764e-08, -3.15720588e-03,
366
-
367
- -3.12075019e-04, -2.89858592e-09, -1.04356368e-09,
368
-
369
- -9.27777196e-11, -3.54295915e-10]])
370
-
371
- message: '`xtol` termination condition is satisfied.'
372
-
373
- nfev: 16
374
-
375
- njev: 1
376
-
377
- optimality: 758.9214282478794
378
-
379
- status: 3
380
-
381
- success: True
382
-
383
- x: array([5.01059643e-04, 5.01059643e-07, 1.15000000e+00, 9.00000000e+05,
384
-
385
- 5.00000000e+07, 3.90000000e+08, 5.00000000e+09, 2.63922107e-02,
386
-
387
- 1.66082493e-03, 3.00547243e+04, 2.07852793e+04, 9.98215810e-01,
388
-
389
- 4.61445490e-01, 5.91196391e+04, 1.19743194e+06, 2.20037645e+04,
390
-
391
- 1.62156618e+05])
157
+
158
+
159
+ ###エラー文
160
+
161
+ ```ここに言語を入力
162
+
163
+ ValueError Traceback (most recent call last)
164
+
165
+ <ipython-input-32-a4a40533ee77> in <module>
166
+
167
+ 68 (g[6],g[6]*100),(g[7],g[7]*100),(g[8],g[8]*100),(g[9],g[9]*100)]
168
+
169
+ 69
170
+
171
+ ---> 70 j = basinhopping(ODE_resid,g)
172
+
173
+ 71 # jj = j.x
174
+
175
+ 72 # df = pd.DataFrame(jj)
176
+
177
+
178
+
179
+ ~\Anaconda3.7\lib\site-packages\scipy\optimize\_basinhopping.py in basinhopping(func, x0, niter, T, stepsize, minimizer_kwargs, take_step, accept_test, callback, interval, disp, niter_success, seed)
180
+
181
+ 667
182
+
183
+ 668 bh = BasinHoppingRunner(x0, wrapped_minimizer, take_step_wrapped,
184
+
185
+ --> 669 accept_tests, disp=disp)
186
+
187
+ 670
188
+
189
+ 671 # start main iteration loop
190
+
191
+
192
+
193
+ ~\Anaconda3.7\lib\site-packages\scipy\optimize\_basinhopping.py in __init__(self, x0, minimizer, step_taking, accept_tests, disp)
194
+
195
+ 72
196
+
197
+ 73 # do initial minimization
198
+
199
+ ---> 74 minres = minimizer(self.x)
200
+
201
+ 75 if not minres.success:
202
+
203
+ 76 self.res.minimization_failures += 1
204
+
205
+
206
+
207
+ ~\Anaconda3.7\lib\site-packages\scipy\optimize\_basinhopping.py in __call__(self, x0)
208
+
209
+ 284 return self.minimizer(x0, **self.kwargs)
210
+
211
+ 285 else:
212
+
213
+ --> 286 return self.minimizer(self.func, x0, **self.kwargs)
214
+
215
+ 287
216
+
217
+ 288
218
+
219
+
220
+
221
+ ~\Anaconda3.7\lib\site-packages\scipy\optimize\_minimize.py in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
222
+
223
+ 593 return _minimize_cg(fun, x0, args, jac, callback, **options)
224
+
225
+ 594 elif meth == 'bfgs':
226
+
227
+ --> 595 return _minimize_bfgs(fun, x0, args, jac, callback, **options)
228
+
229
+ 596 elif meth == 'newton-cg':
230
+
231
+ 597 return _minimize_newtoncg(fun, x0, args, jac, hess, hessp, callback,
232
+
233
+
234
+
235
+ ~\Anaconda3.7\lib\site-packages\scipy\optimize\optimize.py in _minimize_bfgs(fun, x0, args, jac, callback, gtol, norm, eps, maxiter, disp, return_all, **unknown_options)
236
+
237
+ 968 else:
238
+
239
+ 969 grad_calls, myfprime = wrap_function(fprime, args)
240
+
241
+ --> 970 gfk = myfprime(x0)
242
+
243
+ 971 k = 0
244
+
245
+ 972 N = len(x0)
246
+
247
+
248
+
249
+ ~\Anaconda3.7\lib\site-packages\scipy\optimize\optimize.py in function_wrapper(*wrapper_args)
250
+
251
+ 298 def function_wrapper(*wrapper_args):
252
+
253
+ 299 ncalls[0] += 1
254
+
255
+ --> 300 return function(*(wrapper_args + args))
256
+
257
+ 301
258
+
259
+ 302 return ncalls, function_wrapper
260
+
261
+
262
+
263
+ ~\Anaconda3.7\lib\site-packages\scipy\optimize\optimize.py in approx_fprime(xk, f, epsilon, *args)
264
+
265
+ 728
266
+
267
+ 729 """
268
+
269
+ --> 730 return _approx_fprime_helper(xk, f, epsilon, args=args)
270
+
271
+ 731
272
+
273
+ 732
274
+
275
+
276
+
277
+ ~\Anaconda3.7\lib\site-packages\scipy\optimize\optimize.py in _approx_fprime_helper(xk, f, epsilon, args, f0)
278
+
279
+ 668 ei[k] = 1.0
280
+
281
+ 669 d = epsilon * ei
282
+
283
+ --> 670 grad[k] = (f(*((xk + d,) + args)) - f0) / d[k]
284
+
285
+ 671 ei[k] = 0.0
286
+
287
+ 672 return grad
288
+
289
+
290
+
291
+ ValueError: setting an array element with a sequence.
292
+
293
+ ```
294
+
295
+
296
+
297
+
392
298
 
393
299
 
394
300
 
395
301
  ### 補足情報（FW/ツールのバージョンなど）
396
302
 
397
- Python 3
303
+ Python 3(Anaconda)
398
-
399
-
400
-
304
+
305
+
306
+
307
+
308
+
401
- TEST.xlsxの中身
309
+ 前データ.xlsxの中身
402
-
310
+
403
- ![イメージ説明](61ee1191cef36ad744d85996d185fb97.png)
311
+ ![イメージ説明](aaea75f0441d03bbc871790cc02935bb.png)