<p>
        </p>
        thanks both of you,<div><br></div><div><span id="goog_1125882228"></span>Josef, the data that I sent is only the first 100 rows of about 1500, there should be sufficient sampling in each plot.</div><div><br></div><div><span id="goog_1125882232"></span><span id="goog_1125882233"></span><span id="goog_1125882234"></span><span id="goog_1125882235"></span>Skipper, I have attempted to deploy your suggestion for not linearizing the data. It seems to work. I'm a little confused at your modification if the getDiam function and I wonder if you could help me understand. The form of the equation that is being fit is:</div><div><br></div><div><span id="goog_1125882236"></span><span id="goog_1125882237"></span><span id="goog_1125882238"></span><span id="goog_1125882239"></span>Y= a*X^b</div><div><br></div><div><span id="goog_1125882240"></span>your version of the detDaim function:<span id="goog_1125882246"></span><span id="goog_1125882247"></span></div><div><blockquote type="cite" style="border-left-style: solid; border-left-color: rgb(204, 204, 204); border-top-width: 1px; border-right-width: 1px; border-bottom-width: 1px; border-left-width: 1px; margin-left: 0px; padding-left: 10px; padding-right: 0px; margin-right: 0px; "><div><div><div><blockquote type="cite" style="border-left-style: solid; border-left-color: rgb(204, 204, 204); border-top-width: 1px; border-right-width: 1px; border-bottom-width: 1px; border-left-width: 1px; margin-left: 0px; padding-left: 10px; padding-right: 0px; margin-right: 0px; "><div><jsseabold@gmail.com><ptittmann@gmail.com><br>def getDiam(ht, *b):<br>&nbsp; &nbsp;return ht[:,0]**b[0] + np.sum(b[1:]*ht[:,1:], axis=1)</ptittmann@gmail.com></jsseabold@gmail.com></div></blockquote></div></div></div></blockquote><div>Im sorry if this is an obvious question but I don't understand how this works as it seems that the "a" coefficient is missing.</div><div><br></div><div><span id="goog_1125882252"></span><span id="goog_1125882253"></span><span id="goog_1125882254"></span><span id="goog_1125882255"></span>Thanks again!&nbsp;</div></div><div><span id="goog_1125882244"></span><span id="goog_1125882245"></span>&nbsp;<br>
        <br>
        <div id="CCA740C6DAFC4A38900F6CA94C003987">--&nbsp;<br>Peter Tittmann<div><br></div></div>
                <br>
                
        <p>On Thursday, December 2, 2010 at 5:03 PM, josef.pktd@gmail.com wrote:</p>
        <blockquote type="cite" style="border-left-style:solid;border-width:1px;margin-left:0px;padding-left:10px;">
            <div><div><div>On Thu, Dec 2, 2010 at 7:57 PM, Skipper Seabold &lt;<a href="mailto:jsseabold@gmail.com">jsseabold@gmail.com</a>&gt; wrote:<br><blockquote type="cite"><div> On Thu, Dec 2, 2010 at 7:11 PM, Skipper Seabold <jsseabold@gmail.com> wrote:<br>&gt; On Thu, Dec 2, 2010 at 5:59 PM, Peter Tittmann <ptittmann@gmail.com> wrote:<br>&gt;&gt; getDiam is a predictor to get dbh from height. It works with curve_fit to<br>&gt;&gt; find coefficients a and b given datasetset of known dbh/height pairs. You<br>&gt;&gt; are right, what I want is dummy variables for each plot. I'll see if I can<br>&gt;&gt; get that worked out by revising getDiam..<br>&gt;&gt; Thanks again<br>&gt;&gt;<br>&gt;<br>&gt; I think it would be easier to create your dummy variables before you pass it in.<br>&gt;<br>&gt; You might find some of the tools in statsmodels to be helpful here.<br>&gt; We don't yet have an ANCOVA model, but you could definitely do<br>&gt; something like the following. &nbsp;Not sure if it's exactly what you want,<br>&gt; but it should give you an idea.<br>&gt;<br>&gt; import numpy as np<br>&gt; import scikits.statsmodels as sm<br>&gt;<br>&gt; dta = np.genfromtxt('./db_out.csv', delimiter=",", names=True, dtype=None)<br>&gt; plot_dummies, col_map = sm.tools.categorical(dta['plot'], drop=True,<br>&gt; dictnames=True)<br>&gt;<br>&gt; plot_dummies will be dummy variables for all of the "plot" categories,<br>&gt; and col_map is a map from the column number to the plot just so you<br>&gt; can be sure you know what's what.<br>&gt;<br>&gt; I don't see how to use your objective function though with dummy<br>&gt; variables. &nbsp;What happens if the effect of one of the plots is<br>&gt; negative, then you run into 0 ** -1.5 == inf.<br>&gt;<br><br> If you want to do NLLS and not linearize then something like this<br> might work and still keep the dummy variables as shift parameters<br><br> def getDiam(ht, *b):<br> &nbsp; &nbsp;return ht[:,0]**b[0] + np.sum(b[1:]*ht[:,1:], axis=1)<br><br> X = np.column_stack((indHtPlot, plot_dummies))<br> Y = depDbh<br> coefs, cov = optimize.curve_fit(getDiam, X, Y, p0= [0.]*X.shape[1])<br></ptittmann@gmail.com></jsseabold@gmail.com></div></blockquote><br>In the sample file there are 11 levels of the `plot` that have only a<br>single observation each. I tried to use onewaygls, but statsmodels.OLS<br>doesn't work if y is a scalar.<br><br>I don't know whether curvefit or optimize.leastsq will converge in<br>this case, good starting values might be necessary.<br><br>Josef<br><br><br><blockquote type="cite"><div><br><br>&gt; You could linearize your objective function to be<br>&gt;<br>&gt; b*ln(ht)<br>&gt;<br>&gt; and do something like<br>&gt;<br>&gt; indHtPlot = dta['height']<br>&gt; depDbh = dta['dbh']<br>&gt; X = np.column_stack((np.log(indHtPlot), plot_dummies))<br>&gt; Y = np.log(depDbh)<br>&gt; res = sm.OLS(Y,X).fit()<br>&gt; res.params<br>&gt; array([ 0.98933264, -1.35239293, -1.0623305 , -0.99155293, -1.33675099,<br>&gt; &nbsp; &nbsp; &nbsp; -1.30657011, -1.50933751, -1.28744779, -1.43937358, -1.33805883,<br>&gt; &nbsp; &nbsp; &nbsp; -1.32744257, -1.42672539, -1.35239293, -1.60585046, -1.45239093,<br>&gt; &nbsp; &nbsp; &nbsp; -1.45695112, -1.34811186, -1.32658794, -1.21721715, -1.32853084,<br>&gt; &nbsp; &nbsp; &nbsp; -1.45775017, -1.44460388, -2.19065236, -1.3303631 , -1.20509831,<br>&gt; &nbsp; &nbsp; &nbsp; -1.37341535, -1.25746105, -1.33954972, -1.33922709, -1.247304 &nbsp;])<br>&gt;<br>&gt; Note that your coefficient on height is now an elasticity. &nbsp;I'm sure<br>&gt; I'm missing something here, but that might help you along the way.<br>&gt;<br>&gt; Skipper<br>&gt;<br> _______________________________________________<br> SciPy-User mailing list<br> <a href="mailto:SciPy-User@scipy.org">SciPy-User@scipy.org</a><br> <a href="http://mail.scipy.org/mailman/listinfo/scipy-user">http://mail.scipy.org/mailman/listinfo/scipy-user</a><br><br></div></blockquote>_______________________________________________<br>SciPy-User mailing list<br><a href="mailto:SciPy-User@scipy.org">SciPy-User@scipy.org</a><br>http://mail.scipy.org/mailman/listinfo/scipy-user<br></div></div></div>
                        
                        
                        
                        
        </blockquote>
                
                <br>
    

</div>