我正在用Fipy解一个偏微分方程系统,它包含两个参数或常数,所以我想知道在Fipy中是否也可以估计这些参数,或者其他库中的哪些更适合于这一点。你知道吗
注意:我知道scipy有一些功能(优化.最小化对于MLE),但我不确定将它们应用于Fipy代码是否足够。你知道吗
更新:对于下面的PDE系统,我想估计两个未知参数:“Beta”和“m”
在Fipy中求解这个偏微分方程的函数如下:
import scipy as sci
import fipy as fipy
import numpy as np
from fipy import *
# Grid
nx = 100
ny = 100
dx = 1.
dy = dx
mesh = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy)
x = mesh.cellCenters[0]
y = mesh.cellCenters[1]
# Setting variable of results and adding inicial conditions
u = CellVariable(name="Individual 1", mesh=mesh, value=0.)
u.setValue(1., where=(50. < x) & (70. > x) & (50. < y) & (70. > y))
v = CellVariable(name="Individual 2", mesh=mesh, value=0.)
v.setValue(1., where=(40. < x) & (60. > x) & (40. < y) & (60. > y))
p = CellVariable(name= "Marks Individual 1", mesh=mesh, value=0.)
p.setValue(1., where=(50. < x) & (70. > x) & (50. < y) & (70. > y))
q = CellVariable(name= "Marks Individual 2", mesh=mesh, value=0.)
q.setValue(1., where=(40. < x) & (60. > x) & (40. < y) & (60. > y))
# Plotting inicial conditions
if __name__ == '__main__':
viewer = Viewer(u, v, datamin=0., datamax=1.)
viewer.plot()
# Setting PDE
def HRMLE(params):
m = params[0]
beta = params[1]
D = 1.
CU = CellVariable(mesh=mesh, rank=1)
CU[:]= 1.
CU.setValue(-1., where = (x > 60.) * [[[1], [0]]])
CU.setValue(-1., where = (y > 60.) * [[[0], [1]]])
CV = CellVariable(mesh=mesh, rank=1)
CV[:]=1.
CV.setValue(-1., where = (x > 50.) * [[[1], [0]]])
CV.setValue(-1., where = (y > 50.) * [[[0], [1]]])
# Transient formulation
eqU = TransientTerm() == DiffusionTerm(coeff=D) -\
ConvectionTerm(coeff=CU*q.value*beta)
eqV = TransientTerm() == DiffusionTerm(coeff=D) -\
ConvectionTerm(coeff=CV*p.value*beta)
eqP = TransientTerm() == u*(1 + m*q) - p
eqQ = TransientTerm() == v*(1 + m*p) - q
# Solving Transient term
timeStepDuration = 1.
steps = 50
t = timeStepDuration * steps
for step in range(steps):
eqU.solve(var=u, dt=timeStepDuration)
eqV.solve(var=v, dt=timeStepDuration)
eqP.solve(var=p, dt=timeStepDuration)
eqQ.solve(var=q, dt=timeStepDuration)
# Plotting results
#if __name__ == '__main__':
# vieweru = Viewer(u, datamin=0., datamax=1.)
# viewerv = Viewer(v, datamin=0., datamax=1.)
# vieweru.plot()
# viewerv.plot()
loglink = np.sum(np.log(u.value)) + np.sum(np.log(v.value))
return(loglink)
最后,对于最大可能性标准,我想最大化:
设置初始值,并使用scipy
mb = [0., .5]
mb
results = sci.optimize.minimize(HRMLE, mb, method='Nelder-Mead')
results
结果显示,参数的值总是接近我的初始值,这就是我不确定我的过程是否正确的原因。如有任何建议,我们将不胜感激。你知道吗
好吧,我已经意识到,为了使函数最大化,必须将函数的输出乘以-1,然后将其最小化。所以我的函数的输出应该是:
另一方面,实际上最大化只针对网格中的点列表,而不是所有的u1和u2值。你知道吗
相关问题 更多 >
编程相关推荐