[Rcpp-devel] Using sugar expressions for evaluating deviance residuals

[Douglas Bates]

On Mon, Aug 16, 2010 at 11:47 AM, Romain Francois
<romain at r-enthusiasts.com> wrote:
> Le 16/08/10 15:33, Douglas Bates a écrit :
>>
>> On Mon, Aug 16, 2010 at 8:10 AM, Romain Francois
>> <romain at r-enthusiasts.com>  wrote:
>>>
>>> Nothing really conclusive here ...
>>>
>>> I've added a preliminary version of SUGAR_BLOCK_3 (only for the case
>>> where
>>> all inputs are vectors) so that one can write :
>>>
>>> incl8<- '
>>> inline double y_log_y(double y, double mu) {
>>>    return y ? y*log(y/mu) : 0.;
>>> }
>>> inline double yMu(double y, double mu) {
>>>    return 2. * (y_log_y(y, mu) + y_log_y(1.-y,1.-mu));
>>> }
>>> double yEta(double y, double eta, double w) {
>>>    return w * yMu(y, 1./(1. + exp(-eta)));
>>> }
>>> SUGAR_BLOCK_3(yeta, ::yEta)
>>> '
>>> code8<- '
>>>    NumericVector ans = yeta(NumericVector(yy), NumericVector(ee),
>>> NumericVector(ww) );
>>>    return ans ;
>>> '
>>>
>>> I get :
>>>
>>> $ Rscript ylogy.R
>>> Le chargement a nécessité le package : methods
>>>             test replications elapsed user.self
>>> 1 fx0(yy, ee, ww)          500   0.953     0.947
>>> 2 fx1(yy, ee, ww)          500   2.576     2.521
>>> 3 fx2(yy, ee, ww)          500   2.887     2.843
>>> 4 fx3(yy, ee, ww)          500   2.898     2.854
>>> 5 fx4(yy, ee, ww)          500   1.273     1.241
>>> 6 fx5(yy, ee, ww)          500   1.502     1.450
>>> 7 fx6(yy, ee, ww)          500   1.408     1.356
>>> 8 fx7(yy, ee, ww)          500   1.303     1.273
>>> 9 fx8(yy, ee, ww)          500   1.452     1.389
>>>
>>> So slightly better than fx5 but not as good as fx4
>>>
>>>
>>> I guess there needs to be some investigation on why fx1, fx2 and fx3
>>> don't
>>> perform "well".
>>
>> To me the interesting question is why fx0 does better than fx4 when my
>> guess is that it should not be as fast.  Is it possible that the level
>> of optimization for the C compiler when compiling R is different from
>> the level of optimization when inline calls the C++ compiler?
>
> Apparently not, with the attached file, which essentially is the code that
> lives in family.c, I can actually get better performance than the code in R:
>
> $ Rscript ylogy9.R
> Le chargement a nécessité le package : methods
>             test replications elapsed user.self
> 1 fx0(yy, ee, ww)          500   0.954     0.947
> 2 fx9(yy, ee, ww)          500   0.943     0.936
>
> Probably just because I cache the symbols rather than let .Call retrieve
> them each time.
>
> So, is there something in family.c you don't do in fx4 ?
>
> Is it this perhaps:
>
> for (i = 0; i < n; i++) {
>        double etai = reta[i], tmp;
>        tmp = (etai < MTHRESH) ? DOUBLE_EPS :
>            ((etai > THRESH) ? INVEPS : exp(etai));
>        rans[i] = x_d_opx(tmp);
>    }

That part is to truncate the range of the eta variable so that the
corresponding mu doesn't get to 0 or 1.  When you get to the extremes
of the interval [0,1] for mu (which here represents a probability) the
increment calculation breaks down on a division by zero.  These days
we skip the check there and apply it later when calculating the
variance.

>
> --
> Romain Francois
> Professional R Enthusiast
> +33(0) 6 28 91 30 30
> http://romainfrancois.blog.free.fr
> |- http://bit.ly/bzoWrs : Rcpp svn revision 2000
> |- http://bit.ly/b8VNE2 : Rcpp at LondonR, oct 5th
> `- http://bit.ly/aAyra4 : highlight 0.2-2
>
>