Shulou(Shulou.com)06/01 Report--

This article mainly introduces "how to use Python to achieve MCMC model". In daily operation, I believe many people have doubts about how to use Python to achieve MCMC model. The editor consulted all kinds of materials and sorted out simple and easy-to-use operation methods. I hope it will be helpful for you to answer the doubts about "how to use Python to achieve MCMC model". Next, please follow the editor to study!

Introduction

Real-life data is never perfect, but we can still extract valuable information from noise data through the right model.

Sleep data

On the graph, each data node is marked as a point, and the strength of the point indicates the number of observations at a specified time. My watch only records the time it takes to fall asleep, so in order to expand the amount of data, I add nodes to each minute on both sides of the time. If my watch records that I fell asleep at 10:05 in the evening, every minute before is represented as 0 (awake), and every minute after that is represented as 1 (falling asleep). This will expand observations for about 60 nights to 11340 data nodes.

As you can see, I always fall asleep after 10:00 at night, and we want to create a model that captures the transition from wake to sleep probability. When the model changes from waking (0) to falling asleep (1) at a precise time, we can use a simple ladder function for the model, because I don't fall asleep at the same time every night. so we need a function to model the transformation process as a gradual process. The best choice for a given data is a logical function that converts smoothly between 0 and 1. The following is the logical equation of the probability of falling asleep as a function of time.

β (beta) and α (alpha) are the model parameters that we must learn in the MCMC process. Logical functions with different parameters are shown below:

Logical functions are suitable for these data because the probability of sleep transition is gradual, and the changes in my sleep samples are obtained at the same time. We want to be able to add a time variable t to the function and find out the probability of falling asleep, which must be between 0 and 1. In order to create this model, we use a classification technique as MCMC and use the data to find the best α and β parameters.

Markov chain Monte Carlo (MCMC)

MCMC is a kind of method of sampling from probability distribution to construct the most likely distribution. We cannot calculate the logical distribution directly, so we generate thousands of values called samples for the parameters of the function (α and β) to create an approximation of the distribution. The idea behind MCMC is that as more samples are generated, the approximations will get closer and closer to the actual distribution.

The MCMC method has two parts. Monte Carlo is a general technique that uses repeated random samples to obtain numerical results. Monte Carlo can be thought of as conducting many experiments, each time changing the variables in the model and observing the reaction. By selecting random values, we can study a large part of the parameter space, the range of possible variable values. Parameter space, as shown in the following figure:

Obviously, we can't try every point in the graph, but by randomly sampling from areas with a higher probability (red), we can create the most likely model.

Markov chain (Markov Chain)

Markov chain is a process in which the current state is dependent on the next state. Markov chains do not record the state, because only the current state is important, and it does not matter how it reaches that state. If this is a little difficult to understand, then consider a daily phenomenon, the weather. If we want to predict tomorrow's weather, we can only make a reasonable forecast based on today's weather. If it snows today, let's look at the historical data of the weather distribution the next day after it snows to predict the probability of what the weather will be like tomorrow. The concept of Markov chain is that we do not need to know the entire historical data of one process to predict the next output, and similar phenomena can be well predicted in many real situations.

Combining the ideas of Markov chain and Monte Carlo, MCMC is a method of repeatedly extracting random values based on the parameters of the current value pair distribution. The samples of each value are random, but the selection of these values is limited by the assumed prior distribution of the current state and parameters. MCMC can be regarded as a random walk that gradually converges to the actual distribution.

In order to extract the random values of α and β, we need to assume the prior distribution of these values. Since we do not presuppose the parameters in advance, we can use a normal distribution. The normal distribution or Gaussian distribution is defined by the average, which indicates the location, variance and distribution range of the data. The following figure shows several normal distributions with different averages and ranges:

The specific MCMC algorithm we use is called Metropolis Hastings. In order to associate the observed data with the model, each time a set of random values is drawn, the algorithm is evaluated according to the data. If these random values are inconsistent with the data, they are rejected and the model remains in the current state. Conversely, if consistent with the data, these random values are assigned to the parameter and become the current state. This process lasts a certain number of iterations, so the accuracy of the model will improve with the increase of the number of iterations.

To sum up, the basic steps to solve the problem with MCMC are as follows:

(1) Select α and β and the parameter initial value set of logic function.

(2) randomly assign new values to α and β according to the current state.

(3) check whether the new random value is consistent with the observed value. If it does not match, reject these random values and return to the previous state, otherwise, if they match, accept and update to the new current state

(4) repeat steps 2 and 3 if iteration is required.

The algorithm returns the values of all alpha and beta generated by it. We can then use the average of these values as the most likely final values of α and β in the logical function. MCMC cannot return a "True" value, but rather an approximate value of a distribution. The final model of the sleep probability obtained from the existing data is a logical function with α and β averages.

The realization of Python

To implement MCMC in Python, we will use the Bayesian inference library PyMC3, which abstracts most of the details, allowing us to create models instead of empty theories.

The following code creates a complete model with parameters α and β, probability p, and observations. The step variable refers to a specific algorithm, and the variable sleep_trace holds all the parameter values generated by the model.

Withpm.Model () as sleep_model:# Create the alpha and beta parameters# Assume a normal distributionalpha=pm.Normal ('alpha', mu=0.0, tau=0.05, testval=0.0) beta=pm.Normal (' beta', mu=0.0, tau=0.05, testval=0.0) # The sleep probability is modeled as a logistic functionp=pm.Deterministic ('The sleep probability is modeled as a logistic functionp=pm.Deterministic, 1. / (1. + tt.exp (beta * time + alpha) # Create the bernoulli parameter which uses observed data to inform the algorithmobserved=pm.Bernoulli (' obs', p Observed=sleep_obs) # Using Metropolis Hastings Samplingstep=pm.Metropolis () # Draw the specified number of samplessleep_trace=pm.sample (N_SAMPLES, step=step)

To better understand the operation of the code, we can take a look at the alpha and beta values generated during the operation of all the models.

The above picture is called a trajectory map. We can see that each state is related to the previous Markov chain, but its value oscillates obviously in the Monte Carlo sampling.

In MCMC, it is common to discard up to 90 per cent of the tracks. The algorithm does not immediately converge to the real distribution, and the initial value is often not accurate. Later parameter values are usually better, which means that they should be used to build the model. We used 10000 samples and discarded the previous 50%, but hundreds or millions of samples may be used in enterprise applications.

MCMC converges to true values, but it can be difficult to evaluate convergence. This is an important factor if we want the most accurate results. The PyMC3 library has created functions to evaluate the quality of the model, including trajectory maps and autocorrelation maps.

Pm.traceplot (sleep_trace, ['alpha',' beta']) pm.autocorrplot (sleep_trace, ['alpha',' beta'])

Trajectory map (top) and autocorrelation graph (bottom)

Sleep model

After finally building and running the model, it is time to use it. We take the average of the last 5000 α and β samples as the most likely parameter values, which allows us to create a single graph to simulate the probability of falling asleep after a specified time:

The model can well represent the data, but also capture the inherent changes in my sleep pattern, it does not give a result, but a probability. For example, you can query the model to find out how likely I am to fall asleep at a specified time, and to find a point in time with a probability of more than 50%:

The probability of falling asleep at 09:30 at night: 4.80%.

The probability of falling asleep at 09:30 at night: 27.44%.

The probability of falling asleep at 10:00 at night: 73.91%.

At 10:14 at night, the probability of falling asleep increased to more than 50 per cent.

You can see that the average time I fall asleep is around 10:14 in the evening.

These values are the most likely estimates for a given data. However, there will be these probability-related uncertainties because the model is approximate. To express this uncertainty, we can use all the samples of α and β to predict the probability of falling asleep at a given point in time, and then draw a column based on the results:

The above results give a better index that can be achieved by the MCMC model. This method does not find a correct answer, but a sample of possible values. Bayesian reasoning is of practical use because it represents the predicted results in the form of probability. We can say that we have the most likely answer, but it is more accurate to say that any prediction is a range of values.

Wake-up model

I can use the data I woke up in the morning to find a similar model. I usually get up at 6: 00 in the morning, and the following picture shows the final model from sleeping to waking up.

The model can be used to find the probability of falling asleep at a given time and the time when I am most likely to wake up.

The probability of waking up at 05:30: 14.10%.

The probability of waking up at 6: 00 in the morning: 37.94%.

The probability of waking up at 06:30: 69.49%.

The chance of waking up at 06:11 is more than 50%.

Duration of sleep

The last thing I want to create is a sleep time model. First, we need to find a function to simulate the distribution of the data, but only by examining the data.

A normal distribution can be achieved, but it does not get the outliers on the right. You can use two independent normal distributions to represent these two patterns, but I will use a skewed distribution. Skewness distribution has three parameters, average, variance and deviation, and these three parameters must be learned from the MCMC algorithm. The following code creates the model and implements Metropolis Hastings sampling:

Withpm.Model () asduration_model:# Three parameters to samplealpha_skew=pm.Normal ('alpha_skew', mu=0, tau=0.5, testval=3.0) mu_ = pm.Normal (' mu', mu=0, tau=0.5, testval=7.4) tau_ = pm.Normal ('tau', mu=0, tau=0.5, testval=1.0) # Duration is a deterministic variableduration_ = pm.SkewNormal (' duration', alpha=alpha_skew, mu= mu_, sd=1/tau_, observed= duration) # Metropolis Hastings for samplingstep=pm.Metropolis () duration_trace=pm.sample (N_SAMPLES, step=step)

Now we can use the average of the three parameters to construct the most likely distribution. The following figure shows the final skewed distribution of the data:

The above figure seems to be very realistic, and through the query model, we can find that I get at least a certain duration of sleep, and the probability of the most likely range of falling asleep time:

The probability of falling asleep for at least 6.5 hours is 99.16%.

The probability of falling asleep for at least 8 hours is 44.53%.

The probability of falling asleep for at least 9 hours is 10.94%.

The most likely duration of falling asleep is 7.67 hours.

At this point, the study on "how to use Python to implement the MCMC model" is over. I hope to be able to solve your doubts. The collocation of theory and practice can better help you learn, go and try it! If you want to continue to learn more related knowledge, please continue to follow the website, the editor will continue to work hard to bring you more practical articles!

Welcome to subscribe "Shulou Technology Information " to get latest news, interesting things and hot topics in the IT industry, and controls the hottest and latest Internet news, technology news and IT industry trends.