bboxes.shape => (32,125)
x_cell.shape => (32,)
y_cell.shape => (32,)
output.shape => (13,13,32,125) 
     
for i in range(32):
    output[x_cell[i], y_cell[i], i, :] = bboxes[i]

#eg shapes of indexing: output[(1,), (1,), (1,), :] = (125,) 
     
output[x_cell,y_cell,np.arange(32)[:,None],:] = bboxes[:,None,:] 
     
output[x_cell,y_cell,np.arange(32)[:,None]] = bboxes[:,None] 
     
In [53]: # Setup
    ...: np.random.seed(0)
    ...: bboxes = np.random.rand(32,125)
    ...: x_cell = np.random.randint(0,13,(32,13))
    ...: y_cell = np.random.randint(0,13,(32,13))
    ...: 
    ...: output = np.zeros((13,13,32,125))
    ...: for i in range(32):
    ...:     output[x_cell[i], y_cell[i], i, :] = bboxes[i]
    ...: 
    ...: out = np.zeros((13,13,32,125))
    ...: out[x_cell,y_cell,np.arange(32)[:,None]] = bboxes[:,None]
    ...: 
    ...: print np.allclose(out, output)
True