【问题标题】:Equation of the parabola enclosing a scatter plot包含散点图的抛物线方程
【发布时间】:2018-11-11 07:46:45
【问题描述】:

我将下图的 x 和 y 坐标作为两个不同的列表 s 和 k。如果我想确定最能包围这些点的抛物线方程,我该怎么做?

谢谢!

我不确定如何链接到外部数据集。

What is the best way to add your dataset to a question?

这是我能找到的最好的。所以,我在这里建立了一个样本数据集。不过,我稍后会删除它。

                K             S
0    9.500000e-01 -6.500000e-01
1    8.500000e-01 -6.000000e-01
2    9.000000e-01 -6.000000e-01
3    9.500000e-01 -6.000000e-01
4    7.000000e-01 -5.500000e-01
5    7.500000e-01 -5.500000e-01
6    8.000000e-01 -5.500000e-01
7    8.500000e-01 -5.500000e-01
8    9.000000e-01 -5.500000e-01
9    9.500000e-01 -5.500000e-01
10   6.000000e-01 -5.000000e-01
11   6.500000e-01 -5.000000e-01
12   7.000000e-01 -5.000000e-01
13   7.500000e-01 -5.000000e-01
14   8.000000e-01 -5.000000e-01
15   8.500000e-01 -5.000000e-01
16   9.000000e-01 -5.000000e-01
17   9.500000e-01 -5.000000e-01
18   5.000000e-01 -4.500000e-01
19   5.500000e-01 -4.500000e-01
20   6.000000e-01 -4.500000e-01
21   6.500000e-01 -4.500000e-01
22   7.000000e-01 -4.500000e-01
23   7.500000e-01 -4.500000e-01
24   8.000000e-01 -4.500000e-01
25   8.500000e-01 -4.500000e-01
26   9.000000e-01 -4.500000e-01
27   9.500000e-01 -4.500000e-01
28   4.000000e-01 -4.000000e-01
29   4.500000e-01 -4.000000e-01
30   5.000000e-01 -4.000000e-01
31   5.500000e-01 -4.000000e-01
32   6.000000e-01 -4.000000e-01
33   6.500000e-01 -4.000000e-01
34   7.000000e-01 -4.000000e-01
35   7.500000e-01 -4.000000e-01
36   8.000000e-01 -4.000000e-01
37   8.500000e-01 -4.000000e-01
38   9.000000e-01 -4.000000e-01
39   9.500000e-01 -4.000000e-01
40   3.000000e-01 -3.500000e-01
41   3.500000e-01 -3.500000e-01
42   4.000000e-01 -3.500000e-01
43   4.500000e-01 -3.500000e-01
44   5.000000e-01 -3.500000e-01
45   5.500000e-01 -3.500000e-01
46   6.000000e-01 -3.500000e-01
47   6.500000e-01 -3.500000e-01
48   7.000000e-01 -3.500000e-01
49   7.500000e-01 -3.500000e-01
50   8.000000e-01 -3.500000e-01
51   8.500000e-01 -3.500000e-01
52   9.000000e-01 -3.500000e-01
53   9.500000e-01 -3.500000e-01
54   2.000000e-01 -3.000000e-01
55   2.500000e-01 -3.000000e-01
56   3.000000e-01 -3.000000e-01
57   3.500000e-01 -3.000000e-01
58   4.000000e-01 -3.000000e-01
59   4.500000e-01 -3.000000e-01
60   5.000000e-01 -3.000000e-01
61   5.500000e-01 -3.000000e-01
62   6.000000e-01 -3.000000e-01
63   6.500000e-01 -3.000000e-01
64   7.000000e-01 -3.000000e-01
65   7.500000e-01 -3.000000e-01
66   8.000000e-01 -3.000000e-01
67   8.500000e-01 -3.000000e-01
68   9.000000e-01 -3.000000e-01
69   9.500000e-01 -3.000000e-01
70   1.000000e-01 -2.500000e-01
71   1.500000e-01 -2.500000e-01
72   2.000000e-01 -2.500000e-01
73   2.500000e-01 -2.500000e-01
74   3.000000e-01 -2.500000e-01
75   3.500000e-01 -2.500000e-01
76   4.000000e-01 -2.500000e-01
77   4.500000e-01 -2.500000e-01
78   5.000000e-01 -2.500000e-01
79   5.500000e-01 -2.500000e-01
80   6.000000e-01 -2.500000e-01
81   6.500000e-01 -2.500000e-01
82   7.000000e-01 -2.500000e-01
83   7.500000e-01 -2.500000e-01
84   8.000000e-01 -2.500000e-01
85   8.500000e-01 -2.500000e-01
86   9.000000e-01 -2.500000e-01
87   9.500000e-01 -2.500000e-01
88   5.000000e-02 -2.000000e-01
89   1.000000e-01 -2.000000e-01
90   1.500000e-01 -2.000000e-01
91   2.000000e-01 -2.000000e-01
92   2.500000e-01 -2.000000e-01
93   3.000000e-01 -2.000000e-01
94   3.500000e-01 -2.000000e-01
95   4.000000e-01 -2.000000e-01
96   4.500000e-01 -2.000000e-01
97   5.000000e-01 -2.000000e-01
98   5.500000e-01 -2.000000e-01
99   6.000000e-01 -2.000000e-01
100  6.500000e-01 -2.000000e-01
101  7.000000e-01 -2.000000e-01
102  7.500000e-01 -2.000000e-01
103  8.000000e-01 -2.000000e-01
104  8.500000e-01 -2.000000e-01
105  9.000000e-01 -2.000000e-01
106  9.500000e-01 -2.000000e-01
107  8.881784e-16 -1.500000e-01
108  5.000000e-02 -1.500000e-01
109  1.000000e-01 -1.500000e-01
110  1.500000e-01 -1.500000e-01
111  2.000000e-01 -1.500000e-01
112  2.500000e-01 -1.500000e-01
113  3.000000e-01 -1.500000e-01
114  3.500000e-01 -1.500000e-01
115  4.000000e-01 -1.500000e-01
116  4.500000e-01 -1.500000e-01
117  5.000000e-01 -1.500000e-01
118  5.500000e-01 -1.500000e-01
119  6.000000e-01 -1.500000e-01
120  6.500000e-01 -1.500000e-01
121  7.000000e-01 -1.500000e-01
122  7.500000e-01 -1.500000e-01
123  8.000000e-01 -1.500000e-01
124  8.500000e-01 -1.500000e-01
125  9.000000e-01 -1.500000e-01
126  9.500000e-01 -1.500000e-01
127 -5.000000e-02 -1.000000e-01
128  8.881784e-16 -1.000000e-01
129  5.000000e-02 -1.000000e-01
130  1.000000e-01 -1.000000e-01
131  1.500000e-01 -1.000000e-01
132  2.000000e-01 -1.000000e-01
133  2.500000e-01 -1.000000e-01
134  3.000000e-01 -1.000000e-01
135  3.500000e-01 -1.000000e-01
136  4.000000e-01 -1.000000e-01
137  4.500000e-01 -1.000000e-01
138  5.000000e-01 -1.000000e-01
139  5.500000e-01 -1.000000e-01
140  6.000000e-01 -1.000000e-01
141  6.500000e-01 -1.000000e-01
142  7.000000e-01 -1.000000e-01
143  7.500000e-01 -1.000000e-01
144  8.000000e-01 -1.000000e-01
145  8.500000e-01 -1.000000e-01
146  9.000000e-01 -1.000000e-01
147  9.500000e-01 -1.000000e-01
148 -1.000000e-01 -5.000000e-02
149 -5.000000e-02 -5.000000e-02
150  8.881784e-16 -5.000000e-02
151  5.000000e-02 -5.000000e-02
152  1.000000e-01 -5.000000e-02
153  1.500000e-01 -5.000000e-02
154  2.000000e-01 -5.000000e-02
155  2.500000e-01 -5.000000e-02
156  3.000000e-01 -5.000000e-02
157  3.500000e-01 -5.000000e-02
158  4.000000e-01 -5.000000e-02
159  4.500000e-01 -5.000000e-02
160  5.000000e-01 -5.000000e-02
161  5.500000e-01 -5.000000e-02
162  6.000000e-01 -5.000000e-02
163  6.500000e-01 -5.000000e-02
164  7.000000e-01 -5.000000e-02
165  7.500000e-01 -5.000000e-02
166  8.000000e-01 -5.000000e-02
167  8.500000e-01 -5.000000e-02
168  9.000000e-01 -5.000000e-02
169  9.500000e-01 -5.000000e-02
170 -1.500000e-01  8.881784e-16
171 -1.000000e-01  8.881784e-16
172 -5.000000e-02  8.881784e-16
173  8.881784e-16  8.881784e-16
174  5.000000e-02  8.881784e-16
175  1.000000e-01  8.881784e-16
176  1.500000e-01  8.881784e-16
177  2.000000e-01  8.881784e-16
178  2.500000e-01  8.881784e-16
179  3.000000e-01  8.881784e-16
180  3.500000e-01  8.881784e-16
181  4.000000e-01  8.881784e-16
182  4.500000e-01  8.881784e-16
183  5.000000e-01  8.881784e-16
184  5.500000e-01  8.881784e-16
185  6.000000e-01  8.881784e-16
186  6.500000e-01  8.881784e-16
187  7.000000e-01  8.881784e-16
188  7.500000e-01  8.881784e-16
189  8.000000e-01  8.881784e-16
190  8.500000e-01  8.881784e-16
191  9.000000e-01  8.881784e-16
192  9.500000e-01  8.881784e-16
193 -2.000000e-01  5.000000e-02
194 -1.500000e-01  5.000000e-02
195 -1.000000e-01  5.000000e-02
196 -5.000000e-02  5.000000e-02
197  8.881784e-16  5.000000e-02
198  5.000000e-02  5.000000e-02
199  1.000000e-01  5.000000e-02
200  1.500000e-01  5.000000e-02
201  2.000000e-01  5.000000e-02
202  2.500000e-01  5.000000e-02
203  3.000000e-01  5.000000e-02
204  3.500000e-01  5.000000e-02
205  4.000000e-01  5.000000e-02
206  4.500000e-01  5.000000e-02
207  5.000000e-01  5.000000e-02
208  5.500000e-01  5.000000e-02
209  6.000000e-01  5.000000e-02
210  6.500000e-01  5.000000e-02
211  7.000000e-01  5.000000e-02
212  7.500000e-01  5.000000e-02
213  8.000000e-01  5.000000e-02
214  8.500000e-01  5.000000e-02
215  9.000000e-01  5.000000e-02
216  9.500000e-01  5.000000e-02
217 -2.000000e-01  1.000000e-01
218 -1.500000e-01  1.000000e-01
219 -1.000000e-01  1.000000e-01
220 -5.000000e-02  1.000000e-01
221  8.881784e-16  1.000000e-01
222  5.000000e-02  1.000000e-01
223  1.000000e-01  1.000000e-01
224  1.500000e-01  1.000000e-01
225  2.000000e-01  1.000000e-01
226  2.500000e-01  1.000000e-01
227  3.000000e-01  1.000000e-01
228  3.500000e-01  1.000000e-01
229  4.000000e-01  1.000000e-01
230  4.500000e-01  1.000000e-01
231  5.000000e-01  1.000000e-01
232  5.500000e-01  1.000000e-01
233  6.000000e-01  1.000000e-01
234  6.500000e-01  1.000000e-01
235  7.000000e-01  1.000000e-01
236  7.500000e-01  1.000000e-01
237  8.000000e-01  1.000000e-01
238  8.500000e-01  1.000000e-01
239  9.000000e-01  1.000000e-01
240  9.500000e-01  1.000000e-01
241 -1.500000e-01  1.500000e-01
242 -1.000000e-01  1.500000e-01
243 -5.000000e-02  1.500000e-01
244  8.881784e-16  1.500000e-01
245  5.000000e-02  1.500000e-01
246  1.000000e-01  1.500000e-01
247  1.500000e-01  1.500000e-01
248  2.000000e-01  1.500000e-01
249  2.500000e-01  1.500000e-01
250  3.000000e-01  1.500000e-01
251  3.500000e-01  1.500000e-01
252  4.000000e-01  1.500000e-01
253  4.500000e-01  1.500000e-01
254  5.000000e-01  1.500000e-01
255  5.500000e-01  1.500000e-01
256  6.000000e-01  1.500000e-01
257  6.500000e-01  1.500000e-01
258  7.000000e-01  1.500000e-01
259  7.500000e-01  1.500000e-01
260  8.000000e-01  1.500000e-01
261  8.500000e-01  1.500000e-01
262  9.000000e-01  1.500000e-01
263  9.500000e-01  1.500000e-01
264 -1.000000e-01  2.000000e-01
265 -5.000000e-02  2.000000e-01
266  8.881784e-16  2.000000e-01
267  5.000000e-02  2.000000e-01
268  1.000000e-01  2.000000e-01
269  1.500000e-01  2.000000e-01
270  2.000000e-01  2.000000e-01
271  2.500000e-01  2.000000e-01
272  3.000000e-01  2.000000e-01
273  3.500000e-01  2.000000e-01
274  4.000000e-01  2.000000e-01
275  4.500000e-01  2.000000e-01
276  5.000000e-01  2.000000e-01
277  5.500000e-01  2.000000e-01
278  6.000000e-01  2.000000e-01
279  6.500000e-01  2.000000e-01
280  7.000000e-01  2.000000e-01
281  7.500000e-01  2.000000e-01
282  8.000000e-01  2.000000e-01
283  8.500000e-01  2.000000e-01
284  9.000000e-01  2.000000e-01
285  9.500000e-01  2.000000e-01
286 -5.000000e-02  2.500000e-01
287  8.881784e-16  2.500000e-01
288  5.000000e-02  2.500000e-01
289  1.000000e-01  2.500000e-01
290  1.500000e-01  2.500000e-01
291  2.000000e-01  2.500000e-01
292  2.500000e-01  2.500000e-01
293  3.000000e-01  2.500000e-01
294  3.500000e-01  2.500000e-01
295  4.000000e-01  2.500000e-01
296  4.500000e-01  2.500000e-01
297  5.000000e-01  2.500000e-01
298  5.500000e-01  2.500000e-01
299  6.000000e-01  2.500000e-01
300  6.500000e-01  2.500000e-01
301  7.000000e-01  2.500000e-01
302  7.500000e-01  2.500000e-01
303  8.000000e-01  2.500000e-01
304  8.500000e-01  2.500000e-01
305  9.000000e-01  2.500000e-01
306  9.500000e-01  2.500000e-01
307  8.881784e-16  3.000000e-01
308  5.000000e-02  3.000000e-01
309  1.000000e-01  3.000000e-01
310  1.500000e-01  3.000000e-01
311  2.000000e-01  3.000000e-01
312  2.500000e-01  3.000000e-01
313  3.000000e-01  3.000000e-01
314  3.500000e-01  3.000000e-01
315  4.000000e-01  3.000000e-01
316  4.500000e-01  3.000000e-01
317  5.000000e-01  3.000000e-01
318  5.500000e-01  3.000000e-01
319  6.000000e-01  3.000000e-01
320  6.500000e-01  3.000000e-01
321  7.000000e-01  3.000000e-01
322  7.500000e-01  3.000000e-01
323  8.000000e-01  3.000000e-01
324  8.500000e-01  3.000000e-01
325  9.000000e-01  3.000000e-01
326  9.500000e-01  3.000000e-01
327  5.000000e-02  3.500000e-01
328  1.000000e-01  3.500000e-01
329  1.500000e-01  3.500000e-01
330  2.000000e-01  3.500000e-01
331  2.500000e-01  3.500000e-01
332  3.000000e-01  3.500000e-01
333  3.500000e-01  3.500000e-01
334  4.000000e-01  3.500000e-01
335  4.500000e-01  3.500000e-01
336  5.000000e-01  3.500000e-01
337  5.500000e-01  3.500000e-01
338  6.000000e-01  3.500000e-01
339  6.500000e-01  3.500000e-01
340  7.000000e-01  3.500000e-01
341  7.500000e-01  3.500000e-01
342  8.000000e-01  3.500000e-01
343  8.500000e-01  3.500000e-01
344  9.000000e-01  3.500000e-01
345  9.500000e-01  3.500000e-01
346  1.500000e-01  4.000000e-01
347  2.000000e-01  4.000000e-01
348  2.500000e-01  4.000000e-01
349  3.000000e-01  4.000000e-01
350  3.500000e-01  4.000000e-01
351  4.000000e-01  4.000000e-01
352  4.500000e-01  4.000000e-01
353  5.000000e-01  4.000000e-01
354  5.500000e-01  4.000000e-01
355  6.000000e-01  4.000000e-01
356  6.500000e-01  4.000000e-01
357  7.000000e-01  4.000000e-01
358  7.500000e-01  4.000000e-01
359  8.000000e-01  4.000000e-01
360  8.500000e-01  4.000000e-01
361  9.000000e-01  4.000000e-01
362  9.500000e-01  4.000000e-01
363  2.000000e-01  4.500000e-01
364  2.500000e-01  4.500000e-01
365  3.000000e-01  4.500000e-01
366  3.500000e-01  4.500000e-01
367  4.000000e-01  4.500000e-01
368  4.500000e-01  4.500000e-01
369  5.000000e-01  4.500000e-01
370  5.500000e-01  4.500000e-01
371  6.000000e-01  4.500000e-01
372  6.500000e-01  4.500000e-01
373  7.000000e-01  4.500000e-01
374  7.500000e-01  4.500000e-01
375  8.000000e-01  4.500000e-01
376  8.500000e-01  4.500000e-01
377  9.000000e-01  4.500000e-01
378  9.500000e-01  4.500000e-01
379  3.000000e-01  5.000000e-01
380  3.500000e-01  5.000000e-01
381  4.000000e-01  5.000000e-01
382  4.500000e-01  5.000000e-01
383  5.000000e-01  5.000000e-01
384  5.500000e-01  5.000000e-01
385  6.000000e-01  5.000000e-01
386  6.500000e-01  5.000000e-01
387  7.000000e-01  5.000000e-01
388  7.500000e-01  5.000000e-01
389  8.000000e-01  5.000000e-01
390  8.500000e-01  5.000000e-01
391  9.000000e-01  5.000000e-01
392  9.500000e-01  5.000000e-01
393  4.000000e-01  5.500000e-01
394  4.500000e-01  5.500000e-01
395  5.000000e-01  5.500000e-01
396  5.500000e-01  5.500000e-01
397  6.000000e-01  5.500000e-01
398  6.500000e-01  5.500000e-01
399  7.000000e-01  5.500000e-01
400  7.500000e-01  5.500000e-01
401  8.000000e-01  5.500000e-01
402  8.500000e-01  5.500000e-01
403  9.000000e-01  5.500000e-01
404  9.500000e-01  5.500000e-01
405  5.000000e-01  6.000000e-01
406  5.500000e-01  6.000000e-01
407  6.000000e-01  6.000000e-01
408  6.500000e-01  6.000000e-01
409  7.000000e-01  6.000000e-01
410  7.500000e-01  6.000000e-01
411  8.000000e-01  6.000000e-01
412  8.500000e-01  6.000000e-01
413  9.000000e-01  6.000000e-01
414  9.500000e-01  6.000000e-01
415  6.000000e-01  6.500000e-01
416  6.500000e-01  6.500000e-01
417  7.000000e-01  6.500000e-01
418  7.500000e-01  6.500000e-01
419  8.000000e-01  6.500000e-01
420  8.500000e-01  6.500000e-01
421  9.000000e-01  6.500000e-01
422  9.500000e-01  6.500000e-01
423  7.500000e-01  7.000000e-01
424  8.000000e-01  7.000000e-01
425  8.500000e-01  7.000000e-01
426  9.000000e-01  7.000000e-01
427  9.500000e-01  7.000000e-01
428  8.500000e-01  7.500000e-01
429  9.000000e-01  7.500000e-01
430  9.500000e-01  7.500000e-01

【问题讨论】:

  • 这只是一个理论,但你可以在抛物线上取离散点,首先让它变得非常平坦。的,增加梯度untik一定数量的点在一定的边界内。
  • 您能提供一个指向您的数据的链接吗?
  • 我已经建立了一个样本数据集。没有更简单的方法吗?必须有。
  • 当然。您可以链接到 dtpaste 之类的网站。无论如何,我正在使用我自己的示例数据来解决您的问题。
  • 困难的部分实际上是定义“最佳”。我看不出有一个自然的定义。此外,根据您的数据源,您还可以考虑放宽“封闭”部分(例如,来自物理测量的数据可能会产生异常值)。

标签: python python-3.x curve-fitting


【解决方案1】:

试试这个代码。它通过反复试验逼近抛物线。

import matplotlib.pyplot as plt
import numpy as np
s = list([0,1,2,3,4,5,6,7,8])
k = list([0,1,4,9,16,25,36,49,64]) #Insert Data Here(positive side only)

def parab(x):
    for f in range(0,len(s)): 
        v = f**2/(x/2)
        yield v

def parab2():
    t = 201
    for i in range(0,200):
        g = 0
        t-=1
        foo = np.fromiter(parab(t), dtype = 'int')
        for bar in range(0,len(s)):
            if foo[bar] < k[bar]+0.5 and foo[bar] > k[bar]-0.5:
                g += 1
        if g >= len(s)/2:
            print(foo)
            print("x**2/",t/2)

parab2()

试试看。如果它近似于您,请尝试摆弄 x^2 除以的内容。

此外,随着数据集大小的增加,请尝试增加该行中的值 if g&gt;=4:

您可以决定哪种尺寸最适合您。

编辑更新代码:效果更好。

【讨论】:

    【解决方案2】:

    此解决方案需要一些迭代,但似乎大部分时间都有效。然而,在同样的情况下,它会被抛物线卡在错误的一侧,有时开口似乎太宽了。但是,如果输入数据与 OP 图像所暗示的一样密集(与我的随机输入数据相比),这可能是可以避免的。

    该示例使用a + b * x**2 类型的抛物线,但添加线性项没有问题。

    代码如下:

    import matplotlib.pyplot as plt
    import numpy as np
    from random import random
    from scipy.optimize import leastsq, curve_fit, fmin
    from scipy.spatial import ConvexHull
    import warnings
    np.set_printoptions( precision = 5 )
    
    def parabola( x, a, b, c ):
        return a + b * x + c * x**2
    
    ### find the shortest distance between a given point and a parabola
    ### actually give x such as ( x, f(x) ) is closest to ( x0, y0 )
    def orth_dist( x0, y0, a, b, c ):
        f = lambda x: np.sqrt( ( x - x0 )**2 + ( parabola( x, a, b, c ) - y0 )**2 )
        sol = fmin( f, x0, disp=False )
        return sol[ 0 ]
    
    ### cost function: standard distance if point is above parabola
    ### or extra cost if point is below
    def myDist( r, dy , off=4e0 ):
        out =  abs( r )
        if dy <= 0:
            out = max( [out, off + off * abs( r)  ] )
        return out
    
    def residuals( p, data, off=0 ): # data is of type [ [x0,y0], [x1,y1],... ] 
        xList, yList = zip( *data )
        ### a standard fit uses the distance in y direction between curve and 
        ### data point. Here it is better to use the shortest distance.
        ### To find the shortest distance some extra numerics is required.
        xopt = [ orth_dist( x0, y0, *p ) for x0, y0 in zip( xList, yList ) ]
        yopt = [ parabola( x, *p ) for x in xopt ]
        yFitList = [ parabola( x, *p ) for x in xList ]
        rList = [ np.sqrt( ( x - x0 )**2 + (y - y0 )**2 ) for x, x0, y , y0 in zip( xList, xopt, yList, yopt ) ]
        signList = [ y - yTh  for y, yTh in zip( yList, yFitList ) ]
        out = [ myDist( r, dy , off=off ) for r, dy in zip( rList, signList ) ]
        return out
    
    ### removing the upper points of the convex hull
    ### as it is convex, it can be at most a straight line 
    ### from the most left to the most right point
    ### if a point of the hull is on or above this line, dump it
    def remove_upper( data ):
        sortedByX = sorted( data)
        xList, yList = zip( *sortedByX )
        m = ( yList[-1] - yList[0] ) / ( xList[-1] - xList[0] )
        out = [ sortedByX[0] ]
        for p in range( 1, len(data ) ):
            x = xList[ p ]
            y = yList[ p ]
            if y < m * ( x - xList[0]) + yList[0]:
                out += [ sortedByX[ p ] ]
        out += [ sortedByX[ -1 ] ]
        return out
    
    ### my choice of 'true' parameters
    myParameters = ( .45, -1.3, 2.33)
    
    ### just some random data for testing
    scatterData = list()
    for i in range(180):
        x = 3*( 2 * random() - 1) + 1.1
        y = 8 * random()
        if y > parabola( x, *myParameters ):
            scatterData += [ [ x, y ] ]
    
    ### to fit the parabula we do not need inner points, 
    ### so let's take the convex hull
    hull = ConvexHull(scatterData )
    scatterDataHull = list()
    for index in hull.vertices:
        scatterDataHull += [ scatterData[ index ] ]
    
    ### the upper points of the convex hull are also useless 
    #### if not disturbing, so get rid of them as well 
    scatterDataHullR = remove_upper( scatterDataHull )
    
    ### exact data for plotting
    xxList = np.linspace( -4, 4, 300 )
    yyList = [ parabola( x, *myParameters ) for x in xxList ]
    
    scatterX, scatterY =zip(*scatterDataHullR)
    
    start, pcov = curve_fit( parabola,  *zip(*scatterDataHullR) )
    
    fit = list( start )
    counter = 0
    test = True
    myOff = .1
    while test:
        counter += 1
        ### "simple" fit with orthogonal distance
        fit, err = leastsq( residuals, x0=fit, args=( scatterDataHullR, myOff ) )
        ### as points below the parabola are forbidden, we give those a high
        ### cost in the fitting. Starting with very high costs, however, 
        ### results in bad convergence. So we start with a rather low cost,
        ### check if some points are below and---if so---increase the cost.
        ### We then fit again with the previous result as starting point.
        ### (Note sometimes this fails, but I did not have time to
        ### investigate this cases. Moreover I hope that this is due 
        ### to my quite random points, while the points in the OP look 
        ### by for more regular.)
        print "fit",fit 
        xList, yList = zip( *scatterDataHullR )
        yFitList = [ parabola( x, *fit ) for x in xList ]
        signList = [ y - yTh  for y, yTh in zip( yList, yFitList ) ]
        print counter, myOff
        print signList
        if min( signList ) < 0:
            myOff *= 2
        else:
            test = False
        if counter > 50: # maxiter
            test = False
            warnings.warn( "Max iterations reached.", UserWarning )
    
    startList = [ parabola( x, *start ) for x in xxList ]
    fitList = [ parabola( x, *fit ) for x in xxList ]
    
    fig = plt.figure()
    ax = fig.add_subplot( 1, 1, 1 )
    ax.plot( xxList, yyList , label='true parabula')
    ax.plot( xxList, startList, label='simple fit' )
    ax.plot( xxList, fitList, label='result')
    ax.scatter( *zip( *scatterData ), s=60, label='full data')
    ax.scatter( *zip( *scatterDataHull ), s=40, label='convex hull')
    ax.scatter( *zip( *scatterDataHullR ), s=20, label='lower hull' )
    ax.legend( loc=0 )
    ax.set_xlim( [ -2, 5 ] )
    ax.set_ylim( [ 0, 10 ] )
    plt.show()
    

    输出是:

    >>fit [ 0.62572 -1.75611  2.81445]
    >>1 0.1
    >>[0.6213711303862084, 0.18684681979766804, -0.30967646606163246, -0.146185783279136, 1.5115384666142262e-08, 0.04091503213347902, -0.22456258154469477, -0.10366901944547369]
    >>fit [ 0.6484  -1.80986  2.75176]
    >>2 0.2
    >>[0.62199186259362, 0.17515271586535786, -0.34273136096659584, -0.17923227341075365, 9.718855165097295e-10, 0.057645756995909014, -0.10363664705794662, 0.13933840419257315]
    >>fit [ 0.6484  -1.80986  2.75176]
    >>3 0.4
    >>[0.62199186259362, 0.17515271586535786, -0.34273136096659584, -0.17923227341075365, 9.718855165097295e-10, 0.057645756995909014, -0.10363664705794662, 0.13933840419257315]
    >>fit [ 0.53489 -1.56037  2.54129]
    >>4 0.8
    >>[1.3184522612014575, 0.7674472388911582, -0.021666294928952112, 0.025470161835318805, 0.056401815047789794, 0.10025191012633528, 2.6445583856116173e-07, 0.4250810523461581]
    >>fit [ 0.51566 -1.53254  2.53166]
    >>5 1.6
    >>[1.3836247162551256, 0.8254553204503852, 0.016315167888327764, 0.0536843248273029, 0.06793421458085641, 0.10848235358916569, 2.42328743738085e-08, 0.42471967128564003]
    

    【讨论】:

    • 您能再解释一下您的代码吗?它看起来确实很整洁。我第一次看到很多东西。
    • @hegdep 我在代码中添加了几个 cmets。我希望这会有所帮助。如果您还有其他问题,请告诉我。
    • 我现在更好地遵循曲线拟合部分。但是,如何确定抛物线的初始方程为parabola( x, .45, 2.33 )?就我而言,我有两个列表,我想确定抛物线的初始方程。
    • @hegdep 好的,就我而言,这只是我的随机选择,我坚持下去。在您的情况下,我选择了使抛物线独一无二的 3 个点;说最左边,最右边和最低点(希望这些条件不重合)。由此你可以确定一条起始抛物线。
    • @hegdep 你也可以从减少的凸包中取出最左边和最右边的第三个。我实际上会对减少的凸包做一个简单的curve_fit。如果值在 1 的范围内,这似乎是这种情况,它应该可以在没有初始猜测的情况下工作。然后,您可以将该拟合的结果作为迭代的起点。
    猜你喜欢
    • 2011-07-02
    • 1970-01-01
    • 2015-08-17
    • 1970-01-01
    • 1970-01-01
    • 2021-01-28
    • 2020-08-02
    • 1970-01-01
    • 2011-07-14
    相关资源
    最近更新 更多