bayes.py 文件源码

python
阅读 25 收藏 0 点赞 0 评论 0

项目:nmmn 作者: rsnemmen 项目源码 文件源码
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
评论列表
文章目录


问题


面经


文章

微信
公众号

扫码关注公众号