# Show with $\chi^2$ and/or QQ plots that data is Poisson distributed (+ realization in Python)

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)

• 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. Nov 14, 2016 at 8:02