In [1]:
import random
import numpy as np
import matplotlib.pylab as plt
import matplotlib.lines as mlines
In [6]:
## sample and plot stochastic block model
np.random.seed(123)
plt.figure(1)
n=100
B=np.array([[0.5, 0.2], [0.2, 0.5]])
pi=0.5
Z=(np.random.random((n,))>pi)*1
# print Z
A=np.zeros((n,n))
for i in range(n):
for j in range(n):
A[i,j]=np.random.random()<B[Z[i],Z[j]]
A=A*1
I=np.argsort(Z)
# print I
# B=deg[I]
plt.subplot(121)
plt.spy(A, marker='o', markersize=0.5)
plt.subplot(122)
plt.spy(A[np.ix_(I,I)], marker='o', markersize=0.5)
plt.show()
## estimate B & pi
pihat=np.sum(1.*Z)/n
n1=np.sum(Z==1)
# print n1
# print (A[np.ix_(Z==0,Z==0)])
Bhat=np.zeros((2,2))
for i in range(2):
for j in range(2):
Bhat[i,j]=np.sum(A[np.ix_(Z==i,Z==j)])/n1**2.
print pihat
print Bhat
print pi
print B
In [ ]: