restart with(Statistics): XG := h -> RandomVariable(Normal(0, h)) Zio2I0kiaEc2IkYlNiRJKW9wZXJhdG9yR0YlSSZhcnJvd0dGJUYlLV9JK1N0YXRpc3RpY3NHNiQlKnByb3RlY3RlZEdJKF9zeXNsaWJHRiVJL1JhbmRvbVZhcmlhYmxlR0YlNiMtSSdOb3JtYWxHRiU2JCIiITkkRiVGJUYl MaxIt := 1000 IiUrNQ== randomize(): # Be careful: h must be strictly positive h := evalf([seq(1-k/MaxIt, k=0..MaxIt-1)]): h,h[-1]; NiQkIiIiIiIhJCIrKysrKzUhIzc= # do you really need to print the results ? # # for i from 1 to MaxIt do # printf("%a\134n%a\134n", Sample(XG(h[i]), 5), i); # end do: # Maybe this could suit you? # Stage 1 ; generate a random matrix by sampling Normal(0, 1) M := Sample(XG(1), [MaxIt, 5], method = [envelope, range = -1 .. 1]): # Watch out: stage 2 is computationnaly expensive. # It is presented here to help understanding the method # # Stage 2: rescale M according to the standard deviations 1, 0.999, ...0.001 S := LinearAlgebra:-DiagonalMatrix(h): # Stage 3: do the product S.M Res := S.M LUknTWF0cml4RzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliRzYiNiMvSSQlaWRHRiciNXklKSlRO3lTdVklPQ== # The same thing in a more efficient way M := CodeTools:-Usage( Sample(XG(1), [MaxIt, 5], method = [envelope, range = -1 .. 1]) ): Res := CodeTools:-Usage( < seq(M[n, ..] *~ h[n], n=1..MaxIt) > ); memory used=0.73MiB, alloc change=0 bytes, cpu time=35.00ms, real time=11.00ms, gc time=0ns memory used=0.60MiB, alloc change=0 bytes, cpu time=9.00ms, real time=3.00ms, gc time=0ns LUknTWF0cml4RzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliRzYiNiMvSSQlaWRHRiciNU1CQ18jeVN1WSU9 for i from 1 to MaxIt do # printf("%3d : %a\134n", i, [entries(Res[i,..], nolist)]); end do: