tBuLi / symfit

Symbolic Fitting; fitting as it should be.
http://symfit.readthedocs.org
MIT License
233 stars 17 forks source link

How to maximize likelihood with multiple parameters? #276

Closed lesshaste closed 4 years ago

lesshaste commented 4 years ago

I am trying to maximize a likelihood function which has several parameters. For example:

x = Variable('x')
ps0 = Parameter('ps0')
ps1 = Parameter('ps1')
ps2 = Parameter('ps2')
ps3 = Parameter('ps3')
lambdas0 = Parameter('lambdas0')
lambdas1 = Parameter('lambdas1')
lambdas2 = Parameter('lambdas2')
lambdas3 = Parameter('lambdas3')
model = ps0*lambdas0*exp(-lambdas0*x) + ps1*lambdas1*exp(-lambdas1*x) +ps2*lambdas2*exp(-lambdas2*x)+ps3*lambdas3*exp(-lambdas3*x)

xs = scipy.stats.weibull_min.rvs(loc=0, c=0.5, size=1000)

constraints = [Eq(ps0+ps1+ps2+ps3, 1), Ge(ps0, 0.01), Ge(ps1, 0.01), Ge(ps2, 0.01), Ge(ps3, 0.01), Ge(lambdas0, 0.01), Ge(lambdas1, 0.01), Ge(lambdas2, 0.01), Ge(lambdas3, 0.01)]

fit = Fit(model, xs, constraints=constraints, objective=LogLikelihood)

This works OK but is increasingly awkward to write as the number of ps and lambdas variables increases. How can this code be generalized to work with an arbitrary number of parameters? The number of ps and lambdas variables is always equal.

pckroon commented 4 years ago
lambdas = parameters(', '.join(['lambdas{}'.format(idx) for idx in range(4)]))
ps = parameters(', '.join(['ps{}'.format(idx) for idx in range(4)]))

model = 0
for idx in range(4):
    ps[idx].min = 0.01
    lambdas[idx].min = 0.01  # Better to use bounds than constraints where able.
    model += ps[idx] * lambdas[idx]*exp(-lambdas[idx]*x)

Something like this, you should be able to figure out the rest from here. There's probably an even shorter and better sympy way of writing this, but meh. As a last note, don't do Ge(ps2, 0.01) to bound the value of your parameters. Use bounds instead, those are much easier to solve for numerically.

lesshaste commented 4 years ago

Thanks so much. I really didn't realise string joins would be the answer!

Thanks also for the advice about constraints versus bounds. If that isn't documented in symfit it might be worth adding it as it's not at all obvious why it should make a difference.

pckroon commented 4 years ago

Good idea, we'll definitely look into it :)