Integral in GAUSS

There are two things that need to change in your code:

1. Your first input to integrate1d is not correct. It needs to be a pointer to a procedure and can be made like this:

    Answer = integrate1d(&integrand, __INFN, __INFp, yd,xd,thetaNd,hpower);

or this:

    // Create pointer to procedure 'integrand'
    integ = &integrand;

    Answer = integrate1d(integ, __INFN, __INFp, yd,xd,thetaNd,hpower);

1. To create a vector of results, you need to first create a vector that is the size of your final output. The 'zeros' function is a good option for this.

Here is the full code with these changes:

// Minimization Alg. to find theta_hat
proc ( 1 ) = QNminimize(y,x);
    local n,thetaN;
    n = rows(y);
    thetaN = sumc(y .* exp(x) ./ (1+exp(x)))/sumc((exp(x) ./ (1+exp(x)))^2);
    retp(thetaN);
endp;

// TN_integral value
proc ( 1 ) = integrand(u,y,x,thetaN,hpower);
    local n,gam,arg,ker,gxtheta,T_N,T_Ninteg;
    n = rows(y);
    gam = 1/n^hpower;
    arg = (x-u)/gam;
    ker = pdfn(arg);
    gxtheta = thetaN*exp(x) ./ (1+exp(x));
    T_N = ker .* (y-gxtheta);
    T_Ninteg = sumc(T_N)^2 * exp(-u^2/2);
    retp(T_Ninteg);
endp;

/* main program */
library gauss, pgraph;
pqgwin many;

seed1 = 123;
iter = 500; @ # iterations @
nn = 100; @ sample size selections @

hpower = 1/3; @ value of bandwidth @

// Pre-allocate vector to hold the results
Answer = zeros(iter, 1);

k = 1;
do until k > iter; @ loop for simulation @
    xd = rndns(nn,1,seed1);
    yd = rndns(nn,1,seed1);

    // minimization
    thetaNd = QNminimize(yd,xd);

    // Integrand
    Answer[k] = integrate1d(&integrand, __INFN, __INFp, yd,xd,thetaNd,hpower);

    k = k+1;
endo;

I'm not sure if this is exactly what you want, but it should show you the things you need to know to make it work:

X = { 0.2841436,
      1.2426602,
     -0.2392296 };

Y = { -0.08043074,
      -0.25865713,
      -0.82383213 };

/*
** Integrate the procedure 'kernel' from negative to positive infinity    
** __INFN is negative infinity and __INFP is positive infinity
*/
area = integrate1d(&kernel, __INFN, __INFp, X, Y);

proc (1) = kernel(u, X, Y);
    local result;

    result = sumc(pdfn(X - u) .* Y.^2);

    retp(result);
endp;

Thank you for the response. It was very helpful!

I also have an additional question. How can I use this integral inside of a loop where X and Y are different per each iteration.

I know that I should put the proc(1)=kernel(u,X,Y) outside of the loop. Then, how can I put the line called "area" inside of the loop that it could be used for different data set? Please see the small example below, where I have a bit complicated case with minimization step and an additional function inside of "integrand" rather than the kernel function that you used above.

From this example, I need to have 500 Answers of Integrals as output. Thanks in advance for the help!

// Minimization Alg. to find theta_hat
proc ( 1 ) = QNminimize(y,x);
    local n,thetaN;
    n = rows(y);
    thetaN = sumc(y .* exp(x) ./ (1+exp(x)))/sumc((exp(x) ./ (1+exp(x)))^2);
    retp(thetaN);
endp;

// TN_integral value
proc ( 1 ) = integrand(u,y,x,thetaN,hpower);
    local n,gam,arg,ker,gxtheta,T_N,T_Ninteg;
    n = rows(y);
    gam = 1/n^hpower;
    arg = (x-u)/gam;
    ker = pdfn(arg);
    gxtheta = thetaN*exp(x) ./ (1+exp(x));
    T_N = ker .* (y-gxtheta);
    T_Ninteg = sumc(T_N)^2 * exp(-u^2/2);
    retp(T_Ninteg);
endp;

/* main program */
library gauss, pgraph;
pqgwin many;

seed1 = 123;
iter = 500; @ # iterations @
nn = 100; @ sample size selections @

hpower = 1/3; @ value of bandwidth @

k = 1;
do until k > iter; @ loop for simulation @
    xd = rndns(nn,1,seed1);
    yd = rndns(nn,1,seed1);

    //minimization
    thetaNd = QNminimize(yd,xd);

    //Integrand
    Integ = integrand(u,yd,xd,thetaNd,hpower);
    Answer = integrate1d(&Integ, __INFN, __INFp, yd,xd,thetaNd,hpower);

    k = k+1;
endo;