This way is easiest...just get your arrays in, I have created some for demo purposes
>>> a = np.zeros(9,dtype=int).reshape((3,3))
>>> b = np.ones_like(a) # all ones/True
>>> c = b*2
>>> d = b*3
>>> stack = np.array([a,b,c,d])
>>> stack
array([[[0, 0, 0],
[0, 0, 0],
[0, 0, 0]],
[[1, 1, 1],
[1, 1, 1],
[1, 1, 1]],
[[2, 2, 2],
[2, 2, 2],
[2, 2, 2]],
[[3, 3, 3],
[3, 3, 3],
[3, 3, 3]]])
>>> per = 20
>>> np.percentile(stack,per,axis=0)
array([[ 0.6, 0.6, 0.6],
[ 0.6, 0.6, 0.6],
[ 0.6, 0.6, 0.6]])
>>> per = 50
>>> np.percentile(stack,per,axis=0)
array([[ 1.5, 1.5, 1.5],
[ 1.5, 1.5, 1.5],
[ 1.5, 1.5, 1.5]])
>>>
You will only be limited by memory as to whether you can read all the rasters in as arrays all at once.