def plot_density_map(x, y, xbins, ybins, Nlevels=4, cbar=True, weights=None):
Z = np.histogram2d(x, y, bins=(xbins, ybins), weights=weights)[0].astype(float).T
# central values
lt = get_centers_from_bins(xbins)
lm = get_centers_from_bins(ybins)
cX, cY = np.meshgrid(lt, lm)
X, Y = np.meshgrid(xbins, ybins)
im = plt.pcolor(X, Y, Z, cmap=plt.cm.Blues)
plt.contour(cX, cY, Z, levels=nice_levels(Z, Nlevels), cmap=plt.cm.Greys_r)
if cbar:
cb = plt.colorbar(im)
else:
cb = None
plt.xlim(xbins[0], xbins[-1])
plt.ylim(ybins[0], ybins[-1])
try:
plt.tight_layout()
except Exception as e:
print(e)
return plt.gca(), cb
python类contour()的实例源码
def joint_density(X, Y, bounds=None):
"""
Plots joint distribution of variables.
Inherited from method in src/graphics.py module in project
git://github.com/aflaxman/pymc-example-tfr-hdi.git
"""
if bounds:
X_min, X_max, Y_min, Y_max = bounds
else:
X_min = X.min()
X_max = X.max()
Y_min = Y.min()
Y_max = Y.max()
pylab.plot(X, Y, linestyle='none', marker='o', color='green', mec='green', alpha=.2, zorder=-99)
gkde = scipy.stats.gaussian_kde([X, Y])
x,y = pylab.mgrid[X_min:X_max:(X_max-X_min)/25.,Y_min:Y_max:(Y_max-Y_min)/25.]
z = pylab.array(gkde.evaluate([x.flatten(), y.flatten()])).reshape(x.shape)
pylab.contour(x, y, z, linewidths=2)
pylab.axis([X_min, X_max, Y_min, Y_max])
def scatter_plot(x, y, ellipse=False, levels=[0.99, 0.95, 0.68], color='w', ax=None, **kwargs):
if ax is None:
ax = plt.gca()
if faststats is not None:
im, e = faststats.fastkde.fastkde(x, y, (50, 50), adjust=2.)
V = im.max() * np.asarray(levels)
plt.contour(im.T, levels=V, origin='lower', extent=e, linewidths=[1, 2, 3], colors=color)
ax.plot(x, y, 'b,', alpha=0.3, zorder=-1, rasterized=True)
if ellipse is True:
data = np.vstack([x, y])
mu = np.mean(data, axis=1)
cov = np.cov(data)
error_ellipse(mu, cov, ax=plt.gca(), edgecolor="g", ls="dashed", lw=4, zorder=2)
def plot_reg_2D_stoc(X,stoc_vector):
deter_vec = np.invert(stoc_vector)
dom_max = np.amax(X[stoc_vector,:]) + 1
A = np.zeros((dom_max,dom_max))
stoc_indexs = np.arange(0,X.shape[0],1)[stoc_vector].astype(int)
for i,deter_element in enumerate(deter_vec):
if deter_element == True:
A[X[int(stoc_indexs[0]),:].astype(int), X[int(stoc_indexs[1]),:].astype(int)] = X[i,:]
pl.figure(i)
#ax = fig.gca(projection='3d')
#surf = ax.plot_surface(X[int(stoc_indexs[0]),:].astype(int), X[int(stoc_indexs[1]),:].astype(int),X[i,:], rstride=1, cstride=1,
#cmap=cm.coolwarm,linewidth=0, antialiased=False)
pl.contour(A,X[i,:])
#ax.zaxis.set_major_locator(LinearLocator(10))
#ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
#fig.colorbar(surf, shrink=0.5, aspect=5)
pl.show()
def plot_reg_2D_stoc(X,stoc_vector):
deter_vec = np.invert(stoc_vector)
dom_max = np.amax(X[stoc_vector,:]) + 1
A = np.zeros((dom_max,dom_max))
stoc_indexs = np.arange(0,X.shape[0],1)[stoc_vector].astype(int)
for i,deter_element in enumerate(deter_vec):
if deter_element == True:
A[X[int(stoc_indexs[0]),:].astype(int), X[int(stoc_indexs[1]),:].astype(int)] = X[i,:]
pl.figure(i)
#ax = fig.gca(projection='3d')
#surf = ax.plot_surface(X[int(stoc_indexs[0]),:].astype(int), X[int(stoc_indexs[1]),:].astype(int),X[i,:], rstride=1, cstride=1,
#cmap=cm.coolwarm,linewidth=0, antialiased=False)
pl.contour(A,X[i,:])
#ax.zaxis.set_major_locator(LinearLocator(10))
#ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
#fig.colorbar(surf, shrink=0.5, aspect=5)
pl.show()
def plot_2D_contour(states,p,labels,inter=False):
import pylab as pl
from pyme.statistics import expectation as EXP
exp = EXP((states,p))
X = np.unique(states[0,:])
Y = np.unique(states[1,:])
X_len = len(X)
Y_len = len(Y)
Z = np.zeros((X.max()+1,Y.max()+1))
for i in range(len(p)):
Z[states[0,i],states[1,i]] = p[i]
Z = np.where(Z < 1e-8,0.0,Z)
pl.clf()
XX, YY = np.meshgrid(X,Y)
pl.contour(range(X.max()+1),range(Y.max()+1),Z.T)
pl.axhline(y=exp[1])
pl.axvline(x=exp[0])
pl.xlabel(labels[0])
pl.ylabel(labels[1])
if inter == True:
pl.draw()
else:
pl.show()
def contour_area(self,edges,resolution=0.1,**limits):
"""Return the area of each contour interval within the limits.
edges list containing the contour edges (1D array)
resolution resolution of the spline map used to calculate the area
TODO: better handling of limits
"""
lim = self.limits.copy()
lim.update(limits)
XX,YY = numpy.mgrid[lim['minNMP']:lim['maxNMP']:resolution,
lim['minLID']:lim['maxLID']:resolution]
A = self.F(XX,YY)
e = numpy.asarray(edges)
area = numpy.zeros(len(e)-1,dtype=int)
area[:] = [len(A[numpy.logical_and(lo <= A, A < hi)]) for lo,hi in zip(e[:-1],e[1:])]
return area * resolution**2
def contour(self, *args, **kwargs):
defaults = {'origin': 'lower', 'cmap': plt.cm.Greys,
'levels': np.sort(self.nice_levels())}
defaults.update(**kwargs)
ax = kwargs.pop('ax', plt.gca())
return ax.contour(self.im.T, *args, extent=self.e, **defaults)
def plot(self, contour={}, scatter={}, **kwargs):
# levels = np.linspace(self.im.min(), self.im.max(), 10)[1:]
levels = self.nice_levels()
c_defaults = {'origin': 'lower', 'cmap': plt.cm.Greys_r, 'levels':
levels}
c_defaults.update(**contour)
c = self.contourf(**c_defaults)
lvls = np.sort(c.levels)
s_defaults = {'c': '0.0', 'edgecolor':'None', 's':2}
s_defaults.update(**scatter)
self.scatter(lvl=[lvls], **s_defaults)
def autoplot(self,mode='observable',runtype=None,figdir=None, **plotargs):
"""Make the standard plot of the data.
autoplot(AngleprojectedObservable, mode=MODE,runtype='RUNTYPE',figdir=None,**plotargs)
plotargs (see AngleprojectedObservable.plot()):
title None: generate title for each individual plot
using the observables legend
string: use provided title for _all_ plots
"": no title
min_contour |
max_contour +- boundaries for the contour plot
cmap pylab colormap instance or None for the default [cm.jet]
alpha alpha for the filled contours [1.0]
overlay function that is called with no arguments, eg
overlay = lambda : coRMSD.DeltaRMSD.plot(with_colorbar=False,clf=False,with_filled=False,linewidth=3,colors='white')
"""
import pylab
qualifier = self.qualifier or "" # gives just trailing slash if no self.qualifier
figdir = figdir or os.path.join(config.basedir,'figs',self.observable_type,qualifier)
filebasename = self.filebasename(runtype)
was_interactive = pylab.matplotlib.is_interactive()
pylab.matplotlib.interactive(False)
self._auto_plots(mode,filebasename,figdir,plotargs)
pylab.matplotlib.interactive(was_interactive) # revert to previous state
return figdir
def plot_all(self,runtype=None,figdir=None, **plotargs):
"""Make the 6 standard plots of the observable data.
plot_all(runtype='RUNTYPE',figdir=None,**plotargs)
plotargs (see AngleprojectedObservable.plot()):
title None: generate title for each individual plot
using the observables legend
string: use provided title for _all_ plots
"": no title
min_contour |
max_contour +- boundaries for the contour plot
cmap pylab colormap instance or None for the default [cm.jet]
alpha alpha for the filled contours [1.0]
overlay function that is called with no arguments, eg
overlay = lambda : coRMSD.DeltaRMSD.plot(with_colorbar=False,clf=False,with_filled=False,linewidth=3,colors='white')
"""
kwargs = dict(runtype=runtype,figdir=figdir,**plotargs)
self.autoplot('observable',**kwargs)
self.autoplot('stdev',**kwargs)
figdir = self.autoplot('counts',**kwargs)
print "-- Figures in %r." % figdir
return figdir
def autoplot(self,mode='observable',runtype=None,figdir=None, **plotargs):
"""Make the standard plot of the data.
autoplot(AngleprojectedObservable, mode=MODE,runtype='RUNTYPE',figdir=None,**plotargs)
plotargs (see AngleprojectedObservable.plot()):
title None: generate title for each individual plot
using the observables legend
string: use provided title for _all_ plots
"": no title
min_contour |
max_contour +- boundaries for the contour plot
cmap pylab colormap instance or None for the default [cm.jet]
alpha alpha for the filled contours [1.0]
overlay function that is called with no arguments, eg
overlay = lambda : coRMSD.DeltaRMSD.plot(with_colorbar=False,clf=False,with_filled=False,linewidth=3,colors='white')
"""
import pylab
qualifier = self.qualifier or "" # gives just trailing slash if no self.qualifier
figdir = figdir or os.path.join(config.basedir,'figs',self.observable_type,qualifier)
filebasename = self.filebasename(runtype)
was_interactive = pylab.matplotlib.is_interactive()
pylab.matplotlib.interactive(False)
self._auto_plots(mode,filebasename,figdir,plotargs)
pylab.matplotlib.interactive(was_interactive) # revert to previous state
return figdir
def plot_contours(self):
"""filled contour plot with lines"""
try:
from pylab import contour,contourf
except ImportError:
print "pylab is required to plot"
return None
contour(self.x,self.y,self.histogram.transpose())
contourf(self.x,self.y,self.histogram.transpose())
def _test(npoints=1000):
"""test Histogram2D by binning npoints random points which are
Gaussian dsitributed"""
xmin, xmax = -20, 20
ymin, ymax = -3, 5
from random import gauss
from pylab import plot,contour,contourf
xmu, xsigma = 0, 6.0
ymu, ysigma = 0, 1.0
points = [(gauss(xmu,xsigma), gauss(ymu,ysigma)) for i in range(0,npoints)]
h = Histogram2D(points,xlimits={'min':xmin,'max':xmax,'nbins':20},
ylimits={'min':ymin,'max':ymax,'nbins':10})
try:
import pylab
h.plot_points()
h.plot_contours()
pylab.show
except:
pass
return h
def allplot(xb,yb,bins=30,fig=1,xlabel='x',ylabel='y'):
"""
Input:
X,Y : objects referring to the variables produced by PyMC that you want
to analyze. Example: X=M.theta, Y=M.slope.
Inherited from Tommy LE BLANC's code at astroplotlib|STSCI.
"""
#X,Y=xb.trace(),yb.trace()
X,Y=xb,yb
#pylab.rcParams.update({'font.size': fontsize})
fig=pylab.figure(fig)
pylab.clf()
gs = pylab.GridSpec(2, 2, width_ratios=[3,1], height_ratios=[1,3], wspace=0.07, hspace=0.07)
scat=pylab.subplot(gs[2])
histx=pylab.subplot(gs[0], sharex=scat)
histy=pylab.subplot(gs[3], sharey=scat)
#scat=fig.add_subplot(2,2,3)
#histx=fig.add_subplot(2,2,1, sharex=scat)
#histy=fig.add_subplot(2,2,4, sharey=scat)
# Scatter plot
scat.plot(X, Y,linestyle='none', marker='o', color='green', mec='green',alpha=.2, zorder=-99)
gkde = scipy.stats.gaussian_kde([X, Y])
x,y = numpy.mgrid[X.min():X.max():(X.max()-X.min())/25.,Y.min():Y.max():(Y.max()-Y.min())/25.]
z = numpy.array(gkde.evaluate([x.flatten(), y.flatten()])).reshape(x.shape)
scat.contour(x, y, z, linewidths=2)
scat.set_xlabel(xlabel)
scat.set_ylabel(ylabel)
# X-axis histogram
histx.hist(X, bins, histtype='stepfilled')
pylab.setp(histx.get_xticklabels(), visible=False) # no X label
#histx.xaxis.set_major_formatter(pylab.NullFormatter()) # no X label
# Y-axis histogram
histy.hist(Y, bins, histtype='stepfilled', orientation='horizontal')
pylab.setp(histy.get_yticklabels(), visible=False) # no Y label
#histy.yaxis.set_major_formatter(pylab.NullFormatter()) # no Y label
#pylab.minorticks_on()
#pylab.subplots_adjust(hspace=0.1)
#pylab.subplots_adjust(wspace=0.1)
pylab.draw()
pylab.show()
def jointplotx(X,Y,xlabel=None,ylabel=None,binsim=40,binsh=20,binscon=15):
"""
Plots the joint distribution of posteriors for X1 and X2, including the 1D
histograms showing the median and standard deviations. Uses simple method
for drawing the confidence contours compared to jointplot (which is wrong).
The work that went in creating this method is shown, step by step, in
the ipython notebook "error contours.ipynb". Sources of inspiration:
- http://python4mpia.github.io/intro/quick-tour.html
Usage:
>>> jointplot(M.rtr.trace(),M.mdot.trace(),xlabel='$\log \ r_{\\rm tr}$', ylabel='$\log \ \dot{m}$')
"""
# Generates 2D histogram for image
histt, xt, yt = numpy.histogram2d(X, Y, bins=[binsim,binsim], normed=False)
histt = numpy.transpose(histt) # Beware: numpy switches axes, so switch back.
# assigns correct proportions to subplots
fig=pylab.figure()
gs = pylab.GridSpec(2, 2, width_ratios=[3,1], height_ratios=[1,3], wspace=0.001, hspace=0.001)
con=pylab.subplot(gs[2])
histx=pylab.subplot(gs[0], sharex=con)
histy=pylab.subplot(gs[3], sharey=con)
# Image
con.imshow(histt,extent=[xt[0],xt[-1], yt[0],yt[-1]],origin='lower',cmap=pylab.cm.gray_r,aspect='auto')
# Overplot with error contours 1,2 sigma
# Contour plot
histdata, x, y = numpy.histogram2d(X, Y, bins=[binscon,binscon], normed=False)
histdata = numpy.transpose(histdata) # Beware: numpy switches axes, so switch back.
pmax = histdata.max()
cs=con.contour(histdata, levels=[0.68*pmax,0.05*pmax], extent=[x[0],x[-1], y[0],y[-1]], colors=['black','blue'])
# use dictionary in order to assign your own labels to the contours.
#fmtdict = {s[0]:r'$1\sigma$',s[1]:r'$2\sigma$'}
#con.clabel(cs, fmt=fmtdict, inline=True, fontsize=20)
if xlabel!=None: con.set_xlabel(xlabel)
if ylabel!=None: con.set_ylabel(ylabel)
# X-axis histogram
histx.hist(X, binsh, histtype='stepfilled',facecolor='lightblue')
pylab.setp(histx.get_xticklabels(), visible=False) # no X label
pylab.setp(histx.get_yticklabels(), visible=False) # no Y label
# Vertical lines with median and 1sigma confidence
yax=histx.set_ylim()
histx.plot([numpy.median(X),numpy.median(X)],[yax[0],yax[1]],'k-',linewidth=2) # median
xsd=scipy.stats.scoreatpercentile(X, [15.87,84.13])
histx.plot([xsd[0],xsd[0]],[yax[0],yax[1]],'k--') # -1sd
histx.plot([xsd[-1],xsd[-1]],[yax[0],yax[1]],'k--') # +1sd
# Y-axis histogram
histy.hist(Y, binsh, histtype='stepfilled', orientation='horizontal',facecolor='lightyellow')
pylab.setp(histy.get_yticklabels(), visible=False) # no Y label
pylab.setp(histy.get_xticklabels(), visible=False) # no X label
# Vertical lines with median and 1sigma confidence
xax=histy.set_xlim()
histy.plot([xax[0],xax[1]],[numpy.median(Y),numpy.median(Y)],'k-',linewidth=2) # median
ysd=scipy.stats.scoreatpercentile(Y, [15.87,84.13])
histy.plot([xax[0],xax[1]],[ysd[0],ysd[0]],'k--') # -1sd
histy.plot([xax[0],xax[1]],[ysd[-1],ysd[-1]],'k--') # +1sd
