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