Why does Phase-locking value show spurious effects at the beginning & end of a signal?
What was previously posted as an answer should have been an extended comment -- sorry for that. There simply was not enough space for what I thought I'd write about how the Octave code seemed to behave. Now that the code is cleared, what should have been the answer in the first place (given the rather obvious signs), is now the answer: spectral leakage.
The Hilbert transform is calculated by zeroing half the FFT of the signal and then recomputing the analytical signal through the IFFT. Brute-force zeroing will always lead to some leakage due to the nature of the limited length of the data. This can be easily tested with your code. s1 is a simple sin(), so its Hilbert transform will be a -cos(). Simply plot the difference between -cos()-hx.imag, and you'll see the leakage. The same goes for the second signal.
This is the changed code in Python:
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
f = 10
t = np.linspace(0, 100/f, 1000)
N = 50f = 10
nt = 1000
t = np.linspace(0, 100/f, nt)
N = 50
hor = np.ones(nt)
ver = np.ones(N)
r = np.array( [ np.random.uniform(low=0, high=2*np.pi) for i in range(N) ] )
s1 = np.array( [ np.sin( 2*np.pi*f*t ) for i in range(N) ] )
s2 = np.array( np.sin( 2*np.pi*( f*ver[:,None]*t + hor*r[:,None] ) ) )
c1 = -np.array( [ np.cos(2*np.pi*f*t) for i in range(N) ] )
c2 = -np.array( np.cos( 2*np.pi*( f*ver[:,None]*t + hor*r[:,None] ) ) )
hx = signal.hilbert(s1)
hy = signal.hilbert(s2)
nx = hx/np.absolute(hx)
ny = hy/np.absolute(hy)
p = np.absolute( np.mean( nx*ny.conj(), axis=0 ) )
plt.plot(p)
plt.show()

and this is the test (c2 will show non-overlapping plots due to the randomness):
plt.plot(c1.transpose() - hx.imag.transpose())
plt.show()


The Octave code, for checking:
f = 10
t = linspace(0, 100/f, 1000);
N = 50;
c = ones(N, 1);
ft = f*t.*c;
s1 = sin(2*pi*ft);
c1 = -cos(2*pi*ft);
r = rand(N,1).*c;
s2 = sin(2*pi*(ft + r));
c2 = -cos(2*pi*(ft + r));
hx = hilbert(s1');
hy = hilbert(s2');
nx = hx./abs(hx);
ny = hy./abs(hy);
p = abs(mean((nx.*conj(ny))'));
grid on
subplot(2, 1, 1)
plot(p)
subplot(2, 1, 2)
plot(c1' - imag(hx))