SciPy.optimize.least_squares() 객관적인 함수 질문

Question

3 개의 알 수없는 매개 변수 a, b 및 c0을 최적화하여 고도의 비선형 함수를 최소화하려고합니다. http://www.dewtronics.com/tutorials/roulette/documents/Roulette_Physik.pdf

내가 참조됩니다 식 (35) 및 (40)에 : 나는 여기

파이썬 3에 카지노 룰렛 공의 일부 지배 방정식을 복제하려고 해요 것은 연구 논문에 대한 링크입니다 종이. 기본적으로 나는 바퀴에서 회전하는 룰렛 공의 스톱 워치 랩 측정을합니다. 각각의 연속 무릎에 대해 비 보수적 인 마찰 세력의 모멘텀 손실로 인해 랩타임이 증가 할 것입니다. 그런 다음이 시간 측정을하고 방정식 (40)의 Levenberg-Marquardt 최소 제곱 법을 사용하여 방정식 (35)을 맞 춥니 다.

내 질문은 두 가지입니다. (1) scipy.optimize.least_squares() 메서드 = 'lm'을 사용하고 있는데 목적 함수를 작성하는 방법을 모르겠습니다.

def fall_time(k,a,b,c0): 

    F = (1/(a * b)) * (c0 - np.arcsinh(c0) * np.exp(a * k * 2 * np.pi)) 

    return F 

def parameter_estimation_function(x0,tk): 

    a = x0[0] 
    b = x0[1] 
    c0 = x0[2] 

    S = 0 

    for i,t in enumerate(tk): 

     k = i + 1 

     S += (t - fall_time(k,a,b,c0))**2 

    return [S,1,1] 

sol = least_squares(parameter_estimation_function,[0.1,0.8,-0.1],args=([tk1]),method='lm',jac='2-point',max_nfev=2000) 

print(sol) 

지금, 문서의 예제에서, 나는 그것을이 방법을 쓴 목적 함수를 본 적이 : 용지에서와 같이 지금은 정확히 기록하는 기능을 가지고있다. 문서에서 목적 함수는 항상 잔차의 제곱이 아닌 잔차를 반환합니다. 또한 문서에서 그들은 절대 합계를 사용하지 않습니다! 그래서 sum과 square가 least_squares()이라는 두포에서 자동으로 처리되는지 궁금합니다.

(2) 아마도 두 번째 질문은 객관적인 기능을 작성하는 방법을 이해하지 못한 결과 일 것입니다. 하지만 여하튼, 알고리즘을 최소로 수렴하는 데 문제가 있습니다. 나는이 점을 잘 알고 있습니다. 왜냐하면 levenberg 알고리즘은 "탐욕스럽고"가장 가까운 미니 마 근처에서 멈추기 때문입니다. 그러나 나는 다른 초기 추측을 할 때 거의 동일한 결과를 수렴 할 수있을 것이라고 생각했습니다. 초기 추측에 약간의 변경을 가하면 다른 징후로 매개 변수 결과가 나타납니다. 또한, 나는 아직 모호하게 수렴 할 수있는 초기 추측의 조합을 찾지 못했습니다! 해결책을 찾기 전에 항상 시간이 초과됩니다. 나는 기능 평가의 양을 10,000으로 늘려서 그것이 가능한지 알아 보았다. 아무 소용이 없다!

아마도 누군가 내 실수를 밝힐 수 있습니다! 나는 아직도 파이썬과 scipy 라이브러리에 비교적 새로운 편이다! 여기

내가 여기 비디오에서 자신을 측정 한 그 tk에 대한 몇 가지 샘플 데이터입니다 : https://www.youtube.com/watch?v=0Zj_9ypBnzg

tk = [0.52,1.28,2.04,3.17,4.53,6.22] 
tk1 = [0.51,1.4,2.09,3,4.42,6.17] 
tk2 = [0.63,1.35,2.19,3.02,4.57,6.29] 
tk3 = [0.63,1.39,2.23,3.28,4.70,6.32] 
tk4 = [0.57,1.4,2.1,3.06,4.53,6.17] 

감사

Answer 1

1) 예, 합 잔차의 제곱을 의심으로 자동으로 처리됩니다.

2) 문제에 깊이 관여하지 않아서 (예 : 얼마나 많은 로컬 미니 마가 존재하는지, '합리적인'결과를 구성하는지 등) 말하기 어렵습니다. 나중에 더 조사 할 수 있습니다.

하지만 킥을 쳐서 어떤 결과가 나오는지 몇 마디 쳐다 보았습니다. 예를 들어, 상수를 독립 변수 b_inv으로 바꿀 수 있습니다. 결과가 상당히 안정적으로 나타납니다. 결과를 확인하는 데 사용한 코드는 다음과 같습니다. (간결성을 위해 목적 함수를 다시 작성 했으므로 전체 결과를 변경하지 않고 numpy 배열의 요소 별 연산을 활용합니다.) 콘솔 출력

import numpy as np 
from scipy.optimize import least_squares 

def fall_time(k,a,b_inv,c0): 
    return (b_inv/a) * (c0 - np.arcsinh(c0) * np.exp(a * k * 2 * np.pi)) 

def parameter_estimation_function(x,tk): 
    return np.asarray(tk) - fall_time(k=np.arange(1,len(tk)+1), a=x[0],b_inv=x[1],c0=x[2]) 

tk_samples = [ 
    [0.52,1.28,2.04,3.17,4.53,6.22], 
    [0.51,1.4,2.09,3,4.42,6.17], 
    [0.63,1.35,2.19,3.02,4.57,6.29], 
    [0.63,1.39,2.23,3.28,4.70,6.32], 
    [0.57,1.4,2.1,3.06,4.53,6.17] 
    ] 

for i in range(len(tk_samples)): 
    sol = least_squares(parameter_estimation_function,[0.1,1.25,-0.1], 
     args=(tk_samples[i],),method='lm',jac='2-point',max_nfev=2000) 
    print(sol.x) 

:

[ 0.03621789 0.64201913 -0.12072879] 
[ 3.59319972e-02 1.17129458e+01 -6.53358716e-03] 
[ 3.55516005e-02 1.48491493e+01 -5.31098257e-03] 
[ 3.18068316e-02 1.11828091e+01 -7.75329834e-03] 
[ 3.43920725e-02 1.25160378e+01 -6.36307506e-03]