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

0

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