[R] [ESS] How to get correct integration in C for step function?
Thomas Lumley
tlumley at u.washington.edu
Mon Jan 22 18:42:01 CET 2007
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;
}
More information about the R-help
mailing list