注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

一粒浮尘

飘渺虚无

 
 
 

日志

 
 

a python script to plot contours  

2014-04-08 21:30:29|  分类: python |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
from numpy import *
from matplotlib.pylab import *
from scipy.stats import gaussian_kde
import kdestats as kde

chain = loadtxt('unbias/test_with_MB_no_error_int_fix_H0_sim_0.txt.filtered')
chainx= loadtxt('unbias/test_with_MB_no_error_int_fix_H0_sim_0.txt')

#    thin the original chain
chainx= chainx[chainx[:,7] > -29, :]
chainx= chainx[chainx[:,7] < -28, :]
ID    = chainx[:,0] > 0    #    create bool ID
for i in range(chainx.shape[0]):
    ran = rand()
    if ran < 0.1:
        ID[i] = True
    else:
        ID[i] = False

chainx    = chainx[ID,:]

#    plot 2d contours
N = 7
for i in range(N-1):
    for j in range(i+1, N):
        print 'subplot(%d, %d)'%(i,j)
        fig = subplot(N,N,N*j+i+1)

        kdehist     = kde.kdehist2(chain[:,i], chain[:,j], [40,40])
        kdehistx     = kde.kdehist2(chainx[:,i], chainx[:,j], [40,40])
        clevels     = kde.confmap(kdehist[0], [.6827,.9545])
        clevelsx     = kde.confmap(kdehistx[0], [.6827,.9545])

        fig.contour(kdehist[1], kdehist[2], kdehist[0], clevels, linestyles='solid')
        fig.contour(kdehistx[1], kdehistx[2], kdehistx[0], clevelsx, linestyles='dotted')

        xmin = chainx[:,i].min()
        xmax = chainx[:,i].max()
        ymin = chainx[:,j].min()
        ymax = chainx[:,j].max()

        X,Y         = mgrid[xmin:xmax:100j, ymin:ymax:100j]
        position    = vstack([X.ravel(), Y.ravel()])
        values         = vstack([chainx[:,i].T,chainx[:,j].T])
        kernel        = gaussian_kde(values)
        Z            = reshape(kernel.evaluate(position).T, X.shape)
        fig.imshow(rot90(Z), interpolation='gaussian', cmap=cm.Greys, extent=[xmin,xmax,ymin,ymax], aspect='auto')

        if j != (N-1):
            fig.axes.get_xaxis().set_visible(False)

        if i != 0:
            fig.axes.get_yaxis().set_visible(False)

        fig.set_xlim([xmin, xmax])
        fig.set_ylim([ymin, ymax])
        fig.xaxis.set_major_locator(MaxNLocator(3))
        fig.yaxis.set_major_locator(MaxNLocator(3))  

#    plot margionalized 1d distribution
for i in range(N):
    print '%d th 1d-density'%i
    fig = subplot(N,N,N*i+i+1)
    density    = gaussian_kde(chain[:,i])
    densityx= gaussian_kde(chainx[:,i])
    t = linspace(chain[:,i].min(), chain[:,i].max(), 50)
    tx= linspace(chainx[:,i].min(), chainx[:,i].max(), 50)
    plot(t, density(t), 'b-')
    plot(tx, densityx(tx), 'r:')
    
    xmin = chainx[:,i].min()
    xmax = chainx[:,i].max()
    fig.set_xlim([xmin, xmax])

    fig.xaxis.set_major_locator(MaxNLocator(3))
    fig.yaxis.set_major_locator(MaxNLocator(3))

    if i != (N-1):
        fig.axes.get_xaxis().set_visible(False)

    if i != 0:
        fig.axes.get_yaxis().set_visible(False)

show()
  评论这张
 
阅读(416)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017