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

3 Answers



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

aptech

1,773


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,773


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;

Your Answer

3 Answers

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

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;

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;


You must login to post answers.

Have a Specific Question?

Get a real answer from a real person

Need Support?

Get help from our friendly experts.