Integration over an interval

Suppose f(.) is the standard normal PDF and x follows a truncated standard normal distribution over [0,1]. I want to find E(x) = int x f(x)/int f(x).  So I used the integration functions in GAUSS 16 to do it but ran into problems.  A minimal program to replicate the problem can be found below.  Obviously GAUSS can integrate f(x) over [0,1], the denominator. But GAUSS seems to have trouble with the numerator int x f(x) over [0,1] (if you uncomment the four commands in the program and run it you will see the error messages). Mathematically I can't see any problem with the numerator (it's not like we have 0 in the denominator or something like that). Matlab can handle it. Any pointer to a solution? Is it something wrong in my program? Thanks!

new;

e_x_theory1 = integrate1d(&stdnormal,0,1);
e_x_theory2 = intquad1(&stdnormal,1|0);

"Theoretical value based on intergrate1d" e_x_theory1;
"Theoretical value based on intquad1" e_x_theory2;

/*
e_x_theory3 = integrate1d(&xstdn,0,1);
e_x_theory4 = intquad1(&xstdn,1|0);

"Theoretical value based on intergrate1d" e_x_theory3;
"Theoretical value based on intquad1" e_x_theory4;
*/

end;

proc (1) = stdnormal(x);
retp(1/(sqrt(2*pi))*exp(-0.5*x^2));
endp;

proc (1) = xstdn(x);
retp(1/(sqrt(2*pi))*exp(-0.5*x^2)*x);
endp;

2 Answers



1



accepted

When I uncomment the second set of calculations in the code that you posted, I get the error:

G0036 : Matrices are not conformable

The problem I see is that both integrate1d and intquad require the function to be able to take in a vector of input values and return a vector of output values. Your stdnormal procedure is capable of this, for example:

x = { 0, 1, 2 };

print stdnormal(x);

proc (1) = stdnormal(x);
    retp(1/(sqrt(2*pi))*exp(-0.5*x^2));
endp;

will return:

0.39894228 
0.24197072 
0.053990967

because all of the multiplications are either (scalar * vector) or (scalar * scalar). However, your procedure xstdn multiplies two vectors, using the '*' operator which indicates a matrix multiply. If you change the final multiply to a '.*', then the code will run.

x = { 0, 1, 2 };

print xstdn(x);

proc (1) = xstdn(x);
    //Notice the '.*' in the next line
    retp(1/(sqrt(2*pi))*exp(-0.5*x^2).*x);
endp;

Should return to us:

0.0000000 
0.24197072 
0.10798193

aptech

1,323


0



Yes, you are right. Because x^2 is actually an element by element square, I forgot to .* x in the second function. Thanks!

Your Answer


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.

Try GAUSS for 30 days for FREE

See what GAUSS can do for your data

© Aptech Systems, Inc. All rights reserved.

Privacy Policy