In [2]:
import numpy as np
import emcee
In [3]:
def lnprob(x, mu, icov):
diff = x-mu
return -np.dot(diff,np.dot(icov,diff))/2.0
In [4]:
ndim = 50
means = np.random.rand(ndim)
cov = 0.5 - np.random.rand(ndim ** 2).reshape((ndim, ndim))
cov = np.triu(cov)
cov += cov.T - np.diag(cov.diagonal())
cov = np.dot(cov,cov)
In [5]:
icov = np.linalg.inv(cov)
icov
Out[5]:
array([[ 4569.32721243, 1649.02791658, -226.22457095, ...,
-652.00371283, -1578.85027295, 815.42949042],
[ 1649.02791658, 601.38602773, -81.04443964, ...,
-236.44785992, -571.30516305, 293.15943256],
[ -226.22457095, -81.04443964, 13.48812574, ...,
31.63774473, 79.38909582, -40.12141849],
...,
[ -652.00371283, -236.44785992, 31.63774473, ...,
98.5656679 , 224.94293113, -116.1402248 ],
[-1578.85027295, -571.30516305, 79.38909582, ...,
224.94293113, 548.11589098, -280.58151691],
[ 815.42949042, 293.15943256, -40.12141849, ...,
-116.1402248 , -280.58151691, 147.09798135]])
In [6]:
nwalkers = 250
p0 = np.random.rand(ndim * nwalkers).reshape((nwalkers, ndim))
In [7]:
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=[means, icov])
In [9]:
pos, prob, state = sampler.run_mcmc(p0, 100)
sampler.reset()
In [10]:
sampler.run_mcmc(pos, 1000)
Out[10]:
(array([[-4.38252517, -0.12927011, -0.65799899, ..., -0.75047415,
-1.50135091, -3.00523218],
[-0.47170324, 1.00096748, 1.31032215, ..., -1.64176611,
1.09574613, 1.42114428],
[ 2.84330073, 2.53204789, -1.40906236, ..., -1.635766 ,
5.58791383, 0.17650337],
...,
[-1.31801599, 2.6150006 , -3.25900568, ..., 1.20001899,
0.279695 , 2.12647992],
[-3.43968572, -1.50001115, 0.87406067, ..., -1.29098399,
0.69146265, 0.77521614],
[ 0.92668429, 0.14075794, 2.9250658 , ..., 2.15020666,
-2.34730894, 2.21119108]]),
array([-20.44016204, -15.93018203, -34.92443127, -32.53149756,
-29.35767343, -22.26642686, -26.59751834, -24.00759601,
-24.8145144 , -28.00052289, -27.29279042, -25.04097188,
-20.65887922, -28.32471941, -32.11380855, -22.32181377,
-30.97467267, -29.14523669, -32.64269937, -27.09043614,
-23.90436578, -28.24021607, -19.60994852, -17.52146353,
-19.87020662, -30.41581278, -31.81145085, -21.33532887,
-25.34169436, -26.50937885, -20.03365777, -23.16369494,
-25.24327148, -22.0445248 , -28.68941436, -24.27552636,
-31.25163736, -21.67854787, -26.57880513, -24.32108559,
-24.68187058, -29.18666122, -29.20979881, -20.38908801,
-26.71402771, -31.09816957, -20.97794585, -24.81643505,
-20.14746714, -16.33282516, -25.91389792, -29.8157881 ,
-32.31793551, -26.86539664, -25.864152 , -22.52932283,
-34.10743797, -20.67392235, -24.59768755, -32.28294861,
-20.59347301, -29.06290322, -25.28186815, -25.18455821,
-21.29891949, -23.42375366, -27.36278021, -28.40806364,
-17.28918509, -17.57542193, -30.57365779, -18.79335454,
-23.82379626, -35.84060987, -23.82445496, -21.46633662,
-25.24533054, -19.82725075, -24.6181942 , -20.60525109,
-36.06130757, -22.57090093, -25.43453495, -19.93423803,
-18.20979507, -26.37560006, -18.81928201, -17.61924258,
-20.30185229, -32.02818837, -20.90667502, -21.5568702 ,
-25.77900797, -26.63289171, -18.27807174, -21.90528288,
-27.19899213, -12.87128419, -22.4906286 , -22.74304009,
-23.89060472, -27.99559168, -24.54248591, -16.13122174,
-23.42979725, -24.71039283, -21.77078752, -25.12587394,
-26.09736149, -26.12573716, -26.76015823, -18.82988597,
-25.93691881, -23.43146718, -32.47529329, -34.50694545,
-24.83248053, -18.68260086, -26.3948937 , -26.71988499,
-22.57930953, -24.36239428, -36.732451 , -30.10683592,
-24.39281009, -19.56239939, -31.77086101, -19.39290212,
-24.34127802, -26.07550911, -29.94075693, -18.8282449 ,
-18.3163301 , -15.1416788 , -23.99695323, -29.87081338,
-28.81629178, -25.56989336, -22.62276468, -25.43924686,
-29.83583414, -31.37249925, -27.3561364 , -25.00816829,
-24.04311643, -26.84256855, -17.54878769, -32.46031677,
-28.88613985, -20.00769654, -26.82246078, -37.27304501,
-27.6564122 , -21.1905752 , -19.11489522, -21.17920531,
-19.63917538, -36.18955528, -23.44647365, -23.68230174,
-24.2113476 , -25.87055411, -22.4986743 , -40.08119738,
-22.1868643 , -27.99145022, -14.76592707, -32.46633228,
-25.32327606, -23.08186251, -31.40882566, -21.35587879,
-21.86766496, -32.65813026, -22.52584938, -19.0272046 ,
-30.01046291, -12.14271259, -25.15314206, -21.9542314 ,
-28.79191929, -28.74173605, -27.69596951, -23.27154949,
-16.23009467, -23.29300307, -23.79263334, -15.049948 ,
-15.14439066, -29.8642319 , -22.94990865, -21.12020107,
-25.58498475, -28.84859763, -27.84254598, -27.75429851,
-21.55343511, -27.74468185, -25.33416124, -33.93985454,
-29.69582755, -20.22113391, -19.90662096, -28.27892018,
-25.60788694, -26.97315761, -28.0511117 , -31.25324056,
-18.97104912, -22.17657527, -27.22994865, -24.16672909,
-24.61455907, -32.15113618, -24.34384692, -18.75756602,
-23.30051954, -23.7652008 , -20.08217266, -25.33671593,
-21.5650934 , -27.90383389, -32.81923935, -23.16799441,
-24.09414808, -36.46765762, -26.64843839, -25.24659945,
-20.04102056, -18.7566484 , -18.29639532, -32.25305711,
-19.45864788, -19.77602018, -17.0463771 , -35.03621681,
-28.4717697 , -23.60951537, -28.82782546, -22.35800126,
-24.14803114, -22.80726051, -22.02135501, -17.53867186,
-27.74706911, -22.18167168, -23.82342235, -22.59313912,
-19.34012418, -26.4104299 ]),
('MT19937', array([4294253032, 2742946155, 3788714503, 78253922, 2797834069,
4032618511, 2745118386, 2361316848, 2547418610, 1997461169,
2587650482, 771248799, 2994534183, 1324390891, 3786077347,
813846670, 91582724, 3485832874, 171134571, 4106520386,
502394110, 132488865, 2415036623, 2922306338, 1808770452,
1162103820, 1424524937, 853537557, 486748381, 984786319,
3571925405, 2673428819, 3980195266, 1135083245, 3321969629,
800511922, 1855444944, 1970858265, 4122546100, 924446379,
3227614543, 3663225082, 1728950976, 3901609186, 1711681040,
2308447232, 3023432287, 2827866452, 2033774250, 585024501,
249476350, 1798032258, 44933987, 936929107, 1925338671,
3575719851, 1292291576, 2492692104, 2309982139, 3654321734,
3912491797, 40064271, 4033354537, 3623451036, 1982612012,
792777145, 893178646, 3063862693, 2601912550, 2860816748,
290594705, 1974334911, 2881208814, 104942756, 1947339540,
2356261310, 1980872328, 523427896, 2925303993, 3592512117,
3825927416, 3255662568, 3628280426, 1529017728, 958323035,
788313232, 1619892024, 3202409527, 1822666817, 2800138062,
2215190088, 545009342, 150640288, 4150623776, 2493156337,
4130307961, 813704360, 3809089477, 1438118766, 3242583933,
2235745233, 2020336511, 4068868989, 1647473291, 3967891450,
1043843951, 1452200518, 3578459603, 2066660636, 3570593548,
89978921, 1302534492, 644644039, 428266960, 1692692466,
2783000373, 2313488873, 899531025, 32058864, 1459302603,
3085379891, 2308292757, 3770813352, 3918707963, 2088391055,
1636007447, 1309628879, 1699663273, 4267307805, 2348759564,
1326167508, 2870290644, 3783957232, 435184575, 2366491949,
3943232606, 3786559295, 2472679092, 2622601986, 311791359,
2260744215, 2665966919, 1101777404, 1889833219, 2073753998,
3738602425, 4246283916, 1370450325, 2218115840, 2991805807,
2338119848, 1259808171, 4229558315, 1176096208, 1671976954,
2670271469, 2900224472, 4250237457, 4241371945, 2967742730,
3958497955, 3541685773, 3890387595, 3396143287, 3135755617,
2767973578, 4030414755, 3624717496, 1886379792, 2160408373,
1422041627, 1055908283, 1005052828, 4243200104, 1222121597,
2626122950, 2764658651, 3946182077, 962515251, 2244874949,
907530827, 960893914, 71671558, 2305843210, 2805728928,
3693828939, 3314535516, 612186991, 2277819824, 4035202981,
3448271322, 1690401422, 3045438369, 4208726227, 2237681000,
1128195188, 586316349, 607421704, 3354506235, 2840284061,
2753902093, 1193186109, 3126175514, 4178773794, 3128283492,
3507395550, 1094034228, 2369895259, 1706699710, 2351217637,
481948368, 3620060818, 1147308930, 935886246, 767464588,
3629987589, 2156402470, 2936987204, 2238911830, 394210852,
2506312143, 2031138246, 3645727856, 261583387, 2736688512,
1036582662, 432639696, 3748892527, 3267991001, 3861778883,
2502987562, 2687780185, 76576190, 2925691915, 4055630571,
1836969375, 1347886247, 3317807830, 3741171813, 131091152,
1441048883, 2347371913, 3984140232, 2248107609, 3113629325,
2221145039, 3525817007, 2403819698, 3765192963, 986169158,
2267680969, 551229206, 2659687916, 1379311972, 1871804388,
3715984688, 3936493212, 1659025271, 2210563516, 1711879250,
363591566, 1106472074, 3883850677, 3143907926, 428229805,
2840810167, 2235200854, 3893791684, 777844506, 1755655064,
451774848, 3657660, 4096259709, 1503406683, 2402147231,
696324406, 3660747289, 692630762, 975015540, 1401810696,
2643841724, 1051384709, 3146707898, 2006265687, 1909029155,
2239435870, 1955014619, 1338201147, 2484857213, 2084654148,
4194647956, 2238986692, 2425013783, 805596610, 3104361927,
1024343695, 1867420893, 173877899, 1065678361, 2768732398,
1120941714, 2895058597, 1633971437, 2000396444, 964008197,
2294464675, 980660531, 10847195, 309740515, 3206233836,
1083834025, 1203796556, 1269670736, 1330159600, 2795866541,
850749544, 1313094286, 214229714, 2380945322, 775141243,
2799698278, 3741048548, 4282619521, 1293680655, 3199063123,
3204113336, 3451689525, 2225344218, 1557888710, 1005862345,
4257506281, 77341925, 1326125293, 757475434, 1510340782,
1579615791, 883798896, 327682680, 3800144595, 2274496154,
1233225341, 3666344909, 1452095741, 4022862874, 2536789486,
3154441050, 292426748, 4148987056, 2954894233, 4261860659,
1189116051, 3067539072, 1816365736, 3912915826, 1139080356,
532223078, 3112655673, 166465272, 2928902776, 4116863544,
1652643734, 896896575, 2111470494, 2575079547, 2061696900,
1226476348, 1440423299, 361347860, 4027863796, 1959636204,
2662005065, 2745237243, 1569510059, 1690156358, 4030971614,
254775911, 205981275, 3302618561, 2121576132, 1743595501,
4026822129, 1402127878, 715220223, 1091604744, 1503522738,
3410236050, 2990532448, 3939303931, 3087272874, 3664884170,
2377077754, 83311228, 2194266081, 2983448677, 574728806,
4106934418, 4230177757, 2984173702, 3189175601, 2135013083,
1268708441, 3241977338, 1222723522, 640497069, 3451556695,
324710719, 2684397736, 4221751005, 2183469269, 2864478049,
1115977449, 2226970172, 2921444217, 3559646360, 3730773056,
3275546844, 922801, 979419361, 131677171, 4140815299,
2871108555, 1129537665, 1631382704, 2520302358, 1340386787,
50751229, 1591623787, 516745613, 2612562415, 2347111285,
3496131080, 671526414, 2989335696, 3481066304, 2945398474,
3914572454, 3444167693, 2332772632, 3052496683, 3329285444,
2881411707, 1602543061, 1782816360, 4011261914, 3540444867,
2557751800, 2183164581, 729670827, 1984446020, 2439247051,
1682251006, 2204844934, 1910191394, 3684311021, 2783202223,
2543332073, 954444831, 4260908001, 574455460, 270846335,
948135562, 1797777543, 3974534910, 3024710691, 3509390493,
2504164689, 3458935062, 146256678, 1848224385, 3495691509,
1074354715, 963337110, 3123424562, 2449841908, 901124442,
1755676290, 2747487818, 2703704033, 3856543459, 1615190149,
1660898267, 818286002, 1035247329, 2224203613, 1511635111,
2230692419, 2050877706, 3710726494, 3800575392, 2243601312,
1779715095, 477512348, 2028287573, 3521839510, 1835835696,
964644967, 2314186323, 3938721026, 4207242688, 671699961,
3613344227, 1480659831, 3376835260, 3465056784, 627825074,
2156016573, 4140403107, 2005099745, 674096604, 224261478,
3347018690, 3758356434, 1298470202, 1077440761, 3229880788,
3676497288, 2647051617, 2811649003, 3982783181, 4245190689,
3539164744, 3062500828, 3910324189, 590378453, 3577332401,
2948521205, 3937414694, 874371936, 1984595670, 2134048317,
1370658419, 2177318468, 956991263, 255507841, 246935664,
244107291, 534266697, 3082647569, 3974110539, 290258501,
23077744, 1067243386, 3091962918, 3745018503, 576042084,
3234315063, 317099129, 50962975, 3751269971, 1371412593,
3847748263, 3092140383, 337881515, 3743909801, 1153521106,
2554256160, 3342302308, 2050878160, 2259675412, 304960117,
3477178501, 2652796778, 3653559535, 1458512091, 2032889941,
947609055, 1741718739, 3669355889, 1008615690, 1397165331,
121204955, 145605088, 322778945, 2825871901, 1434757418,
571877489, 84276290, 2589793689, 2870989051, 3875630032,
2377272178, 3365998484, 4184537351, 1322307475, 2839363792,
2788496946, 1032639736, 2189903571, 1715261858, 3230840133,
1022667470, 174701639, 2990193529, 3374999811, 749452716,
94864278, 3357373591, 2665633281, 1769154837, 3162146722,
2441196137, 3882252782, 2934229737, 75086702, 3743697063,
2517555464, 1937168755, 4248417510, 510859290, 4167367519,
584607907, 1986902510, 2750609516, 3095212358, 485206545,
15284118, 688738610, 736015505, 4197999190, 333869187,
502731701, 2024768388, 3304007369, 3285839657], dtype=uint32), 542, 0, 0.0))
In [11]:
%matplotlib inline
import matplotlib.pyplot as pl
for i in range(ndim):
pl.figure()
pl.hist(sampler.flatchain[:,i], 100, color="k", histtype="step")
pl.title("Dimension {0:d}".format(i))
pl.show()
/scratch/seismo/bellinger/python/anaconda3-5.3.0/lib/python3.7/site-packages/matplotlib/pyplot.py:522: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
max_open_warning, RuntimeWarning)
In [12]:
print("Mean acceptance fraction: {0:.3f}"
.format(np.mean(sampler.acceptance_fraction)))
Mean acceptance fraction: 0.190
In [ ]:
Content source: earlbellinger/asteroseismology
Similar notebooks: