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

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

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

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

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

```0.0000000
0.24197072
0.10798193``` aptech

1,728

0

Yes, you are right. Because x^2 is actually an element by element square, I forgot to .* x in the second function. Thanks! 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;``````

```0.0000000
0.24197072
0.10798193``` aptech
1,728
0

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

Get a real answer from a real person

### Need Support?

Get help from our friendly experts.