Estimating the parameters of a PDEs system using Fipy

Question

I am solving a system of PDEs using Fipy which involve two parameters or constants, so I would like to know whether is also possible to estimate those parameters in Fipy, or what others libraries would be more adequate for that.

Note: I know that scipy has some functions for that (optimize.minimize for MLE), but I am not sure if is adequate apply them to a code of Fipy.

UPDATE: For the below PDE system I want to estimate the two unknown parameters: "Beta" and "m"

The function for solving this PDE in Fipy would be something like that:

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)

Finally, for Maximun likelihood criteria I would like to maximise:

Setting an initial value, and using scipy

mb = [0., .5]
mb

results = sci.optimize.minimize(HRMLE, mb, method='Nelder-Mead')
results

The results show values that always are close to my initial values for the parameters, that's the reason I am not sure the my procedure is right. Any suggestions will be quite appreciated.

This is an extremely open ended question and would be improved by having a specific example. — wd15, Mar 01 '17 at 15:23
It might be a good idea to view the "Beta" and "m" landscape before embarking on using the optimizer to check that the solution is as you expect and that the landscape looks sensible. Try a very coarse grid over the "Beta", "m" domain depending on the wall clock duration of the simulations. — wd15, Mar 02 '17 at 16:01

score 0 · Accepted Answer · answered Apr 12 '17 at 14:27

Well, I have realised that for maximise the function is necessary multiply by -1 the output of the function and then minimise it. So the output of my function should be:

loglink = np.sum(np.log(u.value)) + np.sum(np.log(v.value))*-1
return(loglink)

On the other hand, actually maximisation is just for a list of points in the grid, and not for all values of u1 and u2.

Estimating the parameters of a PDEs system using Fipy

1 Answers1