[SciPy-User] ancova with optimize.curve_fit
josef.pktd@gmai...
josef.pktd@gmai...
Thu Dec 2 15:01:19 CST 2010
On Thu, Dec 2, 2010 at 2:55 PM, Peter Tittmann <ptittmann@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
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
More information about the SciPy-User
mailing list