2
$\begingroup$

I want to check that the following limiting behaviour for Erdos-Renyi graph with parameter $ER(\lambda/n)$ for $n \to \infty$ :

$$\Delta_n \to Poi \left(\frac{\lambda^3}{6}\right)$$

agrees, where $\Delta_n$ - amount of triangles in the graph.

I made few simulations and stuck with 2 problems:

  1. I do not know how to create QQ plot in python for the discrete( Poi) distribution.( I am new to statistic at all)
  2. I get the terrible $\chi^2$ results.

I have included my trials and code below:

lambda: 5.0,# of experiments: 100, # of nodes: 100
chisquare_value: 27.2365689846, degree: 40

lambda: 5.0,# of experiments: 100, # of nodes: 100
chisquare_value: 27.7760822214, degree: 33

Which seems ok if I compare with Upper-tail critical values from: http://www.itl.nist.gov/div898/handbook/eda/section3/eda3674.htm

where for 33 and 40 degrees of freedom we have: 43.745 and 51.805 for probability $= 0.9$. so, If I tracked the following correct, I go further and increase the number of nodes at the graph ($n$) ( which in theory should increase the similarity of distribution of triangles in a graph and $Poi$ distribution.

lambda: 5.0,# of experiments: 100, # of nodes: 1000
chisquare_value: 36.521541423, degree: 32


lambda: 5.0,# of experiments: 100, # of nodes: 1000
chisquare_value: 34.991726686, degree: 31

I get still reasonable results (again, if I treat it accordingly) as the $\chi^2$ values for the probability $=0.9$ are 42.585 and 41.422 respectively. And for $10000$ nodes it is ok.

lambda: 5.0,# of experiments: 100, # of nodes: 10000
chisquare_value: 38.1744412929, degree: 36


lambda: 5.0,# of experiments: 100, # of nodes: 10000
chisquare_value: 34.1908530453, degree: 28

But! As soon as I increase number of trials (from 100 to 1000) (for $n=1000$) I get:

lambda: 5.0,# of experiments: 1000, # of nodes: 1000
chisquare_value: 4273.20483365, degree: 34
lambda: 5.0,# of experiments: 1000, # of nodes: 1000
chisquare_value: 4252.08983769, degree: 37

and built in $\chi^2$ does not work good for any of those tests:

print scipy.stats.chisquare(test_observed_freq, test_expected_freq, ddof = 1)

my code:

import networkx, numpy, scipy.stats, math, pylab
import matplotlib.pyplot as plt



#initial setup
numb_of_nodes_per_exp = 1000
numb_of_exps = 100
list_of_exps = [numb_of_nodes_per_exp] * numb_of_exps
lambda_parameter = 3.0 # may be different
print "lambda: %s,# of experiments: %s, # of nodes: %s" \
%(lambda_parameter ,numb_of_exps,numb_of_nodes_per_exp)

def generate_tri():
    triangles_per_exp = []
    for n in range(numb_of_exps):
        G = networkx.fast_gnp_random_graph(numb_of_nodes_per_exp,lambda_parameter / numb_of_nodes_per_exp)
        triangles_per_node = networkx.triangles(G).values()
        triangles_per_exp.append(  sum( triangles_per_node ) / 3   )
    return triangles_per_exp

def count_freq_of_tri_apearance(triangles_per_exp):
    list_of_tri_and_freq = []    # amount of tri and freq 
    for k in triangles_per_exp:
        frequency_of_k =  triangles_per_exp.count(k)
        if frequency_of_k != 0: 
            list_of_tri_and_freq.append([k,frequency_of_k])
    return  list_of_tri_and_freq

def expected_freq(tri_and_their_freq):
    probability_list = []
    expected_freq_list = []
    lambda_paramet_star = lambda_parameter**3 / 6

    for each_tri in range(len(tri_and_their_freq)):
        probability_list.append(   lambda_paramet_star ** tri_and_their_freq[each_tri][0] \
        * math.exp(- lambda_paramet_star) / math.factorial(tri_and_their_freq[each_tri][0])    )

    number_of_tri_in_all_exp = sum([x[1] for x in tri_and_their_freq ] )

    for each_tri in range(len(tri_and_their_freq)):
        expected_freq_list.append(   probability_list[each_tri] *  number_of_tri_in_all_exp   )
    return expected_freq_list

def observed_freq(tri_and_their_freq):
    observed_freq = []
    for k in range(len(tri_and_their_freq)):
        observed_freq.append( tri_and_their_freq[k][1] )
    return  observed_freq

def chi_sq(exp_freq,obs_freq):
    a = 0.0
    b = 0.0
    max_bin = max(triangles_per_exp)
    p = 1
    for n in range(len(exp_freq)):
        a = a + ( obs_freq[n] - exp_freq[n])**2 
        b = b + exp_freq[n]
    chisquare  = a/b
    return [chisquare, max_bin - p - 1]


triangles_per_exp =  generate_tri()
#print triangles_per_exp

tri_and_their_freq =  count_freq_of_tri_apearance( triangles_per_exp )
#print tri_and_their_freq

test_expected_freq = expected_freq(tri_and_their_freq)
# print test_expected_freq

test_observed_freq = observed_freq(tri_and_their_freq)
# print test_observed_freq

chisquare_list = chi_sq(test_expected_freq,test_observed_freq)
chisquare_value = chisquare_list[0]
degree = chisquare_list[1]

print "chisquare_value: %s, degree: %s" %(chisquare_value, degree)
$\endgroup$
1
  • 1
    $\begingroup$ There's a "Poissonness plot" discussed in How to know if a data follows a Poisson distribution in R. Such things cannot show something is Poisson (since one can be arbitrarily close to a Poisson without being Poisson). However, they can sometimes make it clear when it isn't Poisson. It's not exactly a QQ plot, but it's somewhat similar. $\endgroup$
    – Glen_b
    Nov 14, 2016 at 8:02

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.