Data analysis with MCMC

Transcript

1. data analysis with Markov chain Monte Carlo Dan Foreman-Mackey CCPP@NYU
2.    dan.iel.fm dfm @exoplaneteer Dan Foreman-Mackey
3. I'm a grad student in Camp Hogg at NYU David W. Hogg
4. I'm an engineer in Camp Hogg at NYU David W. Hogg
5. David W. Hogg I build tools for data analysis in astronomy(mostly)
6. emceethe MCMC Hammer introducing arxiv.org/abs/1202.3665 dan.iel.fm/emcee it's hammer time!
7. I work on the Local Group exoplanets variable stars image modeling calibration
8. I work on the Local Group exoplanets variable stars image modeling calibration not writing papers
9. I work on the Local Group exoplanets variable stars image modeling calibration not writing papers
10. Physics Data The graphical model of my research. a sketch of
11. Physics Data The graphical model of my research. a sketch of the generative view of data analysis
12. 2 = N X n=1 [ yn ( m xn + b )]2 2 n a line ˆ yn = m xn + b + ✏n ✏n ⇠ N(0; 2 n ) synthetic data: noise model: p( { yn, xn, 2 n } | m, b) / exp ✓ 1 2 2 ◆ for example
13. Physics The graphical model of my research. a sketch of inference Data
14. p(data | physics) likelihood function/generative model p(physics | data) / p(physics) p(data | physics) posterior probability
15. but physics is everything. including things we don't care about! good/bad news:
16. a more realistic model: p(D | ✓, ↵) data cool physics garbage
17. what if we underestimated our error bars? N Y n=1 1 p 2 ⇡ ( 2 n + 2 ) exp  1 2 [yn (m xn + b)] 2 2 n + 2 p ({ yn, xn, 2 n } | m, b, 2) = "physics" "jitter" (not physics) (by some additive variance)
18. p(D | ✓) / Z p(↵) p(D | ✓, ↵) d↵ marginalize. Do The Right Thing™
19. Z d 2 p( 2 ) N Y n=1 1 p 2 ⇡ ( 2 n + 2 ) exp  1 2 [yn (m xn + b)] 2 2 n + 2 p ({ yn, xn, n } | m, b ) = in our example BOOM!
20. What do you really want?
21. What do you really want? Ep[f(✓)] = 1 Z Z f(✓) p(✓) p(D | ✓) d✓
22. probably. What do you really want? Ep[f(✓)] = 1 Z Z f(✓) p(✓) p(D | ✓) d✓
23. Ep[f(✓)] = 1 Z Z f(✓) p(✓) p(D | ✓) d✓ p(D | ✓) / Z p(↵) p(D | ✓, ↵) d↵ marginalization expectation
24. Ep[f(✓)] = 1 Z Z f(✓) p(✓) p(D | ✓) d✓ p(D | ✓) / Z p(↵) p(D | ✓, ↵) d↵ marginalization expectation large number of dimensions
25. Ep[f(✓)] = 1 Z Z f(✓) p(✓) p(D | ✓) d✓ p(D | ✓) / Z p(↵) p(D | ✓, ↵) d↵ marginalization expectation large number of dimensions whoa!
26. Ep[f(✓)] = 1 Z Z f(✓) p(✓) p(D | ✓) d✓ p(D | ✓) / Z p(↵) p(D | ✓, ↵) d↵ marginalization expectation This is HARD. (in general) large number of dimensions whoa!
27. OK now that we agree…
28. Z f( x ) p( x ) d x ⇡ 1 N N X n=1 f( xn) xn ⇠ p( x ) where: error: number of independent samples as you learned in middle school / 1 p N0
29. Z f( x ) p( x ) d x ⇡ 1 N N X n=1 f( xn) xn ⇠ p( x ) where: error: number of independent samples as you learned in middle school / 1 p N0
30. MCMC draws samples from a probability function and all you need to be able to do is evaluate the function (up to a constant)
31. Metropolis–Hastings
32. Metropolis–Hastings in an ideal world
33. start here perhaps Metropolis–Hastings in an ideal world
34. propose a new position Metropolis–Hastings in an ideal world x 0 ⇠ q( x 0; x )
35. propose a new position Metropolis–Hastings in an ideal world x 0 ⇠ q( x 0; x ) accept? p(accept) = min ✓ 1, p( x ) p( x 0) q( x ; x 0) q( x 0; x ) ◆
36. propose a new position Metropolis–Hastings in an ideal world x 0 ⇠ q( x 0; x ) definitely. accept? p(accept) = min ✓ 1, p( x ) p( x 0) q( x ; x 0) q( x 0; x ) ◆
37. propose a new position Metropolis–Hastings in an ideal world x 0 ⇠ q( x 0; x ) definitely. accept? p(accept) = min ✓ 1, p( x ) p( x 0) q( x ; x 0) q( x 0; x ) ◆ only relative probabilities
38. Metropolis–Hastings in an ideal world x 0 ⇠ q( x 0; x )
39. Metropolis–Hastings in an ideal world x 0 ⇠ q( x 0; x )
40. Metropolis–Hastings in an ideal world x 0 ⇠ q( x 0; x ) accept? p(accept) = min ✓ 1, p( x ) p( x 0) q( x ; x 0) q( x 0; x ) ◆
41. not this time. Metropolis–Hastings in an ideal world x 0 ⇠ q( x 0; x ) accept? p(accept) = min ✓ 1, p( x ) p( x 0) q( x ; x 0) q( x 0; x ) ◆
42. Metropolis–Hastings in an ideal world
43. double count! Metropolis–Hastings in an ideal world
44. Metropolis–Hastings in an ideal world
45. Metropolis–Hastings in an ideal world
46. Metropolis–Hastings in the real world
47. Metropolis–Hastings in the real world
48. Metropolis–Hastings in the real world
49. Metropolis–Hastings in the real world
50. Metropolis–Hastings in the real world
51. Metropolis–Hastings in the real world the Small Acceptance Fraction problem
52. Metropolis–Hastings in the real world
53. Metropolis–Hastings in the real world
54. Metropolis–Hastings in the real world the Huge Acceptance Fraction problem
55. a brief aside.
56. YOUR LIKELIHOOD FUNCTION IS EXPENSIVE
57. Metropolis–Hastings in the real world
58. YOUR LIKELIHOOD FUNCTION IS EXPENSIVE REMEMBER?
59. Metropolis–Hastings in the real world
60. Metropolis–Hastings in the real world D(D+1)/2 tuning parameters that's
61. Metropolis–Hastings in the real world D(D+1)/2 tuning parameters that's the dimension of your problem
62. Metropolis–Hastings in the real world D(D+1)/2 tuning parameters that's the dimension of your problem Scientific Awesomeness how hard is MCMC Metropolis Hastings how things Should be (~number of parameters)
63. that being said, go code up your own Metropolis–Hastings code right now. (well not right right now)
64. Jonathan Goodman Jonathan Weare (dfm.io/mcmc-gw10) "Ensemble samplers with affine invariance"
65. Ensemble Samplers in the real world
66. Ensemble Samplers in the real world
67. Ensemble samplers with affine invariance
68. y A x + b Affine Transformation hard easy
69. y A x + b Affine Transformation hard easy THE SAME
70. Ensemble Samplers in the real world
71. Ensemble Samplers in the real world
72. choose a helper Ensemble Samplers in the real world
73. choose a helper Ensemble Samplers in the real world
74. choose a helper Ensemble Samplers in the real world propose a new position
75. choose a helper Ensemble Samplers in the real world accept? p(accept) = min ✓ 1, ZD 1 p( x ) p( x 0) ◆ propose a new position
76. Ensemble Samplers in the real world
77. Ensemble Samplers in the real world
78. choose a helper Ensemble Samplers in the real world
79. Ensemble Samplers in the real world
80. go code up your own Ensemble Sampler code right now?
81. emceethe MCMC Hammer introducing arxiv.org/abs/1202.3665 dan.iel.fm/emcee it's hammer time!
82. emceethe MCMC Hammer introducing arxiv.org/abs/1202.3665 dan.iel.fm/emcee it's hammer time! pip install emcee to install:
83. it's hammer time! using emcee is easy let me show you...
84. N Y n=1 1 p 2 ⇡ ( 2 n + 2 ) exp  1 2 [yn (m xn + b)] 2 2 n + 2 p ({ yn, xn, 2 n } | m, b, 2) =
85. import numpy as np import emcee def lnprobfn(p, x, y, yerr): m, b, d2 = p if not 0 <= d2 <= 1: return -np.inf ivar = 1.0 / (yerr ** 2 + d2) chi2 = np.sum((y - m * x - b) ** 2 * ivar) return -0.5 * (chi2 - np.sum(np.log(ivar))) # Load data. # x, y, yerr = ... # Set up sampler and initialize the walkers. nwalkers, ndim = 100, 3 sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprobfn, args=(x, y, yerr)) p0 = [[np.random.rand(), 10 * np.random.rand(), np.random.rand()] for k in range(nwalkers)] # Go. sampler.run_mcmc(p0, 1000) m, b, d2 = sampler.flatchain.T
86. import numpy as np def lnprobfn(p, x, y, yerr): m, b, d2 = p if not 0 <= d2 <= 1: return -np.inf ivar = 1.0 / (yerr ** 2 + d2) chi2 = np.sum((y - m * x - b) ** 2 * ivar) return -0.5 * (chi2 - np.sum(np.log(ivar))) p ( 2 ) = ⇢ 1 , if 0  2  1 0 , otherwise 2 = N X n=1 ( yn m xn + b )2 2 n + 2 ln p ({ xn, yn, n } | m, b, 2) = 1 2 2 1 2 N X n=1 ln 2 n + 2
87. import emcee nwalkers, ndim = 100, 3 sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprobfn, args=(x, y, yerr)) p0 = [[np.random.rand(), 10 * np.random.rand(), np.random.rand()] for k in range(nwalkers)] sampler.run_mcmc(p0, 1000) m, b, d2 = sampler.flatchain.T mk ⇠ U(0, 1) bk ⇠ U(0, 10) 2 k ⇠ U(0, 1) initialize each walker run MCMC for 1000 steps get the samples initialize the sampler
88. burn-in
91. convergence
92. acceptance fraction?
93. acceptance fraction? whoa!
94. acceptance fraction? whoa! that's a lot of significant digits...
95. dfm/acor 
96. displaying results
99. dfm/triangle.py 
100. sharing results you need a number for your abstract?!?
101. 1 sort the samples 2 compute moments/quantiles (0.16, 0.5, 084)
102. 1 sort the samples 2 compute moments/quantiles (0.16, 0.5, 084) X = ¯ X+ + what does it MEAN?
103. 1 sort the samples 2 compute moments/quantiles (0.16, 0.5, 084) X = ¯ X+ + what does it MEAN? 3 machine readable samplings (including priors values)
104. 1 sort the samples 2 compute moments/quantiles (0.16, 0.5, 084) X = ¯ X+ + what does it MEAN? 3 machine readable samplings (including priors values) my pipe dream
105. Brendon Brewer's words of wisdom...
106. emcee isn't always The Right Choice™ Remember: Brendon Brewer's words of wisdom...
107. what about multimodal densities? ✓ p(✓ | D)
108. what about multimodal densities? ✓ p(✓ | D)
109. what about multimodal densities? ✓ p(✓ | D)
110. what about multimodal densities? ✓ p(✓ | D) p(accept)
111. what about multimodal densities? ✓ p(✓ | D) p(accept)
112. what if we want to compare models? p(✓ | D, M) = 1 Z(M) p(✓ | M) p(D | ✓, M) Z = p(D | M) = Z p(D, ✓ | M) d✓
113. what if we want to compare models? p(✓ | D, M) = 1 Z(M) p(✓ | M) p(D | ✓, M) Z = p(D | M) = Z p(D, ✓ | M) d✓
114. what to do?
115. what to do? ✓ p(✓ | D) 1 burn-in & priors
116. what to do? ✓ p(✓ | D) 1 burn-in & priors Z = p(D | M) = Z p(D, ✓ | M) d✓ 2 use a different algorithm (gasp!)
117. what to do? ✓ p(✓ | D) 1 burn-in & priors Z = p(D | M) = Z p(D, ✓ | M) d✓ 2 use a different algorithm (gasp!) github.com/eggplantbren/DNest3
118. sometimes it is useful
119. emcee is community supported.
120. emcee is community supported.
121. emcee has good documentation.
122. emcee has a live support team. danfm@nyu.edu
123. Alex Conley (UC Boulder) Jason Davies (Jason Davies Ltd.) Will Meierjurgen Farr (Northwestern) David W. Hogg (NYU) Dustin Lang (CMU) Phil Marshall (Oxford) Ilya Pashchenko (ASC LPI, Moscow) Adrian Price-Whelan (Columbia) Jeremy Sanders (Cambridge) Joe Zuntz (Oxford) Eric Agol (UW) Jo Bovy (IAS) Jacqueline Chen (MIT) John Gizis (Delaware) Jonathan Goodman (NYU) Marius Millea (UC Davis) Jennifer Piscionere (Vanderbilt) contributors
124. take home. p(D | ✓, ↵) your model:
125. take home. p(D | ✓, ↵) it's hammer time! p(D | ✓) / Z p(↵) p(D | ✓, ↵) d↵ your model:
126. take home. p(D | ✓, ↵) ✓ = it's hammer time! p(D | ✓) / Z p(↵) p(D | ✓, ↵) d↵ your model:
127. Alex Conley (UC Boulder) Jason Davies (Jason Davies Ltd.) Will Meierjurgen Farr (Northwestern) David W. Hogg (NYU) Dustin Lang (CMU) Phil Marshall (Oxford) Ilya Pashchenko (ASC LPI, Moscow) Adrian Price-Whelan (Columbia) Jeremy Sanders (Cambridge) Joe Zuntz (Oxford) Eric Agol (UW) Jo Bovy (IAS) Jacqueline Chen (MIT) John Gizis (Delaware) Jonathan Goodman (NYU) Marius Millea (UC Davis) Jennifer Piscionere (Vanderbilt) thanks!