home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1993 #1 / NN_1993_1.iso / spool / sci / math / stat / 2735 < prev    next >
Encoding:
Text File  |  1993-01-08  |  7.1 KB  |  185 lines

  1. Newsgroups: sci.math.stat
  2. Path: sparky!uunet!paladin.american.edu!gatech!rpi!zaphod.mps.ohio-state.edu!saimiri.primate.wisc.edu!sdd.hp.com!news.cs.indiana.edu!umn.edu!charlie
  3. From: charlie@umnstat.stat.umn.edu (Charles Geyer)
  4. Subject: Re: Metropolis sampler
  5. Message-ID: <1993Jan8.232157.1897@news2.cis.umn.edu>
  6. Keywords: Metropolis sampler
  7. Sender: news@news2.cis.umn.edu (Usenet News Administration)
  8. Nntp-Posting-Host: isles.stat.umn.edu
  9. Organization: School of Statistics, University of Minnesota
  10. References: <1993Jan7.120831.11046@cs.kuleuven.ac.be>
  11. Date: Fri, 8 Jan 1993 23:21:57 GMT
  12. Lines: 171
  13.  
  14. In article <1993Jan7.120831.11046@cs.kuleuven.ac.be> maurits@cs.kuleuven.ac.be
  15. (Maurits Malfait) writes:
  16.  
  17. >We are working on an algorithm  that  incorporates  a  Metropolis
  18. >sampler  (*).   Therefore  we  are  interested  in references and
  19. >information concerning:
  20. >
  21. >1. The convergence properties of a Metropolis sampler:
  22. >   - when  is  the  starting  state  of  a   Metropolis   sampler
  23. >     'forgotten'?
  24. >   - how fast and how  good  does  the  sampler  approximate  the
  25. >     Boltzmann distribution?
  26.  
  27. IMHO that is the wrong question to ask.  The right question is how fast
  28. does the sampler mix.  This is discussed in the current issue of Statistical
  29. Science in papers by Gelman and Rubin and by myself with comments by about
  30. 10 discussants.  It should be out real soon now.
  31.  
  32. The output of the sampler is a Markov chain X_1, X_2, ... whose equilibrium
  33. distribution is the distribution of interest, the one you want to sample
  34. from.  Averages with respect to that distribution are approximated by
  35. averages over the chain.
  36.  
  37.                                          1  n
  38.   E g(X)   is approximated by   gbar_n = - sum g(X_i)
  39.                                          n i=1
  40.  
  41. the asymptotic distribution of gbar_n is usually asymptotically normal
  42.  
  43.   sqrt(n) [gbar_n - E g(X)]   converges in distribution to   N(0,sigma^2)
  44.  
  45. where
  46.  
  47.                         infinity
  48.   sigma^2 = gamma_0 + 2   sum    gamma_i                (*)
  49.                           i=1
  50.  
  51. where
  52.  
  53.   gamma_i = Cov[ g(X_0), g(X_i) ]
  54.  
  55. at stationarity, i. e. if X_0 has the equilibrium distribution and X_i
  56. is the i-th step of a chain started at X_0.
  57.  
  58. More precisely, if the sampler is reversible and irreducible, the CLT holds
  59. if the sum (*) is finite (Kipnis and Varadhan, Commun. Math. Phys., 1986)
  60. and in all nice cases (the precise technical condition is Harris recurrence),
  61. the asymptotics do not depend on the starting distribution.
  62.  
  63. Hence asymptotically the starting distribution is irrelevant.
  64.  
  65. In practice one does usually want to throw away a short initial segment of
  66. a run, but there is no reason to run until the initial state is "forgotten"
  67. or the chain "reaches equilibrium".  Strictly speaking, it never does. All
  68. that is necessary is that the chain reach a state that is reasonably likely,
  69. a state that has reasonable probability of occuring in a sample of the
  70. size collected.  If one knows the mode of the distribution, for example,
  71. one could start there with no burn-in.
  72.  
  73. The trick is to prove that (*) is finite and then to estimate it accurately.
  74. This is still an area of active research, but I do give some advice in the
  75. Statistical Science paper mentioned above.
  76.  
  77. >2. Variants of the Metropolis algorithm and recent developments.
  78.  
  79. The main variant is the so-called Metropolis-Hastings algorithm (Hastings,
  80. Biometrika, 1970), which I described here last year
  81.  
  82. > Newsgroups: sci.math.stat
  83. > From: charlie@umnstat.stat.umn.edu (Charles Geyer)
  84. > Subject: Re: Sampling from mixtures of distributions
  85. > Organization: School of Statistics, University of Minnesota
  86. > Date: Thu, 23 Apr 1992 00:07:22 GMT
  87. > In article <20514@castle.ed.ac.uk> efsv02@castle.ed.ac.uk (M J Brewer) writes:
  88. > >I am conducting Gibbs sampling on graphical models.
  89. > >
  90. > >I need to be able to sample from, for example,
  91. > >
  92. > >  f( x )  =  k * exp(-.5*((x-u)/v)^2) * exp(-.5*((g(x)-w)/y)^2)
  93. > >
  94. > >where k is an unknown constant ( => rejection sampling)
  95. > >      u,w are known "means",
  96. > >      v,y are known "variances",  (cf. Normal)
  97. > >      and g(x) is a NON-LINEAR function of x.
  98. > >
  99. > >I'm not sure about the envelope for rejection sampling; I've plotted the
  100. > >above function, and for a quadratic g(), a Normal envelope would seem ok,
  101. > >but I'm not sure about the choice of variance; also I need to be sure
  102. > >that the envelope "blankets" f() in the tails.
  103. > >
  104. > >Does anyone have any ideas?
  105. > Yes.  If the one-dimensional conditionals needed for the Gibbs sampler
  106. > are hard to sample from, DON'T USE the Gibbs sampler.
  107. > There is sure to be a much simpler Metropolis-Hastings algorithm that
  108. > not only will be easier to code but also will run faster.
  109. > For those who may still think that the Gibbs sampler is the name of the
  110. > game, the Metropolis-Hastings algorithm (Hastings, Biometrika, 1970) is
  111. > the following.
  112. > One update step (say just updating one variable) moving from the state x
  113. > to a new state x'
  114. >   1. simulate a realization y from the density q(x, )
  115. >      [the density q is chosen to be easy to simulate from but is
  116. >      otherwise arbitrary]
  117. >   2. calculate the Hastings ratio
  118. >              p(y) q(y,x)
  119. >          r = -----------
  120. >              p(x) q(x,y)
  121. >      [where p is the density of the distribution of interest times
  122. >      an unknown proportionality constant which cancels]
  123. >   3. (Metropolis rejection)
  124. >          if (r >= 1) then
  125. >            x' := y
  126. >          else
  127. >            generate u uniform(0,1)
  128. >            if u < r then
  129. >              x' := y
  130. >            else
  131. >              x' := x
  132. >            fi
  133. >          fi
  134. > The Gibbs sampler is the special case where q is the one-dimensional
  135. > conditional.  It avoids the rejection step (r is always exactly 1),
  136. > but it does so at the cost of complicating step 1 and so is rarely
  137. > worth the trouble unless the sampling from one-dimensional conditionals
  138. > is very fast.
  139. > Almost any q will do (one does have to check that the resulting Markov
  140. > chain is irreducible).  Choose one that makes the update fast.  The best
  141. > way to speed up your sampler is to speed up the update even if that increases
  142. > the sampling variance a bit.  Though more samples are required for the
  143. > same accuracy, they are generated much faster.
  144.  
  145. The Metropolis algorithm requires q(x,y) = q(y,x) for no particular reason.
  146. The more general Hastings update permits a wider range of "proposal
  147. distributions".
  148.  
  149. BTW, note that there is no need for the distribution being sampled to
  150. be a Boltzman distribution (what statisticians call an exponential family
  151. model).
  152.  
  153. >3. Possibilities to accelerate a Metropolis sampler.
  154.  
  155. There are many.  One is the Swendsen-Wang algorithm and related schemes
  156. (Swendsen and Wang, Phys. Rev. Lett., 1987; Wang and Swendsen, Physica A,
  157. 1990).  Another is the Metropolis-coupled chain scheme mentioned by Geyer
  158. (Computing Science and Statistics: Proc. 23rd Symp. Interface, 1991).
  159. These are, however, special cases of the general Metropolis-Hastings
  160. algorithm.
  161.  
  162. The main point is that the Metropolis-Hastings algorithm gives tremendous
  163. scope for experimentation.  There are lots of things to try.
  164.  
  165. -- 
  166. Charles Geyer
  167. School of Statistics
  168. University of Minnesota
  169. charlie@umnstat.stat.umn.edu
  170.