[SciPy-User] ancova with optimize.curve_fit

[Skipper Seabold]

On Thu, Dec 2, 2010 at 4:10 PM, Peter Tittmann <ptittmann at gmail.com> wrote:
> Thanks very much for your reply Josef, here are the tracebacks from both the
> original and from your suggestion:
> with: curve_fit(getDiam, indHtPlot, depDbh)
> In [20]: ancova()
> ------------------------------------------------------------
> Traceback (most recent call last):
>   File "<ipython console>", line 1, in <module>
>   File "<ipython console>", line 5, in ancova
>   File "/usr/local/lib/python2.6/dist-packages/scipy/optimize/minpack.py",
> line 422, in curve_fit
>     res = leastsq(func, p0, args=args, full_output=1, **kw)
>   File "/usr/local/lib/python2.6/dist-packages/scipy/optimize/minpack.py",
> line 273, in leastsq
>     m = check_func(func,x0,args,n)[0]
>   File "/usr/local/lib/python2.6/dist-packages/scipy/optimize/minpack.py",
> line 13, in check_func
>     res = atleast_1d(thefunc(*((x0[:numinputs],)+args)))
>   File "/usr/local/lib/python2.6/dist-packages/scipy/optimize/minpack.py",
> line 343, in _general_function
>     return function(xdata, *params) - ydata
> ValueError: shape mismatch: objects cannot be broadcast to a single shape
>
> with curve_fit(getDiam, indHtPlot.T, depDbh)
> In [22]: ancova()
> Error in sys.excepthook:
> Traceback (most recent call last):
>   File "/usr/lib/pymodules/python2.6/IPython/iplib.py", line 2028, in
> excepthook
>     self.showtraceback((etype,value,tb),tb_offset=0)
>   File "/usr/lib/pymodules/python2.6/IPython/iplib.py", line 1729, in
> showtraceback
>     self.InteractiveTB(etype,value,tb,tb_offset=tb_offset)
>   File "/usr/lib/pymodules/python2.6/IPython/ultraTB.py", line 998, in
> __call__
>     print >> out, self.text(etype, evalue, etb)
>   File "/usr/lib/pymodules/python2.6/IPython/ultraTB.py", line 1012, in text
>     return FormattedTB.text(self,etype,value,tb,context=5,mode=mode)
>   File "/usr/lib/pymodules/python2.6/IPython/ultraTB.py", line 937, in text
>     if len(elist) > self.tb_offset:
> TypeError: object of type 'NoneType' has no len()
> Original exception was:
> ValueError: object too deep for desired array
> ------------------------------------------------------------
> Traceback (most recent call last):
>   File "/usr/lib/pymodules/python2.6/IPython/iplib.py", line 2257, in
> runcode
>     exec code_obj in self.user_global_ns, self.user_ns
>   File "<ipython console>", line 1, in <module>
>   File "<ipython console>", line 5, in ancova
>   File "/usr/local/lib/python2.6/dist-packages/scipy/optimize/minpack.py",
> line 422, in curve_fit
>     res = leastsq(func, p0, args=args, full_output=1, **kw)
>   File "/usr/local/lib/python2.6/dist-packages/scipy/optimize/minpack.py",
> line 281, in leastsq
>     maxfev, epsfcn, factor, diag)
> error: Result from function call is not a proper array of floats.
>
> --
> Peter Tittmann
>
> On Thursday, December 2, 2010 at 1:01 PM, josef.pktd at gmail.com wrote:
>
> On Thu, Dec 2, 2010 at 2:55 PM, Peter Tittmann <ptittmann at gmail.com> wrote:
>
> Greetings,
> Im attempting to conduct analysis of covariance (ANCOVA) using a non-linear
> regression with curve_fit in optimize. The doc string states:
> "xdata : An N-length sequence or an (k,N)-shaped array for functions with k
> predictors. The independent variable where the data is measured."
> I am hoping that this means that if I pass the independent variable and a
> categorical variable, the resulting covariance matrix will reflect the
> variability in the equation coefficients among the categorical variables.
> 1. Is tis the case?
> 2. If so, i'm having a problem with the input array for xdata. The following
> extracts data from a relational database (thats the sql). the eqCoeff()
> function works fine, however when I add a second dimension to the xdata in
> th ancova() function (indHtPlot), curve fit produces an error which seems to
> be related to the structure of my input array. I've tried column_stack and
> vstack to form the arrays. Any assistance would be gratefully received.
>
> import birdseye_db as db
> import numpy as np
> from scipy.optimize import curve_fit
> def getDiam(ht, a, b):
>     dbh = a * ht**b
>     return dbh
> def eqCoeff():
>     '''estimates coefficients a and b in dbh= a* h**b using all trees where
> height was measured'''
>     species=[i[0].strip(' ') for i in db.query('select distinct species from
> plots')]
>     res3d=db.query('select dbh, height, species from plots where ht_code=1')
>     indHt=[i[1] for i in res3d]
>     depDbh=[i[0] for i in res3d]
>     estimated_params, err_est = curve_fit(getDiam, indHt, depDbh)
>     return estimated_params, err_est
>
> def ancova():
>     res=db.query('select dbh, height, plot, species from plots where
> ht_code=1')
>     indHtPlot= np.column_stack(([i[1] for i in res],[i[2] for i in res] ))
>     depDbh=[i[0] for i in res]
>     estimated_params, err_est = curve_fit(getDiam, indHtPlot, depDbh)
>     return estimated_params, err_est
>
>
> Can you post the actual traceback?
>
> My first guess, but I have to look it up, is that you need to
> transpose xdata, e.g. indHtPlot.T
>
> curve_fit(getDiam, indHtPlot.T, depDbh)
>
> Josef
>
>
> Thanks in advance
> --
> Peter Tittmann
>

Can you post a small sample of your data to replicate?

Skipper