You may want to use MTEX http://mtex-toolbox.github.io/ . There you can generate, e.g. 1000 random rotations. with a single command

rot = rotation.rand(1000)

As a bonus you can generate random samples with respect to any distribution in orientation space.

Euler space is not metric and a random distribution of Euler angles does not generate a random orientation distribution. Only a very simple example which should show you the problem: If you would rotate phi1,90,phi2 with phi1=-phi2 you would always describe the same orientation. 

The best way is to use quaternions. See e.g. 

Orientations and Rotations:

Computations in Crystallographic Textures

of Morawiec, Adam

 I'm not a specialist in crystals.

I'm engaged in gyroscopes and inertial navigation.

In my area, the problem is formulated as follows: there is a random angular velocity vector (due to noise of gyroscopes). And the question: To what random orientation of the object does this lead?

There is another task: there are random errors of orientation, to what errors of navigation does this lead?

From this experience I can suggest using a random vector of Euler rotation. This is much more convenient, since it is not necessary to calculate random trigonometric functions from random angles, as is the case when three Euler angles or four Rodrig-Hamilnon parameters are used (in quaternionic or matrix formalism).

Best regards.

The question boils down to: How to generate random points on a sphere (for which you will easily find solutions). Then you have to transform the vectors to Euler Angles e.g. using

https://sourceforge.net/projects/orilib/.

You may want to use MTEX http://mtex-toolbox.github.io/ . There you can generate, e.g. 1000 random rotations. with a single command

rot = rotation.rand(1000)

As a bonus you can generate random samples with respect to any distribution in orientation space.