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
Your Answer
2 Answers
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