Jacobi growth process (hasn't been tested this rigorously, but seems to work, up to the answer being transposed and in reverse order in Julia compared to Matlab.)

Author: jiahao

demos/jacobigrowth2.jl

@@ -0,0 +1,13 @@
+#jacobigrowth2.jl
+
+using Distributions
+function jacobigrowth2(m1,m2,p,beta)
+   CJ=eye(p)
+   SJ=zeros(p, p)
+   for k=1:m2
+       Ju = sqrt(rand(Chisq(beta),1,p))
+       Ju *= sqrt( rand(Chisq(beta*p)) / rand(Chisq(beta*(k+m1-p))) ) / norm(Ju)
+       _, _, _, CJ, SJ=svd(CJ,[SJ; Ju])
+   end
+   diag(CJ)
+end