 # Integral in GAUSS

I would like to take an integral where the main function contains data (x_i,y_i). The length of the data is N. The integrand is the (summation of kernel times y_i)^2. The kernel function is denoted by $K[(x_i-u)/h]$ which is the kernel of standard Normal, and h is the known bandwidth, 0.3 (say).

We need to take the integral over u and the bounds of integral are  -infinity and +infinity. The form of integration is as follows:

$\int_{-\infty}^{+\infty}{\sum_{i=1}^{N} K[(x_i-u)/h]*y_i}^2 du$

I do not know how to code this in GAUSS. Any guidance would be appreciated.

I can also put a simple data for (x_i,y_i) here:

Xk        Yk
0.2841436 -0.08043074
1.2426602 -0.25865713
-0.2392296 -0.82383213

The standard Normal kernel is as follows:

fn gauss(r) = 1/(2*sqrt(pi)).*exp(-(r^2)/2);

0

accepted

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

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; aptech

1,728

0

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; aptech

1,728

0

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; 0
accepted

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

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; aptech
1,728
0

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; aptech
1,728
0

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; ### Have a Specific Question?

Get a real answer from a real person

### Need Support?

Get help from our friendly experts.