In [1]:
from numpy import array, arange, argmin, sum, mean, var, size, zeros,\
where, histogram
from numpy.random import normal
from matplotlib.pyplot import figure, plot, hist, bar, xlabel, ylabel,\
title, show, savefig
#x = normal(0, 100, 1e2) # Generate n pseudo-random numbers whit(mu,sigma,n)
x = [4.37,3.87,4.00,4.03,3.50,4.08,2.25,4.70,1.73,4.93,1.73,4.62,\
3.43,4.25,1.68,3.92,3.68,3.10,4.03,1.77,4.08,1.75,3.20,1.85,\
4.62,1.97,4.50,3.92,4.35,2.33,3.83,1.88,4.60,1.80,4.73,1.77,\
4.57,1.85,3.52,4.00,3.70,3.72,4.25,3.58,3.80,3.77,3.75,2.50,\
4.50,4.10,3.70,3.80,3.43,4.00,2.27,4.40,4.05,4.25,3.33,2.00,\
4.33,2.93,4.58,1.90,3.58,3.73,3.73,1.82,4.63,3.50,4.00,3.67,\
1.67,4.60,1.67,4.00,1.80,4.42,1.90,4.63,2.93,3.50,1.97,4.28,\
1.83,4.13,1.83,4.65,4.20,3.93,4.33,1.83,4.53,2.03,4.18,4.43,\
4.07,4.13,3.95,4.10,2.27,4.58,1.90,4.50,1.95,4.83,4.12]
x_max = max(x)
x_min = min(x)
N_MIN = 4 #Minimum number of bins (integer)
N_MAX = 50 #Maximum number of bins (integer)
N = arange(N_MIN, N_MAX) # #of Bins
D = (x_max - x_min) / N #Bin size vector
C = zeros(shape=(size(D), 1))
C = zeros(size(D))
#Computation of the cost function
for i in xrange(size(N)):
ki = histogram(x, bins=N[i])
ki = ki[0]
k = mean(ki) #Mean of event count
v = var(ki) #Variance of event count
C[i] = (2 * k - v) / (D[i]**2) #The cost Function
#Optimal Bin Size Selection
idx = argmin(C)
cmin = C[idx]
optD = D[idx]
print '#', cmin, N[idx], optD
fig = figure()
if 0: # Two ways of plotting histogram.
# http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.hist
hist(x, bins=N[idx])
else:
# http://docs.scipy.org/doc/numpy/reference/generated/numpy.histogram.html
counts, bins = histogram(x, bins=N[idx])
bar(bins[:-1], counts, width=optD)
title("Histogram")
ylabel("Frequency")
xlabel("Value")
show()#savefig('Hist.png')
fig = figure()
plot(D, C, '.b', optD, cmin, '*r')
show()#savefig('Fobj.png')
In [2]:
from numpy import array, arange, argmin, sum, mean, var, size, zeros, where, histogram
from numpy.random import normal
from matplotlib.pyplot import figure, plot, hist, bar, xlabel, ylabel,title, show, savefig
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
x = normal(0, 100, 1e4) # Generate n pseudo-random numbers whit(mu,sigma,n)
y = normal(0, 100, 1e4) # Generate n pseudo-random numbers whit(mu,sigma,n)
x_max = max(x)
x_min = min(x)
y_max = max(y)
y_min = min(y)
Nx_MIN = 1 #Minimum number of bins in x (integer)
Nx_MAX = 100 #Maximum number of bins in x (integer)
Ny_MIN = 1 #Minimum number of bins in y (integer)
Ny_MAX = 100 #Maximum number of bins in y (integer)
Nx = arange(Nx_MIN, Nx_MAX) # #of Bins
Ny = arange(Ny_MIN, Ny_MAX) # #of Bins
Dx = (x_max - x_min) / Nx #Bin size vector
Dy = (y_max - y_min) / Ny #Bin size vector
Dxy=[]
for i in Dx: #Bin size vector
a=[]
for j in Dy: #Bin size vector
a.append((i,j))
Dxy.append(a)
Dxy=array( Dxy, dtype=[('x', float),('y', float)]) #matrix of bin size vector
Cxy=zeros(np.shape(Dxy))
Cxy__Dxy_plot=[] #save data to plot in scatterplot x,y,z
#Computation of the cost function to x and y
for i in xrange(size(Nx)):
for j in xrange(size(Ny)):
ki = np.histogram2d(x,y, bins=(Nx[i],Ny[j]))
ki = ki[0] #The mean and the variance are simply computed from the event counts in all the bins of the 2-dimensional histogram.
k = mean(ki) #Mean of event count
v = var(ki) #Variance of event count
Cxy[i,j] = (2 * k - v) / ( (Dxy[i,j][0]*Dxy[i,j][1])**2 ) #The cost Function
#(Cxy , Dx , Dy)
Cxy__Dxy_plot.append((Cxy[i,j] , Dxy[i,j][0] , Dxy[i,j][1]))#Save result of cost function to scatterplot
Cxy__Dxy_plot = np.array( Cxy__Dxy_plot , dtype=[('Cxy', float),('Dx', float), ('Dy', float)]) #Save result of cost function to scatterplot
#Optimal Bin Size Selection
#combination of i and j that produces the minimum cost function
idx_min_Cxy=np.where(Cxy == np.min(Cxy)) #get the index of the min Cxy
Cxymin=Cxy[idx_min_Cxy[0][0],idx_min_Cxy[1][0]] #value of the min Cxy
print sum(Cxy==Cxymin) #check if there is only one min value
optDxy=Dxy[idx_min_Cxy[0][0],idx_min_Cxy[1][0]]#get the bins size pairs that produces the minimum cost function
optDx=optDxy[0]
optDy=optDxy[1]
idx_Nx=idx_min_Cxy[0][0]#get the index in x that produces the minimum cost function
idx_Ny=idx_min_Cxy[1][0]#get the index in y that produces the minimum cost function
print '#', Cxymin, Nx[idx_Nx], optDx
print '#', Cxymin, Ny[idx_Ny], optDy
#PLOTS
#plot histogram2d
fig = figure()
H, xedges, yedges = np.histogram2d(x, y,bins=[Nx[idx_Nx],Ny[idx_Ny]])
Hmasked = np.ma.masked_where(H==0,H)
plt.imshow( Hmasked.T,extent = [xedges[0], xedges[-1], yedges[0], yedges[-1] ] ,interpolation='nearest',origin='lower',aspect='auto',cmap=plt.cm.Spectral)
plt.ylabel("y")
plt.xlabel("x")
plt.colorbar().set_label('z')
plt.show()
#plot scatterplot3d to Dx,Dy and Cxy
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
x=Cxy__Dxy_plot['Dx']
y=Cxy__Dxy_plot['Dy']
z =Cxy__Dxy_plot['Cxy']
ax.scatter(x, y, z, c=z, marker='o')
ax.set_xlabel('Dx')
ax.set_ylabel('Dy')
ax.set_zlabel('Cxy')
plt.draw()
ax.scatter( [optDx], [optDy],[Cxymin], marker='v', s=150,c="red")
ax.text(optDx, optDy,Cxymin, "Cxy min", color='red')
plt.draw()
plt.show()
In [4]:
# Author: Jake VanderPlas
# License: BSD
# The figure produced by this code is published in the textbook
# "Statistics, Data Mining, and Machine Learning in Astronomy" (2013)
# For more information, see http://astroML.github.com
# To report a bug or issue, use the following forum:
# https://groups.google.com/forum/#!forum/astroml-general
import numpy as np
from matplotlib import pyplot as plt
from scipy import stats
from astroML.plotting import hist
#----------------------------------------------------------------------
# This function adjusts matplotlib settings for a uniform feel in the textbook.
# Note that with usetex=True, fonts are rendered with LaTeX. This may
# result in an error if LaTeX is not installed on your system. In that case,
# you can set usetex to False.
from astroML.plotting import setup_text_plots
setup_text_plots(fontsize=8, usetex=True)
def plot_labeled_histogram(style, data, name,
x, pdf_true, ax=None,
hide_x=False,
hide_y=False):
if ax is not None:
ax = plt.axes(ax)
counts, bins, patches = hist(data, bins=style, ax=ax,
color='k', histtype='step', normed=True)
ax.text(0.95, 0.93, '%s:\n%i bins' % (name, len(counts)),
transform=ax.transAxes,
ha='right', va='top')
ax.fill(x, pdf_true, '-', color='#CCCCCC', zorder=0)
if hide_x:
ax.xaxis.set_major_formatter(plt.NullFormatter())
if hide_y:
ax.yaxis.set_major_formatter(plt.NullFormatter())
ax.set_xlim(-5, 5)
return ax
#------------------------------------------------------------
# Set up distributions:
Npts = 5000
np.random.seed(0)
x = np.linspace(-6, 6, 1000)
# Gaussian distribution
data_G = stats.norm(0, 1).rvs(Npts)
pdf_G = stats.norm(0, 1).pdf(x)
# Non-Gaussian distribution
distributions = [stats.laplace(0, 0.4),
stats.norm(-4.0, 0.2),
stats.norm(4.0, 0.2)]
weights = np.array([0.8, 0.1, 0.1])
weights /= weights.sum()
data_NG = np.hstack(d.rvs(int(w * Npts))
for (d, w) in zip(distributions, weights))
pdf_NG = sum(w * d.pdf(x)
for (d, w) in zip(distributions, weights))
#------------------------------------------------------------
# Plot results
fig = plt.figure(figsize=(5, 2.5))
fig.subplots_adjust(hspace=0, left=0.07, right=0.95, wspace=0.05, bottom=0.15)
ax = [fig.add_subplot(3, 2, i + 1) for i in range(6)]
# first column: Gaussian distribution
plot_labeled_histogram('scotts', data_G, 'Scott\'s Rule', x, pdf_G,
ax=ax[0], hide_x=True, hide_y=True)
plot_labeled_histogram('freedman', data_G, 'Freed.-Diac.', x, pdf_G,
ax=ax[2], hide_x=True, hide_y=True)
plot_labeled_histogram('knuth', data_G, 'Knuth\'s Rule', x, pdf_G,
ax=ax[4], hide_x=False, hide_y=True)
ax[0].set_title('Gaussian distribution')
ax[2].set_ylabel('$p(x)$')
ax[4].set_xlabel('$x$')
# second column: non-gaussian distribution
plot_labeled_histogram('scotts', data_NG, 'Scott\'s Rule', x, pdf_NG,
ax=ax[1], hide_x=True, hide_y=True)
plot_labeled_histogram('freedman', data_NG, 'Freed.-Diac.', x, pdf_NG,
ax=ax[3], hide_x=True, hide_y=True)
plot_labeled_histogram('knuth', data_NG, 'Knuth\'s Rule', x, pdf_NG,
ax=ax[5], hide_x=False, hide_y=True)
ax[1].set_title('non-Gaussian distribution')
ax[5].set_xlabel('$x$')
plt.show()
In [9]:
# Define our test distribution: a mix of Cauchy-distributed variables
import numpy as np
from scipy import stats
import pylab as pl
%pylab inline
np.random.seed(0)
x = np.concatenate([stats.cauchy(-5, 1.8).rvs(500),
stats.cauchy(-4, 0.8).rvs(2000),
stats.cauchy(-1, 0.3).rvs(500),
stats.cauchy(2, 0.8).rvs(1000),
stats.cauchy(4, 1.5).rvs(500)])
# truncate values to a reasonable range
x = x[(x > -15) & (x < 15)]
plt.figure()
pl.hist(x, normed=True)
plt.figure()
pl.hist(x, bins=100, normed=True)
Out[9]:
In [10]:
def bayesian_blocks(t):
"""Bayesian Blocks Implementation
By Jake Vanderplas. License: BSD
Based on algorithm outlined in http://adsabs.harvard.edu/abs/2012arXiv1207.5578S
Parameters
----------
t : ndarray, length N
data to be histogrammed
Returns
-------
bins : ndarray
array containing the (N+1) bin edges
Notes
-----
This is an incomplete implementation: it may fail for some
datasets. Alternate fitness functions and prior forms can
be found in the paper listed above.
"""
# copy and sort the array
t = np.sort(t)
N = t.size
# create length-(N + 1) array of cell edges
edges = np.concatenate([t[:1],
0.5 * (t[1:] + t[:-1]),
t[-1:]])
block_length = t[-1] - edges
# arrays needed for the iteration
nn_vec = np.ones(N)
best = np.zeros(N, dtype=float)
last = np.zeros(N, dtype=int)
#-----------------------------------------------------------------
# Start with first data cell; add one cell at each iteration
#-----------------------------------------------------------------
for K in range(N):
# Compute the width and count of the final bin for all possible
# locations of the K^th changepoint
width = block_length[:K + 1] - block_length[K + 1]
count_vec = np.cumsum(nn_vec[:K + 1][::-1])[::-1]
# evaluate fitness function for these possibilities
fit_vec = count_vec * (np.log(count_vec) - np.log(width))
fit_vec -= 4 # 4 comes from the prior on the number of changepoints
fit_vec[1:] += best[:K]
# find the max of the fitness: this is the K^th changepoint
i_max = np.argmax(fit_vec)
last[K] = i_max
best[K] = fit_vec[i_max]
#-----------------------------------------------------------------
# Recover changepoints by iteratively peeling off the last block
#-----------------------------------------------------------------
change_points = np.zeros(N, dtype=int)
i_cp = N
ind = N
while True:
i_cp -= 1
change_points[i_cp] = ind
if ind == 0:
break
ind = last[ind - 1]
change_points = change_points[i_cp:]
return edges[change_points]
In [11]:
# plot a standard histogram in the background, with alpha transparency
H1 = hist(x, bins=200, histtype='stepfilled',
alpha=0.2, normed=True)
# plot an adaptive-width histogram on top
H2 = hist(x, bins=bayesian_blocks(x), color='black',
histtype='step', normed=True)
In [ ]:
hmqdata_sorted_full.csv
In [ ]:
In [13]:
from numpy import array, arange, argmin, sum, mean, var, size, zeros,\
where, histogram
from numpy.random import normal
from matplotlib.pyplot import figure, plot, hist, bar, xlabel, ylabel,\
title, show, savefig
%matplotlib inline
import astropy.io.ascii as ascii
#x = normal(0, 100, 1e2) # Generate n pseudo-random numbers whit(mu,sigma,n)
#x = [4.37,3.87,4.00,4.03,3.50,4.08,2.25,4.70,1.73,4.93,1.73,4.62,\
#3.43,4.25,1.68,3.92,3.68,3.10,4.03,1.77,4.08,1.75,3.20,1.85,\
#4.62,1.97,4.50,3.92,4.35,2.33,3.83,1.88,4.60,1.80,4.73,1.77,\
#4.57,1.85,3.52,4.00,3.70,3.72,4.25,3.58,3.80,3.77,3.75,2.50,\
#4.50,4.10,3.70,3.80,3.43,4.00,2.27,4.40,4.05,4.25,3.33,2.00,\
#4.33,2.93,4.58,1.90,3.58,3.73,3.73,1.82,4.63,3.50,4.00,3.67,\
#1.67,4.60,1.67,4.00,1.80,4.42,1.90,4.63,2.93,3.50,1.97,4.28,\
#1.83,4.13,1.83,4.65,4.20,3.93,4.33,1.83,4.53,2.03,4.18,4.43,\
#4.07,4.13,3.95,4.10,2.27,4.58,1.90,4.50,1.95,4.83,4.12]
hmqdat=ascii.read("./hmqdata_sorted_full.csv")
x=hmqdat['Z']
x_max = max(x)
x_min = min(x)
N_MIN = 4 #Minimum number of bins (integer)
N_MAX = 50 #Maximum number of bins (integer)
N = arange(N_MIN, N_MAX) # #of Bins
D = (x_max - x_min) / N #Bin size vector
C = zeros(shape=(size(D), 1))
C = zeros(size(D))
#Computation of the cost function
for i in xrange(size(N)):
ki = histogram(x, bins=N[i])
ki = ki[0]
k = mean(ki) #Mean of event count
v = var(ki) #Variance of event count
C[i] = (2 * k - v) / (D[i]**2) #The cost Function
#Optimal Bin Size Selection
idx = argmin(C)
cmin = C[idx]
optD = D[idx]
print '#', cmin, N[idx], optD
fig = figure()
if 0: # Two ways of plotting histogram.
# http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.hist
hist(x, bins=N[idx])
else:
# http://docs.scipy.org/doc/numpy/reference/generated/numpy.histogram.html
counts, bins = histogram(x, bins=N[idx])
bar(bins[:-1], counts, width=optD)
title("Histogram")
ylabel("Frequency")
xlabel("Value")
show()#savefig('Hist.png')
fig = figure()
plot(D, C, '.b', optD, cmin, '*r')
show()#savefig('Fobj.png')
In [14]:
from numpy import array, arange, argmin, sum, mean, var, size, zeros,\
where, histogram
from numpy.random import normal
from matplotlib.pyplot import figure, plot, hist, bar, xlabel, ylabel,\
title, show, savefig
%matplotlib inline
import astropy.io.ascii as ascii
hmqdat=ascii.read("./hmqdata_sorted_full.csv")
x=hmqdat['Z']
x_max = max(x)
x_min = min(x)
N_MIN = 4 #Minimum number of bins (integer)
N_MAX = 15 #Maximum number of bins (integer)
N = arange(N_MIN, N_MAX) # #of Bins
D = (x_max - x_min) / N #Bin size vector
C = zeros(shape=(size(D), 1))
C = zeros(size(D))
#Computation of the cost function
for i in xrange(size(N)):
ki = histogram(x, bins=N[i])
ki = ki[0]
k = mean(ki) #Mean of event count
v = var(ki) #Variance of event count
C[i] = (2 * k - v) / (D[i]**2) #The cost Function
#Optimal Bin Size Selection
idx = argmin(C)
cmin = C[idx]
optD = D[idx]
print '#', cmin, N[idx], optD
fig = figure()
if 0: # Two ways of plotting histogram.
# http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.hist
hist(x, bins=N[idx])
else:
# http://docs.scipy.org/doc/numpy/reference/generated/numpy.histogram.html
counts, bins = histogram(x, bins=N[idx])
bar(bins[:-1], counts, width=optD)
title("Histogram")
ylabel("Frequency")
xlabel("Value")
show()#savefig('Hist.png')
fig = figure()
plot(D, C, '.b', optD, cmin, '*r')
show()#savefig('Fobj.png')
In [15]:
len(x)
Out[15]:
In [ ]: