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


0.51
[[ 0.46482122  0.19992311]
 [ 0.1987697   0.49480969]]
0.5
[[ 0.5  0.2]
 [ 0.2  0.5]]

In [ ]: