Average size of dot products, Textbook problem

Question

This problem comes from Gilbert Strang's Linear Algebra 5th Edition pg 21, problem 34:

"Using v=randn(3,1) in MATLAB, create a random unit vector u= v / ||v||. Using V = randn(3,30) create 30 more random unit vectors U_j. What is the average size of the dot products |u$\cdot$U_j|? In calculus, the average is $\int_{0}^{\pi} |cos\theta|d\theta/\pi$ = 2/$\pi$."

I used the free alternative, Octave. Here are two implementations. After running them over and over in a 3rd implementation, I begin to see that the average of the averages hovers around .5 not 2/$\pi$

The answers: http://math.mit.edu/~gs/linearalgebra/ila_sol5_ch01.pdf

Please advise.

# my first attempt to calculate one mean
v=randn(3,1)
u=v/norm(v)
V=randn(3,100)
U=randn(3,100) 
for i=1:100
  U(:,i) = V(:,i)/norm(V(:,i))
endfor
udotU = randn(1,100) 
for i=1:100
  udotU(i) = abs(dot(u,U(:,i)))
endfor
mean(udotU)

#second attempt to calculate one mean with specific code snippets given in book
v=randn(3,1)
u=v/norm(v)
V=randn(3,100)
U=randn(3,100) 
D = sqrt (diag (V'* V ))
for i=1:100
  U(:,i) = V(:,i)/D(i)
endfor
udotU = randn(1,100) 
udotU = abs(u'*U)
mean(udotU)

# calculate 50 runs
means=0
for i=1:50
  v=randn(3,1);
  u=v/norm(v);
  V=randn(3,100);
  U=randn(3,100) ;
  D = sqrt (diag (V'* V ));
  for i=1:100
    U(:,i) = V(:,i)/D(i);
  endfor
  udotU = randn(1,100) ;
  udotU = abs(u'*U);
  means=means+mean(udotU)
endfor
means/50

Answer 1

It seems to me that the mean would be $2/\pi$ in $\Bbb R^2$, but your calculations are in $\Bbb R^3$ (as suggested in the textbook answer key you linked to... perhaps this is a typo?). In $\Bbb R^3$, we would need a double integral, for example in spherical co-ordinates, with the calculus approach. Try changing v = randn(3, 1) to v = randn(2, 1), and similarly elsewhere, and you should get $0.636\ldots \approx 2/\pi$.

As a side note, your code is a bit confusing: first, you are initializing U and udotU to random vectors but then overwriting them immediately. Second, rethink your variable names: instead of u = v/norm(v), how about V = v/norm(v) (in which case you should be careful to rename everything else consequently). Finally, it's a bit strange to run through 50 loops when you could achieve the same thing by making your matrix 50 times bigger; i.e., get rid of the loop and use V = randn(3, 5000) instead of V = randn(3, 100).