In [1]:
import numpy as np
import pylab as plt
%matplotlib inline
In [2]:
#2d gaussian
def lnprob(theta,sigma):
x, y = theta[0], theta[1]
#linear version 1/np.sqrt(2 * np.pi * sigma**2) * np.exp(-(x**2+y**2)/(2. * sigma**2))
return -0.5 * np.sum(((x)**2 + (y)**2)/(2. * sigma**2))
In [3]:
import emcee
In [4]:
# set up number of dimensions, number of walkers, and initial state
ndim, nwalkers = 2, 100
p0 = [10+np.random.rand(ndim) for i in range(nwalkers)]
#print(p0)
In [5]:
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=[0.5])
sampler.run_mcmc(p0, 1000)
Out[5]:
(array([[-0.26947342, -0.21850108],
[ 0.14104332, -0.59469276],
[ 1.16772509, 0.60212308],
[ 1.74297206, -0.4923362 ],
[ 0.96993967, -0.21520631],
[-0.72669737, -0.12321428],
[-0.71983118, 0.51652562],
[ 0.56242918, -0.06220203],
[-0.51979727, 0.09881683],
[ 0.3767379 , -1.46723498],
[ 0.54274494, 1.00169786],
[-0.3129038 , 0.21690838],
[ 0.14815222, 1.1096948 ],
[ 0.15174132, 0.52440467],
[ 0.55006864, 1.19540035],
[ 1.26344948, -0.82170239],
[-1.36511439, 0.30938949],
[-0.08895898, 0.19277736],
[-0.35781709, 0.48521987],
[-0.19825114, 0.85092948],
[ 0.72978842, -0.05551697],
[ 0.2022186 , -0.58143419],
[ 0.60633443, 0.1448309 ],
[-0.08825784, -0.22908696],
[-0.72718663, -0.34689252],
[-0.74877302, -0.56478272],
[ 0.48265009, -0.33981233],
[-0.21672466, -0.21400301],
[ 0.22429966, 1.14625239],
[ 0.4925886 , -0.26893401],
[ 0.30556108, 0.52542065],
[ 0.19251603, 0.54774326],
[ 0.19251415, -0.55864976],
[ 0.35228346, 0.64338629],
[ 0.27014052, -0.80451341],
[ 0.43358257, 0.47773112],
[ 1.06598138, 0.2674888 ],
[-0.72859838, -0.83340589],
[ 0.61331502, -1.25567567],
[ 1.52540358, 0.65800524],
[-1.67743209, -0.49230401],
[ 0.32414545, 0.67300311],
[ 0.11412731, -1.23133586],
[-0.10368923, 0.23500953],
[ 0.73021574, -0.5339627 ],
[ 0.02719722, -0.05180397],
[ 0.774923 , 0.40036739],
[-0.22084093, -0.76377042],
[-1.03317355, 0.34634083],
[ 0.15099627, -0.16738194],
[-0.36871865, 0.5148968 ],
[-1.56559407, -0.6168216 ],
[ 0.68854128, 0.34932394],
[-0.79829192, -0.97127212],
[-0.09448415, 0.67744614],
[ 0.18837253, -0.05437237],
[-0.60617616, 0.50757636],
[-0.35513248, -0.5247915 ],
[-0.36556421, 1.23240278],
[ 0.33420941, 0.02677793],
[ 0.69923625, -0.25075284],
[-0.62383196, -0.67864724],
[-0.14691878, -1.0561328 ],
[-0.29258142, 0.2539456 ],
[-0.02070014, 0.48425666],
[-0.03035314, -0.01183446],
[ 0.33051762, 1.0117574 ],
[ 1.17425521, 0.01786306],
[-0.68402106, 1.50245507],
[ 0.21808119, 0.54461488],
[ 0.14828291, 0.03635819],
[ 0.09987575, 0.84704261],
[ 0.02897883, -0.70417089],
[ 0.46608851, 0.65790296],
[ 0.87743779, 0.21957558],
[ 0.16415359, -0.42459478],
[-0.12939446, 0.37399356],
[-1.01379682, -0.70202772],
[-0.09902164, -0.08477536],
[-0.57787587, 1.08843895],
[ 0.5353784 , -0.02166332],
[ 0.10160048, -0.0779501 ],
[-0.82957523, -0.7327819 ],
[ 0.57965025, 0.00992026],
[-0.01695569, -0.57539799],
[-1.19407263, 0.42493821],
[-0.69127764, 0.27503269],
[-1.260881 , -0.76239769],
[-0.42722828, -0.74061864],
[ 0.11093326, 0.72510112],
[-0.24805775, 0.65726605],
[-0.88999238, 1.04032058],
[-1.41772254, -0.98470291],
[-0.49282503, -0.48352346],
[ 0.11859199, 0.7338821 ],
[ 1.67354802, -0.23727322],
[-0.81793279, 0.1297868 ],
[ 0.01580233, 0.7908139 ],
[-0.94610026, -0.94596928],
[ 0.94359291, 1.9022856 ]]),
array([ -1.20358647e-01, -3.73552694e-01, -1.72613409e+00,
-3.28034652e+00, -9.87096727e-01, -5.43270827e-01,
-7.84955642e-01, -3.20195675e-01, -2.79953968e-01,
-2.29470994e+00, -1.29797067e+00, -1.44958030e-01,
-1.25337163e+00, -2.98025684e-01, -1.73155750e+00,
-2.27149941e+00, -1.95925916e+00, -4.50768108e-02,
-3.63471391e-01, -7.63384493e-01, -5.35673270e-01,
-3.78958082e-01, -3.88617431e-01, -6.02702835e-02,
-6.49134824e-01, -8.79640561e-01, -3.48423530e-01,
-9.27668660e-02, -1.36420489e+00, -3.14969034e-01,
-3.69434437e-01, -3.37085106e-01, -3.49151257e-01,
-5.38049553e-01, -7.20217732e-01, -4.16220869e-01,
-1.20786657e+00, -1.22542098e+00, -1.95287669e+00,
-2.75982700e+00, -3.05614166e+00, -5.58003460e-01,
-1.52921305e+00, -6.59809337e-02, -8.18331181e-01,
-3.42334027e-03, -7.60799699e-01, -6.32115968e-01,
-1.18739956e+00, -5.08165872e-02, -4.01072159e-01,
-2.83155368e+00, -5.96116310e-01, -1.58063952e+00,
-4.67860532e-01, -3.84405632e-02, -6.25083308e-01,
-4.01525196e-01, -1.65245380e+00, -1.12412984e-01,
-5.51808318e-01, -8.49728400e-01, -1.13700163e+00,
-1.50092256e-01, -2.34933005e-01, -1.06136746e-03,
-1.13289493e+00, -1.37919438e+00, -2.72525603e+00,
-3.44164775e-01, -2.33097382e-02, -7.27456353e-01,
-4.96696418e-01, -6.50074805e-01, -8.18110507e-01,
-2.07227126e-01, -1.56614109e-01, -1.52062690e+00,
-1.69921471e-02, -1.51863987e+00, -2.87099335e-01,
-1.63988767e-02, -1.22516438e+00, -3.36092821e-01,
-3.31370342e-01, -1.60638194e+00, -5.53507750e-01,
-2.17107113e+00, -7.31039967e-01, -5.38077815e-01,
-4.93531309e-01, -1.87435334e+00, -2.97957702e+00,
-4.76671452e-01, -5.52646991e-01, -2.85706157e+00,
-6.85858668e-01, -6.25636344e-01, -1.78996356e+00,
-4.50905808e+00]),
('MT19937', array([1159659705, 1339467057, 1721083420, 3526080451, 84074072,
2153743248, 1872770706, 3480158278, 1734991940, 606016387,
2406728178, 594971958, 3461275893, 3551466588, 2512930667,
1095678051, 59185319, 3668614240, 3178053843, 2441618752,
943797906, 2564611594, 1490473649, 4043575317, 3462857997,
1796949804, 2775388337, 1233506658, 1788261737, 3723381686,
2703040593, 3354288910, 2487344272, 246243510, 4254076894,
3532592456, 1714513617, 3195334606, 3491745202, 950612154,
2113748466, 3764894377, 3784449354, 1056358780, 76153220,
1421680570, 2112823832, 3771630645, 1080616219, 2672560933,
291010686, 2660007578, 4053840369, 3086838384, 3017098531,
2815792891, 3764674878, 4057261372, 1072524337, 3335291567,
1078529654, 846855808, 3278043901, 2619560265, 3668228590,
3421731818, 3210304275, 1943520061, 177356622, 3277019294,
4086100115, 638107852, 49446194, 1499921969, 3026923937,
422414923, 1845902977, 1962707008, 1712366549, 2257043799,
4197157430, 3936637922, 2157565133, 2015292223, 2692270840,
74698108, 784940165, 3437839326, 2495154730, 2074975689,
1228369189, 3632862706, 2328663322, 2783251332, 1507709094,
3865671721, 3747750791, 2075444302, 537210543, 2302625751,
3327340610, 2700368470, 327361030, 743452757, 267996025,
1581137187, 4201218646, 3148783500, 3013058234, 1038374744,
2576666612, 220926526, 3680154594, 1194474618, 2697265928,
1954180407, 916859791, 736849588, 509496121, 4213757978,
791259528, 2089935602, 924057727, 1111810469, 4163042055,
1164788881, 2409695729, 2059306058, 3292570800, 456655920,
2612806935, 875264876, 3014231460, 663348339, 3900038840,
2610702330, 2789868869, 1483910097, 1019734471, 3620391045,
3681961799, 3400206220, 117798146, 447832829, 1432382381,
1410106873, 296412957, 2059732311, 3962300950, 2662561885,
136581771, 1831810269, 3883645736, 1672800283, 1635495891,
2878846817, 2991470428, 4111155238, 1584221316, 1409744428,
4027909701, 2045784019, 666198445, 3878401291, 3668196312,
3856996368, 1900586793, 3583499101, 1751685106, 2624633992,
1978605035, 432014020, 3170689313, 3078320327, 4096799075,
1570543007, 2692466851, 86565054, 3131120191, 1874942512,
1661916175, 1561910638, 2245839617, 103649041, 2101551214,
1846642167, 140161088, 1861352757, 1275926755, 3518884519,
3224953507, 422590647, 791491053, 3737081906, 2135503242,
4291918831, 4198513554, 3843825079, 4143426099, 3313417496,
1347716700, 1372026265, 495663198, 1278463339, 3255235768,
1647975816, 4148103204, 2936851445, 3579597817, 1088602402,
2598315331, 2173438853, 2121964710, 196434925, 2430546574,
2065946741, 2211735587, 4030405441, 2411566065, 2136459416,
3873011400, 3586110035, 2803183574, 3897925522, 2104852619,
2523395103, 390827765, 3439692681, 1450377861, 253517019,
229399917, 2673140883, 713496405, 2354584215, 3251228465,
421441315, 4186378321, 2817205607, 4221182419, 716042404,
4053340586, 3536206411, 4156393997, 3284729479, 2099699702,
1427443375, 88689508, 1225303164, 2356719910, 4042855090,
834994682, 550218351, 23639448, 2945410702, 914754050,
3248464007, 2710000493, 3013569128, 2614463170, 2428468161,
3064557040, 3976565593, 635867608, 411457663, 4252027675,
1819731115, 3427824731, 3445593843, 1755400854, 1842729038,
2098712502, 3052017438, 3898260035, 2533950150, 4219967998,
480383599, 1510597457, 901253871, 4074143709, 325440753,
372690993, 3355844626, 3128236581, 3586730122, 527306922,
2037545790, 246067249, 3302343390, 3025187139, 4172799850,
584506894, 3277781832, 2408640884, 3767761170, 125627345,
3570017511, 4130826729, 1615292620, 1246879326, 1035351408,
558520527, 3084276957, 301232803, 3769287268, 1976947424,
1831866502, 3384374471, 671740853, 4023771638, 1396302821,
2678306489, 4096982280, 1264236973, 538396444, 4021852912,
2055491440, 2237931926, 2518744216, 2722101810, 1970572776,
632816861, 2080307234, 1842940400, 334260914, 334866466,
754078556, 3994217984, 1462787727, 2629859174, 2096742851,
1836962930, 180301854, 272548457, 2745879665, 4181215531,
3809884110, 3083694818, 145143934, 2461435347, 3249583206,
538832040, 937188819, 1918577294, 2384008021, 1140672753,
189614592, 4263542061, 899592956, 2181207792, 2267248814,
1963235771, 180149934, 3147551562, 37733444, 321462687,
2454597172, 3339018551, 389106353, 2285743239, 2351497080,
3257160495, 4230890669, 3663231891, 682477544, 284492344,
382872823, 2904771529, 1671824067, 3298812481, 3515124819,
3805810547, 993014737, 2080913113, 1611400687, 566115964,
1197542394, 456221232, 310959073, 1452501269, 3217142686,
3894432062, 3060274224, 3135706583, 2267082865, 1958389294,
327115988, 1891539775, 3456851150, 2214719591, 3416803126,
679602750, 716087890, 1880093704, 2817193542, 2100294893,
2367664073, 1893234990, 3992820414, 1358795339, 856561443,
1308797856, 2389041835, 866779676, 2279041722, 4270188564,
1975667379, 781024296, 3105157762, 3318106183, 1087711264,
668763188, 1342034372, 354857932, 3046028138, 2933423556,
2733472808, 45783868, 4084880232, 343887657, 2670143250,
616724866, 3048253197, 3791227043, 1791113053, 861847905,
4020485112, 3542597966, 491794628, 2875824002, 2996583476,
1752040757, 164289595, 312848029, 2914924138, 3673391765,
366540646, 1941675876, 2600473069, 1553907530, 3315453576,
3707928356, 917814388, 921057515, 1173135393, 927156620,
1489945302, 2289589040, 3546717742, 3560596790, 547383114,
1317163326, 1169747033, 217033752, 877231059, 575091476,
574722402, 43962450, 1334373214, 1658059264, 1715532284,
32240334, 993791392, 2595537884, 3626783063, 3609145679,
132185162, 3669450232, 1831764850, 1646774837, 1836759350,
3460173122, 1823137005, 297933785, 3212452777, 4244388594,
167586316, 2900848683, 1903224330, 3899839740, 850669957,
898953533, 4120874693, 594258725, 3336985596, 3997746803,
2778467162, 1934404194, 4116861911, 624962615, 833746335,
2863907128, 102165527, 3038855647, 2183911521, 2321419864,
2007417883, 272429920, 1922776900, 2808026282, 1418264259,
866314534, 3624395404, 2272620246, 3761483436, 661368998,
708071634, 1890395683, 1352650440, 2246238602, 1772599472,
2459640745, 3165037827, 2082759043, 3664969203, 1165487692,
39299256, 2692651223, 3627653438, 568174685, 1231600027,
1897229245, 3202411667, 1199828201, 2524418309, 1192588485,
1666391747, 724694042, 1597273893, 1564200335, 1438920242,
1555528722, 3429377936, 2195357879, 539848986, 3538663573,
2233495804, 4188094723, 246685994, 2752680369, 2465197035,
2994622137, 1320805426, 913797755, 515248866, 1345778619,
2015861777, 70147203, 4271678818, 1963089032, 3631383366,
1092971791, 2405820801, 3881516705, 3538221870, 2744494817,
3415428486, 2376117367, 1151129041, 3835179714, 1678069509,
2854030433, 813922559, 3209959318, 1965688007, 1247332174,
2728504008, 3029714249, 4168899934, 2904060891, 2593020373,
385698308, 4160246665, 2364790542, 3929924994, 951132696,
2349186020, 2858669308, 3823909569, 580897040, 1735438457,
1427629100, 2158486820, 3658062802, 611389355, 2000261867,
927873487, 4087875064, 1455937473, 4256307559, 2593383933,
3788724052, 3617452286, 1736965832, 182929474, 1123820943,
3610480890, 54152311, 3331162757, 515344560, 253624633,
68556256, 2765881777, 3739380565, 2935554497, 1534113743,
4215815579, 1374286845, 1041371175, 3373380497, 3772335043,
3353864671, 3168033464, 3896721909, 3219657925, 578533742,
2143246658, 458097735, 4287445163, 2446699732, 578422860,
6430587, 1410708115, 3354408793, 2338277805], dtype=uint32), 65, 0, 0.0))
In [6]:
import corner
corner.corner(sampler.flatchain, labels=['$x$', '$y$'],quantiles=[0.16, 0.5, 0.84], truths=[0.,0.])
#, levels=(1-np.exp(-0.5),))
//anaconda/lib/python3.4/site-packages/matplotlib/collections.py:650: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
if self._edgecolors_original != str('face'):
//anaconda/lib/python3.4/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
if self._edgecolors == str('face'):
Out[6]:
In [7]:
sampler.chain.shape
sampler.flatchain.shape
Out[7]:
(100000, 2)
In [8]:
[plt.plot(sampler.lnprobability[i,:]) for i in range(100)]
plt.xlabel('step number')
plt.ylabel('lnP')
Out[8]:
<matplotlib.text.Text at 0x10dc4b978>
In [9]:
[plt.plot(sampler.chain[i,:,0]) for i in range(100)]
#[plt.plot(sampler.chain[i,:,0]) for i in range(sampler.chain.shape[0])]
plt.xlabel('step')
plt.ylabel('x')
Out[9]:
<matplotlib.text.Text at 0x10ed37550>
In [10]:
[plt.scatter(sampler.chain[i,:,0], sampler.lnprobability[i,:]) for i in range(1)]
plt.xlabel('x')
plt.ylabel('lnp')
Out[10]:
<matplotlib.text.Text at 0x10dec4da0>
//anaconda/lib/python3.4/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
if self._edgecolors == str('face'):
In [12]:
sampler.get_autocorr_time()
Out[12]:
array([ 34.88753109, 37.78685551])
In [ ]:
In [ ]:
In [ ]:
Content source: dweisz/ay250_fall2016
Similar notebooks: