[SciPy-User] ODR fitting several equations to the same parameters

josef.pktd@gmai... josef.pktd@gmai...
Thu Nov 12 10:55:40 CST 2009


On Thu, Nov 12, 2009 at 10:04 AM, ms <devicerandom@gmail.com> wrote:
> josef.pktd@gmail.com ha scritto:
>> On Wed, Nov 11, 2009 at 11:26 AM, ms <devicerandom@gmail.com> wrote:
>>> Let's start with a simple example. Imagine I have several linear data
>>> sets y=ax+b which have different b (all of them are known) but that
>>> should fit to the same (unknown) a. To have my best estimate of a, I
>>> would want to fit them all together. In this case it is trivial, you
>>> just subtract the known b from the data set and fit them all at the same
>>> time.
>>>
>>> In my case it is a bit different, in the sense that I have to do
>>> conceptually the same thing but for a highly non-linear equation where
>>> the equivalent of "b" above is not so simple to separate. I wonder
>>> therefore if there is a way to do a simultaneous fit of different
>>> equations differing only in the known parameters and having a single
>>> output, possibly with the help of ODR. Is this possible? And/or what
>>> should be the best thing to do, in general, for this kind of problems?
>>
>> I don't know enough about ODR, but for least squares, optimize.leastsq
>> or curve_fit, it seems you can just substitute any known parameters
>> into your equation.
>>
>> y_i = f(x_i, a, b_i) for each group i
>> plug in values for all b_i, gives reduced f(x_i, a) independent of
>> specific parameters
>> stack equations [y_i for all i] and [f(..) for all i]
>>
>> If you fit this in curve_fit you could also choose the weights, in
>> case the error variance differs by groups.
>>
>> Does this work or am I missing the point?
>
> Probably it's me missing it. Do you just mean to fit them all together
> separately and then make a weighted average of the fitted parameters,
> and using the standard deviation of the mean as the error of the fit? I
> am confused.

I meant stacking all equations into one big estimation problem y =
f(x,a) and minimize squared residual over all equations.
This assumes homoscedastic errors (identical noise in each equation).

an example
(quickly written and not optimized, there are parts I don't remember
about curve_fit, fixed parameters could be better handled by a class)

####################
"""stack equations with different known parameters

I didn't get curve_fit to work with only 1 parameter to estimate

Created on Thu Nov 12 11:17:21 2009
Author: josef-pktd
"""
import numpy as np
from scipy import optimize


def fsingle(a,c,b,x):
    return b*x**a + c

atrue = 1.
ctrue = 10.
b = np.array([[1.]*10, [2.]*10, [3.]*10])
b = np.array([1.,2.,3.])
x = np.random.uniform(size=(3,10))
y = np.hstack([fsingle(atrue, ctrue, b[i], x[i]) for i in range(x.shape[0])])
y += 0.1*np.random.normal(size=y.shape)

def fun(x,a,c):
    #b is taken from enclosing scope
    #print x.shape
    xx=x.reshape((3,10))
    return np.hstack([fsingle(a, c, b[i], xx[i]) for i in range(xx.shape[0])])

res = optimize.curve_fit(fun,x.ravel(),y, p0=np.array([2.,1.]))

print 'true parameters   ', atrue, ctrue
print 'parameter estimate', res[0]
print 'standard deviation', np.sqrt(np.diag(res[1]))
####################




>
> sorry,
> m.
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>


More information about the SciPy-User mailing list