Optimising the generation of a large number of random numbers using python 3
Question: I am wanting to generate eight random numbers within a range (0 to pi/8), add
them together, take the sine of this sum, and after doing this N times, take
the mean result. After scaling this up I get the correct answer, but it is too
slow for `N > 10^6`, especially when I am averaging over N trials `n_t = 25`
more times! I am currently getting this code to run in around _12 seconds_ for
`N = 10^5`, meaning that it will take _20 minutes_ for `N = 10^7`, which
doesn't seem optimal (it may be, I don't know!).

My code is as follows:

    
    
    import random
    import datetime
    from numpy import pi
    from numpy import sin
    import numpy
    t1 = datetime.datetime.now()
    
    def trial(N):
        total = []
        uniform = numpy.random.uniform
        append = total.append
        for j in range(N):
            sum = 0
            for i in range (8):
                sum+= uniform(0, pi/8)
            append(sin(sum))
        return total
    
    N = 1000000
    n_t = 25
    total_squared = 0
    ans = []
    for k in range (n_t):
        total = trial(N)
        f_mean = (numpy.sum(total))/N
        ans.append(f_mean*((pi/8)**8)*1000000)
    sum_square = 0
    for e in ans:
        sum_square += e**2
    sum = numpy.sum(ans)
    mean = sum/n_t
    variance = sum_square/n_t - mean**2
    s_d = variance**0.5
    print (mean, " ± ", s_d)
    t2 = datetime.datetime.now()
    print ("Execution time: %s" % (t2-t1))
    

If anyone can help me optimise this it would be much appreciated!

Thank you :)


Answer: Given your requirement of obtaining the result with this method,
`np.sin(np.random.uniform(0,np.pi/8,size=(8,10**6,25)).sum(axis=0)).mean(axis=0)`
gets you your 25 trials pretty quickly... This is fully vectorised (and
concise which is always a bonus!) so I doubt you could do any better...

Explanation:

You generate a big random 3d array of size `(8 x 10**6 x 25)`. `.sum(axis=0)`
will get you the sum over the first dimension (`8`). `np.sin(...)` applies
elementwise. `.mean(axis=0)` will get you the mean over the first remaining
dimension (`10**6`) and leave you with a 1d array of length (`25`)
corresponding to your trials.

