【发布时间】:2021-12-13 10:59:57
【问题描述】:
我一直试图找到一种有效的方法来最大化四个变量中的以下怪物函数,但程序运行需要很长时间,我什至不确定结果是否正确。谁能帮我用 Python 更好地编写代码? 函数如下:
在哪里
a=[p,q,r,s]。
Y 是在 30 个点采样的测量数据。
这是我的代码。
import numpy as np
import math
Y=Y_t #Y_t is a predefined column vector with 30 entries.
tstep=0.05 #in s
N=30
cov=np.zeros([30,30])
def R(p,q,r,t):
om_D=p*np.sqrt(1-q**2)
return np.pi*r*(np.exp(-q*p*abs(t)))*(np.cos(om_D*t)+(q/(np.sqrt(1-q**2)))*(np.sin(om_D*abs(t))))/(2*q*(p**3))
def I(m,p):
if m==p:
return 1
else:
return 0
def func(a):
a1=a[0] #natural angular frequency bounds=[3,20]
a2=a[1] #damping ratio bounds=[0,1]
a3=a[2] #psd of forcing signal bounds=[300,600]
a4=a[3] #variance of noise bounds=[0,0.0001] in m
#assuming uniform prior for a, we only have to maximise the likelihood function
for i in range(30):
for j in range(30):
cov[i,j]+=R(a1,a2,a3,(j-i)*tstep)+a4*I(i,j)
P=((2*np.pi)**(-N/2)) * ((np.linalg.det(cov))**(-0.5)) * np.exp((-0.5) *np.linalg.multi_dot([np.transpose(Y),np.linalg.inv(cov),Y]))
return (-1)*P[0]
a_start=[5,0.05,100,0.00001]
bnds=((5,20),(0,1),(300,600),(0,0.0001))
result=spo.differential_evolution(func,bounds=bnds)
print(result.x) ```
【问题讨论】:
-
有趣的问题。但是,我建议对符号进行更多详细说明。 Γ 似乎是一个协方差矩阵。什么是 R_x[ (p-m)Δt |一种] ?条件概率?
-
您可能需要阅读 toeplitz 矩阵,因为这就是您的协方差,有特别有效的例程可以对它们进行反转和分解。
-
Γ 是协方差矩阵; |Γ|是它的行列式,R_x 是自相关函数。
-
我建议你让这个例子可以重现。第二件事 - 您是在寻求另一种方法(更好的优化算法)还是只想加快代码速度、提高性能?
-
@dankal444 如果有更好的优化算法,请分享。该程序仅需要 30 个数据点即可运行半小时。在未来的实践中,我必须增加数据点和变量的数量。所以,它不应该给我带来麻烦。
标签: python math matrix optimization scipy