On Mon, Aug 3, 2009 at 2:36 AM, laurent<lgaut...@gmail.com> wrote:
>
>
> On Wed, 2009-07-29 at 20:38 -0400, Skipper Seabold wrote:
>> On Wed, Jul 29, 2009 at 8:17 PM, Skipper Seabold<jsseab...@gmail.com> wrote:
>> > On Thu, Jul 23, 2009 at 9:50 AM, Skipper Seabold<jsseab...@gmail.com>
>> > wrote:
>> >> On Thu, Jul 23, 2009 at 7:12 AM, Peter<rpy-l...@maubp.freeserve.co.uk>
>> >> wrote:
>> >>> On Thu, Jul 23, 2009 at 12:02 PM, Laurent Gautier<lgaut...@gmail.com>
>> >>> wrote:
>> >>>>>
>> >>>>> Since x is a NumPy array these can come in different flavours, and
>> >>>>> perhaps something here is breaking rpy2 or R, e.g.
>> >>>>
>> >>>> The bug report if for rpy-1.0.x, I think.
>> >>>
>> >>
>> >> In [2]: rpy.__version__
>> >> Out[2]: '1.4.0.dev7055'
>> >>
>> >> Is that expected behavior? That is the numpy release that I currently
>> >> have installed. Note that I tested this on two machines however. But
>> >> yes it's rpy-1.0.3 I believe.
>> >>
>> >>> Sorry, you are right. I would still try making a copy of the NumPy array
>> >>> which should give you a continuous C ordered array which owns its
>> >>> own memory.
>> >>>
>> >>> Peter
>> >>>
>> >>
>> >>
>> >> Thanks, but I think I checked this. When I said exactly the same, as
>> >> far as I can tell they are exactly the same.
>> >>
>> >> Say exog is my x that I have imported and x is built as above from rpy
>> >> and numpy.
>> >>
>> >> In [16]: x.flags
>> >> Out[16]:
>> >> C_CONTIGUOUS : False
>> >> F_CONTIGUOUS : True
>> >> OWNDATA : False
>> >> WRITEABLE : True
>> >> ALIGNED : True
>> >> UPDATEIFCOPY : False
>> >>
>> >> In [17]: exog.flags
>> >> Out[17]:
>> >> C_CONTIGUOUS : False
>> >> F_CONTIGUOUS : True
>> >> OWNDATA : False
>> >> WRITEABLE : True
>> >> ALIGNED : True
>> >> UPDATEIFCOPY : False
>> >>
>> >> In [18]: x.flags == exog.flags
>> >> Out[18]: True
>> >>
>> >> In [19]: x.dtype
>> >> Out[19]: dtype('float64')
>> >>
>> >> In [20]: exog.dtype
>> >> Out[20]: dtype('float64')
>> >>
>> >> In [21]: x.shape
>> >> Out[21]: (21, 4)
>> >>
>> >> In [22]: exog.shape
>> >> Out[22]: (21, 4)
>> >>
>> >> I also tried making a copy of exog, so I had a C contiguous array that
>> >> owns its own data, but it was the same result. The weird thing is
>> >> that using exog works for psi.bisquare and psi.huber, but when I try
>> >> psi.hampel I get an error about the conformability of x and the wts,
>> >> so it appears to be some sort of a dimensions issue in the iterated
>> >> weighted least squares steps for RLM. Anyway, I've update my code to
>> >> use the array as built up from rpy for now and it works. When I have
>> >> time I will try to dig deeper. I might also try to reinstall rpy
>> >> tonight, though again this happened on two machines.
>> >>
>> >> Cheers,
>> >>
>> >> Skipper
>> >>
>> >
>> > For posterity's sake, I will post my solution. The problem was in how
>> > Hampel's norm works and there are two workarounds. The first is to
>> > add an axis to the endogenous variable (as the initial exception said.
>> > I just looked in the wrong place.). So in my example above for my
>> > variables not from R/Rpy, I would have to add an axis to the response
>> > variable so y = y[:,numpy.newaxis] or y = y[:,None]. This needs to be
>> > done because Hampel's norm is a hard redescender for large outliers
>> > (big ones get a weight of zero rather than just downweighted), so the
>> > starting values are chosen differently and it was here that my
>> > approach was failing.
>> >
>> > The second and probably "correct" way to do it (the above lead to some
>> > odd behavior in how the results were returned) is to specify how you
>> > want the initial values to be calculated. So my command is then
>> >
>> > results3 = r.rlm(formula, data=frame, psi='psi.hampel', init='lts') #
>> > lts = Least Trimmed Squares
>> >
>> > and everything works as usual.
>> >
>> > Skipper
>> >
>>
>> I spoke too soon. You still need to add the axis, and the results are
>> returned as dicts rather than arrays in my script though nothing else
>> is different. The weirdest thing has been the consistency with which
>> I can reproduce this...
>
>
> Did you have better luck with rpy2 ?
>
I inherited a wrapper around RPy as part of a test suite, so I didn't
bother changing to RPy2. I have been out of town, but I will see if
the results are different for RPy2 some time this week and if it is
worth changing over.
Skipper
------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now. http://p.sf.net/sfu/bobj-july
_______________________________________________
rpy-list mailing list
rpy-list@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rpy-list
