matlab mcmc tutorial

【实例简介】
文档详细讲述了mcmc算法的matlab实现教程，方便理工及经管类的朋友进行蒙特卡洛马尔科夫算法分析时取用。
CONTENTS 4.3 Graphical Model notation 54 4.3.1 Example: Consensus Modeling with gaussian variables 54 5 Approximate Inference in Graphical Models 57 5. 1 Prior predictive distributions 57 5.2 Posterior distributions 5.2.1 Rejection Sampling 5.2.2 CMC Sampling 61 5.2.3 Example: Posterior Inference for the consensus model with normally distributed variables 5.2.4 Example: Posterior Inference for the consensus model with contaminants 66 6 Sequential Monte Carlo 69 6.1 Hidden markov models 70 6. 1.1 Example HMM with discrete outcomes and states 6. 1.2 Viterbi Algorithm 6.2 Bayesian Filtering 73 6.3 Particle filters 74 6.3.1 Sampling Importance Resampling(SIR) 74 6.3.2 Direct Simulation 76 Note to students Eⅹ excises This course book contains a number of exercises in which you are asked to simulate Matlab code, produce new code, as well as produce graphical illustrations and answers to questions The exercises marked with *x are optional exercises that can be skipped when time is limited Organizing answers to exercises It is helpful to maintain a document that organizes all the material related to the exercises Matlab can facilitate part of this organization using the "publish"option. For example, if you have a matlab script that produces a figure you can publish the code as well as the figure produced by the code to a single external document. You can find the publishing option in the Matlab editor under the file menu. You can also use the publish function directly in the command window. You can change the publish configuration (look under the file menu of the editor window) to produce pdfs, Word documents and a number of file formats For this course however, you can hand in your answers in any way you see fit and sometimes, it might be useful to just hand in handwritten answers. I Matlab documentation It will probably happen many times that you will need to find the name of a matlab function or a description of the input and output variables for a given Matlab function. It is strongly recommended to have the matlab documentation running in a separate window for quick consultation. You can access the Matlab documentation by typing doc in the command window For specific help on a given matlab function, such as the function fprintf you can type doc fprint f to get a help screen in the matlab documentation window or help fprintf to get a description in the matlab command window note: in this version, a new chapter 3 was inserted which is not fully completed yet Chapter 1 Sampling from Random variables Probabilistic models proposed by researchers are often too complicated for analytic ap proaches. Increasingly, researchers rely on computational, numerical-based methods when dealing with complex probabilistic models. By using a computational approach, the re- searcher is freed from making unrealistic assumptions required for some analytic techniques (e. g. such as normality and independence The key to most approximation approaches is the ability to sample from distributions Sampling is needed to predict how a particular model will behave under some set of cir cumstances, and to find appropriate values for the latent variables("parameters?")when applying models to experimental data. Most computational sampling approaches turn the problem of sampling from complex distributions into subproblems involving simpler sampling distributions. In this chapter, we will illustrate two sampling approaches: the inverse trans formation method and rejection sampling. These approaches are appropriate mostly for the univariate case where we are dealing with single-valued outcomes. In the next chapter, we discuss markov chain Monte Carlo approaches that can operate efficiently with multivariate distributions 1.1 Standard distributions Some distributions are used so often, that they become part of a standard set of distributions supported by matlab. The Matlab Statistics Toolbox supports a large number of probabil- ity distributions. Using Matlab, it becomes quite easy to calculate the probability density cumulative density of these distributions, and to sample random values from these distribu- tions. Table 1. 1 lists some of the standard distributions supported by matlab. The matlab documentation lists many more distributions that can be simulated with Matlab. Using online resources, it is often easy to find support for a number of other common distributions To illustrate how we can use some of these functions, Listing 1. I shows Matlab code that visualizes the normal(u, o)distribution where ul=100 and o= 15. To make things concrete imagine that this distribution represents the observed variability of IQ coefficients in some population. The code shows how to display the probability density and the cumulative 4 CHAPTER 1. SAMPLING FROM RANDOM VARIABLES Table 1.1: Examples of Matlab functions for evaluating probability density, cumulative den- sity and drawing random numbers Distribution PDF CDE Random Number generation Normal normpdf normcdf norm Uniform(continuous) unifpdf unifcdf unifrnd Beta betapdf betacdf betarnd Exponential exppdf expcdf expand Uniform(discrete unidpdf unidedf unidrnd Binomial binopdf binocdf binornd Multinomial mnpdf morn Poisson poisspdf I poisscdf poissrnd Probability Density Function Cumulative Dersity Function Histogram of random value 巨0015 0 0005 Figure 1.1: Illustration of the Normal(u, a) distribution where u=100 and o-15 density. It also shows how to draw random values from this distribution and how to visualize the distribution of these random samples using the hist function. The code produces the output as shown in Figure 1.1. Similarly, Figure 1.2 visualizes the discrete distribution Binomial(n, e) distribution where N= 10 and 9=0.7. The binomial arises in situations where a researcher counts the number of successes out of a given number of trials. For example, the Binomial(10, 0.7) distribution represents a situation where we have 10 total trials and the probability of success at each trial, 6, equals 0.7 Exercises 1. Adapt the Matlab program in Listing 1. 1 to illustrate the Betala, B)distribution where &=2 and B=3. Similarly, show the Exponential(A)distribution where A=2 2. Adapt the matlab program above to illustrate the Binomial(N, 0) distribution where N=10 and 8=0.7. Produce an illustration that looks similar to Figure 1.2 3. Write a demonstration program to sample 10 values from a Bernoulli(0)distribution with 0=0.3. Note that the Bernoulli distribution is one of the simplest discrete distri CHAPTER 1. SAMPLING FROM RANDOM VARIABLES Listing 1.1 Matlab code to visualize Normal distribution 1 Explore the Normal distribution N( mu, sigma 100;是 the mear 15; the standard deviation 4 xm- 70; minimum x value for pdf and cdf plot 5 xmaX 130; maximum x value for pdf and caf plot 6 I 100; number of points on pdf and cdf plot 7 k 10000;9 number of random draws for histogram 9 create a set of values ranging from xmin to xmax linspace( xmin, max n )i 11 p= normpdf( x, mu, sigma )i calculate the pdf 12 C= normcdf(x, mu, sigma )i calculate the cdf figure(1 )icti create a new figure and clear the contents 16subp1ot(1,3,1); 17 plot(x,p 18 label('x)i ylabel( 'pdf) 19 title( 'Probabi-ity Density Function 21 subplot( 1,3,2) 22 plot( x 23 xlabel(x)i ylabel( ' cdf) 24 tit.le( Cumul at:ve Density Function 26 draw k random numbers from a N( mu, sigma ) distribution normrnd( mu, sigma ,k, 1 29 subplot( 1,3,3 ) 30 hist(y, 20 31 Label( x)yl abel(' frequency 32 title( 'histogram of random values Probability Distribution Cumulative Probability Distribution Histogram 0.35 0. 806 0.15 0.05 0 012345678910 012345678910 0123456 8910 Figure 1.2 Illustration of the Binomial(N, 0) distribution where N=10 and 0=0.7 CHAPTER 1. SAMPLING FROM RANDOM VARIABLES butions to simulate. There are only two possible outcomes, 0 and 1. With probability 6, the outcome is 1, and with probability 1-0, the outcome is 0. In other words P(X=1)=0, and p(X=0)=1-8. This distribution can be used to simulate outcomes in a number of situations, such as head or tail outcomes from a weighted coin, correct/incorrect outcomes from true/false questions, etc. In Matlab, you can simulate the bernoulli distribution using the binomial distribution with N=1. How- ever, for the purpose of this exercise, please write the code needed to sample bernoull distributed values that does not make use of the built-in binomial distribution 4. It is often useful in simulations to ensure that each replication of the simulation gives the exact same result. In Matlab. when drawing random values from distributions the values are different every time you restart the code. There is a simple way to seedthe random number generators to insure that they produce the same sequence Write a Mat lab script that samples two sets of 10 randon values drawn from a, uniform distribution between 0, 1. Use the seeding function between the two sampling steps to demonstrate that the two sets of random values are identical. Your matlab code could use the following line seed=1: rand(state,, seed); randn(state, seed) Note that there are more advanced methods to seed the random number generator in Matlab using RandStream methods. For most applications in this course, such Imet,hods are not, needed 5. Suppose we know from previous research that in a given population, IQ coefficients are normally distributed with a mean of 100 and a standard deviation of 15. Calculate the probability that a randomly drawn person from this population has an iQ greater than 110 but smaller than 130. You can achieve this using one line of matlab code What does this look like? ck 6 The Dirichlet distribution is currently not supported by Matlab. Can you find a matlab function, using online resources, that implements the sampling from a dirichlet distribution? 1.2 Sampling from non-standard distributions Suppose we wish to sample fron a distribution that is not one of the standard distributions that is supported by Matlab. In modeling situations, this situation frequently arises, be cause a researcher can propose new noise processes or combinations of existing distributions Computational methods for solving complex sampling problems often rely on sampling dis tributions that we do know how to sample from efficiently. The random values from these simple distributions can then be transformed or compared to the target distribution. In fact some of the techniques discussed in this section are used by matlab internally to sample from distributions such as the Normal and Exponential distributions CHAPTER 1. SAMPLING FROM RANDOM VARIABLES 1.2.1 Inverse transform sampling with discrete variables Inverse transform sampling(also known as the inverse transform method is a method for generating random numbers from any probability distribution given the inverse of its cumu lative distribution function. The idea is to sanple uniformly distributed random nunbers (between 0 and 1)and then transform these values using the inverse cumulative distribution function. The simplicity of this procedure lies in the fact that the underlying sampling is just based on transformed uniform deviates. This procedure can be used to sample many different kinds of distributions. In fact, this is how matlab implements many of its random number generators It is easiest to illustrate this approach on a discrete distribution where we know the probability of each individual outcome. In this case, the inverse transform method just requires a simple table lookup To give an example of a non-standard discrete distribution, we use some data from ex periments that have looked at how well humans can produce uniform random numbers(e. g Treisman and Faulkner, 1987). In these experiments, subjects produce a large number of ran- dom digits(0,, 9) and investigators tabulate the relative frequencies of each random digit produced. As you inight suspect, subjects do not always produce uniform distributions Table 1.2. 1 shows some typical data. Some of the low and the high numbers are underrepre sented while some specific digits(e. g. 4)are overrepresented. For some reason, the digits 0 and 9 were never generated by the subject(perhaps because the subject misinterpreted the instructions). In any case, this data is fairly typical and demonstrates that humans are not very good are producing uniformly distributed random numbers Table 1.2: Probability of digits observed in human random digit generation experiment. The generated digit is represented by X; p(X)and F(X)are the probability mass and cumulative probabilities respectively. The data was estimated from subject 6, session 1, in experiment by Treisman and Faulkner(1987) Xp()F(X) 00.000.00 10.1000.100 30.0950.285 40.2000.485 0.1750.660 6|0.1900850 70.0500.900 80.1001.000 90.0001.000 Suppose we now want to mimic this process and write an algorithm that samples digits CHAPTER 1. SAMPLING FROM RANDOM VARIABLES according to the probabilities shown in Table 1.2.1. Therefore, the program should produce a 4 with probabilit, y. 2, a 5 with probability. 175, etc. For example, the code in Listing 1.2 implements this process using the built-in matlab function randsample. The code produces the illustration shown in Figure 1.2.1 Instead of using the built-in functions such as randsample or mnrnd, it is helpful to consider how to implement the underlying sampling algorithn using the inverse transform method. We first need to calculate the cumulative probability distribution. In other words we need to know the probability that we observe an outcome equal to or smaller than some particular value. If F(X) represents the cumulative function, we need to calculate F(X X ) For discrete distributions, this can be done using simple summation The cumulative probabilities of our example are shown in the right column of Table 1.2.1 In the inverse transform algorithm, the idea is to sample uniform random deviates (i.e. random numbers between 0 and 1)and to compare each random number against the table of cumulative probabilities. The first outcome for which the random deviate is smaller than (or is equal to) the associated cumulative probability corresponds to the sampled outcome Figure 1. 2.1 shows an example with a uniform random deviate of U=0.8 that leads to a sampled outcome X=6. This process of repeated sampling of uniform deviates and comparing these to the cumulative distribution forms the basis for the inverse transform method for discrete variables. Note that we are applying an inverse function, because we are doing an inverse table lookup Listing 1.2: Matlab code to simulate sampling of random digits 2 Simulate the distribution observed in the 2 human random digit generation task 4 probab=lities for each digit 5 theta=[0.000 号 digit0 0.10C ,各 digit1 0.090;..号 digit2 0.095;...号 digit3 0.200;,,.爷 digit4 0.175;... digit5 0.190 号 digit6 0.050;,,号 digit7 0.100 o digit 8 0.000]..号icit9 15 16 fix the random nurrber generator 17 seed =l; rand( 'state, seed )i 19 et's say we draw k random values K=10000 21 digiset =0: 9i 22 Y = randsample(digitset, k, true, theta)i 【实例截图】
【核心代码】