Re: [R] [ESS] How to get correct integration in C for step function?

From: Thomas Lumley <tlumley_at_u.washington.edu>
Date: Mon 22 Jan 2007 - 17:42:01 GMT

On Mon, 22 Jan 2007, Lynette wrote:

> Well,
>
> I have no idea either. I can get correct answers for continous functions but
> incorrect for step functions.
>

I have just tried using Rdqags from C for the function x>0 and it worked fine (once I had declared all the arguments correctly). The code is below.

         -thomas

#include "Rinternals.h"
#include "R_ext/Applic.h"

double stepfn(double x){

    return (x>0);
}

SEXP call_stepfn(SEXP x){
SEXP answer;
int i,n;

n=length(x);
PROTECT(answer=allocVector(REALSXP,n));
for(i=0;i<n;i++){

    REAL(answer)[i]=stepfn(REAL(x)[i]);
}
UNPROTECT(1);
return answer;
}

void vec_stepfn(double *x, int n, void *ex){ int i;
for (i=0;i<n;i++) x[i]=stepfn(x[i]);
return;
}

void Cvec_stepfn(double *x, double *n){

  vec_stepfn(x, *n, (void *) NULL);
  return;
}

SEXP int_stepfn(SEXP lower, SEXP upper){

SEXP answer;
double result, abserr;
int last, neval, ier;
int lenw;
int *iwork;
double *work;
int limit=100;
double reltol=0.00001;
double abstol=0.00001;

          lenw = 4 * limit;
          iwork =   (int *) R_alloc(limit, sizeof(int));
          work = (double *) R_alloc(lenw,  sizeof(double));

          Rdqags(vec_stepfn, (void *) NULL, REAL(lower), REAL(upper),
                        &abstol,  &reltol,
                       &result,  &abserr,  &neval,  &ier,
                        &limit,  &lenw, &last,
                        iwork, work);

printf("%f %f %d\n", result,abserr, ier);

    PROTECT(answer=allocVector(REALSXP,1));     REAL(answer)[0]=result;
    UNPROTECT(1);     return answer;
}



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Tue Jan 23 04:48:26 2007

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Mon 22 Jan 2007 - 19:31:49 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.