import numpy import scipy.io def z_measure(z, pp): #Computes the probability of a partition pp under the z-measure, but does not divide by (n! * rising factorial). dim_squared = pp.dimension()^2 content_prod = prod([z + c[1] - c[0] for c in pp.cells()])^2 return dim_squared * content_prod def return_lambda1_dist(z = 0.5, n = 30): #Returns a vector of tuples (i,p), where: #i is of the form a/n, and p is the probability that lambda_1 = a, under the z-measure on Y_n. l_probs = [0 for i in xrange(n+1)] for par in Partitions(n): l_probs[par[0]] += z_measure(z,par) n_factorial = factorial(n) pochhammer = rising_factorial(z^2, n) v = [(i * 1.0 / len(l_probs), l_probs[i] / (n_factorial * pochhammer) * 1.0) for i in xrange(len(l_probs))] return v def F_graphs(n = 60): points_list = [] points_deriv = [] for z in [0.25, 0.5, 0.75, 1.5, 2.0, 2.5, 3.0, 3.5]: #The probability distribution of lambda_1/n under the n-th z-measure approximates the #probability distribution of alpha_1 under the z-measure. lambda1_dist = return_lambda1_dist(z, n=n) #Computing commulative distribution function of lambda_1 from the probability distribution function: lambda1_dist_cumsum = numpy.cumsum([t[1] for t in lambda1_dist]).tolist() lambda1_dist_CDF = [(lambda1_dist[i][0], lambda1_dist_cumsum[i]) for i in xrange(len(lambda1_dist))] points_list.append((z, lambda1_dist_CDF)) #Computing discrete derivative of commulative distribution function: points_deriv.append((z,[(lambda1_dist_CDF[i][0], lambda1_dist_CDF[i+1][1]-lambda1_dist_CDF[i][1]) for i in xrange(len(lambda1_dist)-1)])) #Graph of (approximation of) F_z for z = 0.25, z=0.5, z=0.75: g = Graphics() g += plot(points(points_list[0][1], rgbcolor=(0.0,0.0,0.5), pointsize=15, legend_label = 'z=0.25', legend_color = 'black')) g += plot(points(points_list[1][1], rgbcolor=(0.0,0.5,0.8), pointsize=15, legend_label = 'z=0.5', legend_color = 'black')) g += plot(points(points_list[2][1], rgbcolor=(0.0,0.6,0.0), pointsize=15, legend_label = 'z=0.75', legend_color ='black')) show(g) #Graph of (approximation of) F_z for z = 1.5, 2.0, 2.5, 3.0, 3.5: g = Graphics() g += plot(points(points_list[3][1], rgbcolor = (0.0,0.0,0.5), pointsize = 15, legend_label = 'z=1.5', legend_color = 'black')) g += plot(points(points_list[4][1], rgbcolor = (0.0,0.5,0.8), pointsize = 15, legend_label = 'z=2.0', legend_color = 'black')) g += plot(points(points_list[5][1], rgbcolor = (0.0,0.6,0.0), pointsize = 15, legend_label = 'z=2.5', legend_color = 'black')) g += plot(points(points_list[6][1], rgbcolor = (0.6,0.5,0.0), pointsize = 15, legend_label = 'z=3.0', legend_color = 'black')) g += plot(points(points_list[7][1], rgbcolor = (0.8,0.2,0.0), pointsize = 15, legend_label = 'z=3.5', legend_color = 'black')) show(g) #Graph of discrete derivative of (approximation of) F_z for z = 1.5, 2.0, 2.5, 3.0, 3.5: g = Graphics() g += plot(points(points_deriv[3][1][:-1], rgbcolor = (0.0,0.0,0.5), pointsize = 15, legend_label = 'z=1.5', legend_color = 'black')) g += plot(points(points_deriv[4][1][:-1], rgbcolor = (0.0,0.5,0.8), pointsize = 15, legend_label = 'z=2.0', legend_color = 'black')) g += plot(points(points_deriv[5][1][:-1], rgbcolor = (0.0,0.6,0.0), pointsize = 15, legend_label = 'z=2.5', legend_color = 'black')) g += plot(points(points_deriv[6][1][:-1], rgbcolor = (0.6,0.5,0.0), pointsize = 15, legend_label = 'z=3.0', legend_color = 'black')) g += plot(points(points_deriv[7][1][:-1], rgbcolor = (0.8,0.2,0.0), pointsize = 15, legend_label = 'z=3.5', legend_color = 'black')) show(g) def return_G_dist(z = 0.5, n = 35): #Returns a vector of tuples (a/n ,T(n;a)), where: #a ranges from 1 to n, and T(n;a) is defined in Theorem 1.1. lls = [[1 for i in xrange(n+2)]] for i in xrange(1, n+1): ll = map(lambda x: x[1], return_lambda1_dist(z,i)) ll = numpy.cumsum(ll).tolist() lls.append(ll + [ll[-1] for i in xrange(n+2-len(ll))]) #the list lls satisfies lls[i][a] = P_{Y_i}(\lambda_1 <= a), also for a>i. probs = [0] + [0 for i in xrange(n)] for i in xrange(n+1): j= n-i bin1 = falling_factorial(i+z^2-1,i)/factorial(i) bin2 = falling_factorial(j+z^2-1,j)/factorial(j) for a in xrange(1,len(probs)): probs[a] += bin1*bin2*lls[i][a-1]*lls[j][a] return probs def G_graphs(n = 50): probs = return_G_dist(n=n) #We want to normalize T(n;a) by T(n;n). Asymptotically, T(n,n) ~ 1/sqrt(pi *n), #but in fact a simple closed form formula for T(n,n) exists: normalization_fact = binomial(n - 0.5, n) - binomial(n - 0.5, n)^2 #The ratio T(n;a)/T(n;n) is important since it approximates G(a/n): G_approx = [(i * 1.0 / len(probs), probs[i] * 1.0 / normalization_fact) for i in xrange(len(probs))] #Graph of (approximation of) G(s): show(points(G_approx, rgbcolor = (0.0,0.0,0.5), pointsize = 15)) #Graph of (approximation of) derivative of G(s). Defined as the discrete derivatve of T(n;a)/T(n;n): show(points([(G_approx[i][0], G_approx[i+1][1]-G_approx[i][1]) for i in xrange(len(G_approx)-1)], rgbcolor = (0.0,0.0,0.5), pointsize = 15)) def variance_graph(path = 'results.mat', n = 50): probs = return_G_dist(n=n) normalization_fact = binomial(n - 0.5, n) - binomial(n - 0.5, n)^2 G_approx = [(i * 1.0 / len(probs), probs[i] * 1.0 / normalization_fact) for i in xrange(len(probs))] mat = scipy.io.loadmat(path) x_axis = mat['x_axis'] qcount = mat['qcount'] K = 0.76422365 #The Landau-Ramanujan constant G_timesK = [(1-G_approx[i][0], G_approx[i][1] * K) for i in xrange(len(G_approx))] Empirical = [(1-x_axis[i], qcount[i]) for i in xrange(len(x_axis))] g = Graphics() g += plot(spline(G_timesK), 0, 2, rgbcolor = (0.0,0.0,0.5), legend_label = 'Prediction', legend_color = 'black') g += plot(points(Empirical, rgbcolor = (0.0,0.5,0.0), legend_label = 'Data', legend_color = 'black', pointsize = 3)) show(g)