Epidemic disease prediction model constructed based on pedigree information and multi-chain PMCMC and prediction method thereof

Abstract

The invention discloses an epidemic prediction model constructed based on pedigree information and multi-chain PMCMC and a prediction method thereof, and the modeling method comprises the steps: firstly, selecting an SEIR model as a state transition model for epidemic prediction, and determining and initializing a to-be-estimated model parameter; then representing the number of various epidemic patients by using particles, predicting by using multi-chain particle filtering to obtain the weight of each particle, and weighting through pedigree information; calculating likelihood function values based on the current iteration parameter and the candidate parameter of the next iteration based on the weighted weight, further calculating the acceptance rate of the candidate parameter, determining whether to accept the candidate parameter as the optimal parameter of the next iteration of the model according to the acceptance rate, and repeating the above process until the number of iterations is completed. The state transition model based on the optimal parameters at the moment is the constructed epidemic disease prediction model. The method can predict latent features in the epidemic infection process, and the prediction has better stability and higher accuracy.

Description

Technical Field

[0001] The present invention relates to the application of machine learning technology in the field of computational epidemiology, and specifically to an epidemic prediction model based on pedigree information and multi-chain PMCMC and its prediction method. Background Art

[0002] Traditional epidemic models typically rely on basic frameworks such as SEIR (Susceptible-Exposed-Infectious-Recovered). However, because the infectious characteristics of different epidemics vary, for each outbreak, it is necessary to first determine the transmission parameters in the SEIR model and use this to predict the future development of the epidemic.

[0003] The spread of an epidemic is affected by many complex factors, such as weather, environment, population size, external conditions, and epidemic prevention policies. Therefore, the outbreak and spread of an epidemic is actually a multidimensional, non-Gaussian, and nonlinear dynamic process.

[0004] Particle Markov Chain Monte Carlo (PMCMC) is a powerful statistical inference method that is particularly suitable for dealing with situations with high-dimensional parameters and complex dynamic models. Because the PMCMC method is suitable for nonlinear and non-stationary processes, using PMCMC to predict parameters in epidemic models can achieve higher flexibility and accuracy. This method has been embodied in the Chinese patent application with publication number CN117153423A, "A method for predicting the time of outbreak of emerging infectious diseases based on Bayesian reasoning." However, this patent application still has obvious shortcomings, specifically:

[0005] (1) Convergence and sample accuracy are constrained by model assumptions. When the epidemiological model structure is complex or data is insufficient, especially in the early stages of the epidemic, the limitations of case data significantly increase the uncertainty of the inference results, thereby significantly reducing the accuracy of the prediction results.

[0006] (2) It only uses the particle weights in PMCMC as a reference for the adaptability of the parameters to be measured. However, the trend of the epidemic is extremely unstable. The early data scale is small, and the maximum observation data scale is only the number of new cases in a certain area over 30 days. The amplitude of data changes is insufficient to reflect the overall trend of the epidemic.

[0007] (3) When a large number of new case data are used as observation data, the data may fluctuate drastically (e.g., there are 150 new cases at time t1, which increases sharply to 14,000 at time t2, and then drops again to 320 at time t3), which may lead to large errors in data tracking, thereby affecting the accuracy and efficiency of the prediction.

[0008] (4) The single-chain PMCMC method is used as a method for predicting epidemics. However, the computational cost of the PMCMC method is high. Especially when processing large-scale population data, the single-chain computational efficiency is significantly reduced, which prolongs the duration of the entire prediction process. However, epidemics have high requirements for the timeliness of predictions. It is necessary to implement corresponding epidemic prevention and control, isolation and other policies based on timely prediction results, especially before the next outbreak peak. Therefore, the single-chain PMCMC method is obviously unable to meet the prediction needs. Summary of the Invention

[0009] The present invention provides an epidemic prediction model and prediction method based on pedigree information and multi-chain PMCMC to improve the accuracy of epidemic prediction.

[0010] In order to achieve the above technical objectives, the present invention adopts the following technical solutions:

[0011] A method for constructing an epidemic prediction model based on pedigree information and multi-chain PMCMC, comprising:

[0012] Step 1: Select the SEIR model as the state transition model for epidemic prediction and determine the model parameters to be estimated;

[0013] Step 2: Initialize the model parameters to be estimated and use them as the model parameters θ0 of the g=0th iteration;

[0014] Step 3, calculate the likelihood function value based on the current model parameters;

[0015] Step 3.1, initialize the particle set, the state of each particle includes the number of susceptible, exposed, infected and recovered people;

[0016] Step 3.2, divide the current particle set evenly into D parallel chains;

[0017] Step 3.3: Based on the particle set and state transition model at time t-1, predict the state of each particle in each parallel chain at time t, and then calculate the predicted value of the number of newly infected people at time t;

[0018] Step 3.4: Calculate the weight of each particle in each parallel chain at time t based on the predicted and actual number of newly infected people at time t.

[0019] Step 3.5, calculate the probability of a common ancestor event occurring at time t based on the true number of susceptible and infected individuals in several subintervals at time t; Step 3.5 is performed after or simultaneously with Steps 3.3-3.4;

[0020] Step 3.6, multiply the weight of each particle at time t by the probability of the common ancestor event occurring at time t to obtain the weighted weight of each particle at time t;

[0021] Step 3.7, based on the weighted weights of each particle at time t, calculate the likelihood function value at time t;

[0022] Step 3.8, resample the particles in each parallel chain; wherein, the execution order of step 3.8 and step 3.7 can be interchanged or performed simultaneously;

[0023] Step 3.9: Repeat steps 3.3 to 3.8 to calculate the likelihood function value at time t+1 until the likelihood function value at time T is calculated;

[0024] Step 3.10: Select one particle from the particle set and use its state value from time 1 to time T as the model parameter θ of the current g-th iteration. g The sample status value of

[0025] Step 4: Use the MCMC proposal distribution function to obtain the candidate parameters θ for the g+1th iteration of the model g+1 ', and calculate the candidate parameter θ based on step 3 g+1 'The likelihood function value;

[0026] Step 5: Based on the model parameters θ g The likelihood function value and candidate parameters θ g+1 'Likelihood function value, calculate the candidate parameter θ g+1 'Acceptance rate α; then determine whether to accept the candidate parameter θ based on the acceptance rate α g+1 'As the optimal parameter θ for the g+1th iteration of the model g+1 If not accepted, the model parameter θ of the gth iteration g As the optimal parameter θ for the g+1th iteration of the model g+1 ;

[0027] Step 6: Repeat steps 3 to 5 until the model iterations are completed. The state transition model based on the optimal parameters at this time is the constructed epidemic prediction model.

[0028] Furthermore, the SEIR model is selected as the state transition model for epidemic prediction, which is expressed as:

[0029]

[0030] I t+1 =I t +αE(t)-γI t

[0031] R t+1 =R t +γI t

[0032] Where St 、E t , I t and R t Represent the number of susceptible, exposed, infected, and recovered individuals at time t; A is the population size, and A=S t +E t +I t +R t ; ξ is noise, which is a normal random variable with mean zero and variance 1; F is the constant noise amplitude; β is the transmission coefficient; α is the rate at which individuals enter the infectious period from the latent period; γ is the recovery probability; among them, β = R0 / γ, R0 is the basic reproduction number; γ = 1 / DI, DI represents the average infectious period; α = 1 / DE, DE represents the average latent period.

[0033] Furthermore, the observation error of epidemic prediction follows a normal distribution, so when y is observed at time t t The probability of a case is expressed as:

[0034] P(y t |ΔI t )=N(y t |ρΔI t ,ρτΔI t )

[0035] Where ρ is the observation probability, τ is the scaling factor, and ΔI t For the number of new patients;

[0036] The model parameters to be estimated include: basic reproduction number R0, average incubation period DE, average infectious period DI, observation probability ρ, scaling factor τ, and the length of time d from the first infection time to the actual discovery time of the epidemic virus.

[0037] Furthermore, the predicted value of the number of newly infected people at time t is calculated as follows:

[0038] est_y t =ρ*I t

[0039] Where est_y t represents the predicted value of the number of new epidemic cases at time t obtained according to the above observation equation.

[0040] Furthermore, step 3.4 calculates the weight of each particle in each parallel chain at time t, which is expressed as:

[0041]

[0042] Where, is the weight of particle i at time t, normpdf() is the normal probability density function, real trepresents the true value of the number of new epidemic infections at time t, est_y t It represents the predicted number of new epidemic infections at time t, and sqrt() represents square root operation.

[0043] Furthermore, step 3.5 calculates the probability of a common ancestor event occurring at time t, which is expressed as:

[0044]

[0045] Where W2 represents the probability of a common ancestor event occurring at time t, k and m represent the number and index of common ancestor events occurring at time t, respectively. m represents the time when the mth common ancestor event occurs within time t, dt represents the time intervals divided within time t, and λ dt represents the common ancestor rate in the time interval dt under the SEIR model, I dt S represents the number of infected people predicted by the SEIR model within the time interval dt, dt It represents the number of susceptible individuals predicted by the SEIR model within the time interval dt.

[0046] Furthermore, step 3.6 performs resampling in parallel through multiple chains. The higher the weight of a particle, the higher its probability of accepting resampling. After resampling, new particles with equal weights are obtained and the number is the same as before resampling. After resampling, a particle exchange mechanism is introduced between different chains, that is, particles in different chains are randomly exchanged.

[0047] Furthermore, the calculation formula for the likelihood function value in step 3.7 is:

[0048]

[0049] Where likelihood represents the likelihood function, ∑W s (X 1:N ) represents the sum of the weighted weights of all N particles in the particle set at time s.

[0050] The present invention also provides an epidemic prediction method, comprising:

[0051] Build an epidemic prediction model using any of the above methods;

[0052] Obtain the true value of the number of susceptible, exposed, infected, and recovered individuals at a certain time t0; the time t0 is the moment closest to the current time and the number of susceptible, exposed, infected, and recovered individuals is known;

[0053] According to the actual values of the number of susceptible people, exposed people, infected people and recovered people at time t0, the constructed epidemic prediction model is used to make a rolling prediction of the number of susceptible people, exposed people, infected people and recovered people at the time to be predicted.

[0054] Beneficial effects

[0055] The present invention uses future new patients as hidden variables. PMCMC uses particles to sample using the SEIR model as the state transition equation to predict the particle state at the next moment, thereby obtaining the number of newly infected people at the next moment. The particle sampling process fully reflects the characteristics of the epidemic from the initial outbreak, incubation period, infectious period, and recovery period, and uses pedigree information to infer the weight of the predicted state at the current moment. Because this process relies on the SEIR model, the weight also fully reflects the characteristics of the epidemic from the initial outbreak, incubation period, infectious period, and recovery period. Therefore, the sampling process of the present invention conforms to the latent characteristics of the epidemic transmission process and has better stability and higher accuracy than existing methods.

[0056] Existing methods employ single-chain PMCMC. Because particles require centralized processing, the centralized resampling step in these methods is computationally intensive and time-consuming, significantly impacting computational speed. The present invention, however, employs a multi-chain PMCMC approach to predict model parameters and epidemic trends. Specifically, by introducing a local resampling method into the algorithm, N particles are evenly distributed across D running chains, allowing each chain to integrate sampling, weight calculation, resampling, and likelihood estimation, effectively addressing the issue of excessive particle computational density. This improvement significantly improves computational speed during epidemic prediction compared to previous methods, making it more adaptable to real-time application scenarios. Furthermore, to enhance the computational accuracy of multi-chain PMCMC, the present invention introduces a particle exchange mechanism after the resampling step. This mechanism alters the particle transmission channel, allowing particles of different weights to be exchanged during the local resampling process. This not only introduces more randomness into the particles in the system but also avoids the situation where prediction results do not converge due to excessive or insufficient particle weights in a particular chain, further improving the precision and accuracy of the epidemic prediction model. BRIEF DESCRIPTION OF THE DRAWINGS

[0057] Figure 1 The flowchart of the method according to the embodiment of the present invention is shown in FIG.

[0058] Figure 2 This is a diagram of the operation mode of the multi-chain PMCMC described in an embodiment of the present invention.

[0059] Figure 3 This is a diagram of the operation mode of adding a particle exchange mechanism as described in an embodiment of the present invention.

[0060] Figure 4 A complete flowchart modeling the method described in an embodiment of the present invention. DETAILED DESCRIPTION

[0061] The following is a detailed description of an embodiment of the present invention. This embodiment is based on the technical solution of the present invention, provides a detailed implementation method and a specific operation process, and further explains the technical solution of the present invention.

[0062] Example 1

[0063] The embodiment of the present invention provides a method for constructing an epidemic prediction model based on pedigree information and multi-chain PMCMC, and predicts the trend of the outbreak and the number of people according to the outbreak of the new coronavirus from time period A to time period B, thereby verifying the practicality of the method of the present invention. Figure 1 As shown, the specific steps include:

[0064] Step 1: Select the SEIR model as the state transition model for epidemic prediction and determine the model parameters to be estimated.

[0065] Since the nature of epidemiological dynamics is random and affected by multiple factors, this study introduces noise related to the change of transmission rate into the state transition model of the spatial state model (SSMs) using the SEIR model to simulate the dynamic changes of future epidemics. The present invention selects the SEIR model as the state transition model for epidemic prediction, as shown in formula (1):

[0066]

[0067] Where S t 、E t , I t and R t Represent the number of susceptible, exposed, infected, and recovered individuals at time t; A is the population size, and A=S t +E t +I t +R t ; ξ is noise, which is a normal random variable with mean zero and variance 1; F is the constant noise amplitude; β is the transmission coefficient, which measures the probability that an infected person will transmit the disease to a susceptible person in one contact, and β = R0 / γ, R0 is the basic reproduction number, γ is the recovery probability, γ = 1 / DI, DI represents the average infectious period; α is the rate at which individuals enter the infectious period from the latent period, α = 1 / DE, DE represents the average incubation period.

[0068] On the other hand, the observation model in spatial state models (SSMs) links the observed data to the underlying process model. Assuming that only a certain proportion ρ of new infections are detected and the observation error follows a normal distribution, then at a specific time t, y t The probability of a case can be expressed as formula (2):

[0069] P(y t |ΔI t )=N(y t |ρΔI t ,ρτΔI t ) (2)

[0070] Where ρ is the observation probability, τ is the scaling factor, and ΔI t For the number of new patients;

[0071] Assuming that the first case of novel coronavirus infection occurred d days before the actual discovery time, the model parameters to be estimated in the embodiment of the present invention can be determined based on the SSMs model, including: basic reproduction number R0, average incubation period DE, average infectious period DI, observation probability ρ, scaling factor τ, and the duration d of the epidemic virus from the first infection time to the actual discovery time.

[0072] Step 2: Initialize the model parameters to be estimated to obtain R0_initial, DE_initial, DI_initial, d_initial, ρ_initial and τ_initial, and use them as the model parameters θ0 of the current g=0 iteration.

[0073] Step 3: Calculate the likelihood function value based on the current model parameters.

[0074] Step 3.1, initialize the particle set, the total number of particles is N, the weight of each particle is the same, and the status of each particle includes the number of susceptible, exposed, infected and recovered people.

[0075] Step 3.2, divide the current particle set evenly into D parallel chains, such as Figure 2 As shown, the number of particles in each chain is P = N / D.

[0076] In step 3.3, based on the particle set and state transition model at time t-1, predict the state of each particle in each parallel chain at time t, and then calculate the predicted value of the number of newly infected people at time t.

[0077] For example, at time t=1, according to the system state X at time t=0 t-1 =[S t-1 ,E t-1 ,I t-1 ,R t-1] and the above state transition model to update the state X at the current time t=1 t =[S t ,E t ,I t ,R t ].

[0078] Then, based on the state X at the current time t=1 t =[S t ,E t ,I t ,R t ] and the observation equation shown in formula (3) to calculate the predicted value of new cases est_y at time t = 1 t :

[0079] est_y t =ρ*I t (3)

[0080] Where est_y t Represents the predicted value of the number of new epidemic cases at time t obtained based on the observation equation.

[0081] Step 3.4: Count the daily new cases of COVID-19 infection published on the official website, and calculate the particle weight W1 based on the actual new observed data and predicted data. The weight calculation part still divides the P particles into D parallel chains evenly and calculates them according to the time series method, as shown in formula (4):

[0082]

[0083] Where W1 represents the particle weight inferred from the time series, is the weight of particle i at time t, normpdf() is the normal probability density function, real t It represents the true value of the number of new epidemic infections at time t, and sqrt() represents the square root operation.

[0084] Step 3.5: Calculate the probability of a common ancestor event occurring at time t based on the true values of the number of susceptible and infected individuals in several subintervals at time t.

[0085] While calculating the weight W1 at the time series level in step 3.4, the embodiment of the present invention also reconstructs the past population dynamics of the disease from the pedigree information level. In epidemiology, especially in the study of transmission models and population genetics, it is necessary to use an exponential distribution to predict when new variants or key events in the spread of epidemics will appear, which is crucial for formulating public health strategies. This exponential distribution is called the common ancestor probability model, which is a random process used to describe how individuals are related to each other through common ancestors over a period of time. It traces the generations of individuals and establishes a "pedigree tree" between individuals. The tree shows the merging (common ancestor) process experienced from the current individual back to the common ancestor. As shown in formula (5):

[0086] P(T m |λ)=λe -λt (5)

[0087] This is an exponentially decaying function, Tm represents the time when the common ancestor event of the mth lineage branch occurs, λ is the common ancestor rate, which determines the speed at which the common ancestor event occurs, and t is the observation time (1 to T). In order to compare and combine with the inference results of PMCMC, the lineage must also be divided into T unit time intervals according to the observation time, and each unit time interval is divided into multiple sub-intervals (dt) so that the number of infected and susceptible individuals in the sub-interval changes slightly. This experiment divides each unit time interval (each day) into 24 sub-intervals, that is, dt = 1 / 24; Formula 5 describes the possibility of common ancestor. Among them, λ depends on the number of lineages present in the lineage m, the effective population size Ne, and the generation time τ of a disease. Among them, Ne depends on the number of infected hosts I dt According to the SEIR model, the disease generation time is more appropriately defined as the average length of time required for an infected individual to infect a susceptible host. The disease generation time is not constant over time because it depends on the rate at which infection occurs. Therefore, the common ancestor rate λ under the SEIR model can be expressed by formula (6):

[0088]

[0089] After obtaining the probability of common ancestor for each sub-interval, if k common ancestor events occur within a unit time interval, the probability of this time interval can be written as the joint probability of multiple common ancestor intervals. The probability of the family tree within a unit time interval W2 is shown in formula (7):

[0090]

[0091] Where W2 represents the probability of a common ancestor event occurring at time t, k and m represent the number and index of common ancestor events occurring at time t, respectively. mrepresents the time when the mth common ancestor event occurs within time t, dt represents the time intervals divided within time t, and λ dt represents the common ancestor rate in the time interval dt under the SEIR model, I dt S represents the number of infected people predicted by the SEIR model within the time interval dt, dt It represents the number of susceptible individuals predicted by the SEIR model within the time interval dt.

[0092] In step 3.6, multiply the weight of each particle at time t by the probability of a common ancestor event occurring at time t to obtain the weighted weight of each particle at time t.

[0093] The weight W1 of the latent variable inferred by the time series is combined with the possibility W2 of the family tree of the latent variable inferred by the pedigree information within a given time interval to obtain the comprehensive weight W = W1 * W2, which is hereinafter referred to as the weighted weight. At this time, the parameter adaptation rate is jointly determined by W1 and W2.

[0094] Step 3.7, based on the weighted weights of each particle at time t, calculate the likelihood function value at time t; the calculation formula of the likelihood function value is:

[0095]

[0096] Where likelihood represents the likelihood function, ∑W s (X 1:N ) represents the weighted summation of the weighted weights of all N particles in the particle set at time s.

[0097] In step 3.8, resample the particles in each parallel chain.

[0098] This embodiment adopts the local resampling method in the resampling step to evenly distribute N particles into D resampling chains, such as Figure 2 As shown in Figure 1, resampling is performed in parallel using multiple chains. The higher the weight of a particle, the higher its probability of accepting resampling. After resampling, new particles with equal weights are obtained, and the number is the same as before resampling. This multi-chain local resampling method can greatly improve the prediction speed.

[0099] In order to solve the problem of uneven weight distribution between particles during the iteration process, including the problem of decreased accuracy caused by excessive concentration on certain high-weight particles, a particle exchange mechanism can be introduced in a more optimal embodiment. Its function is to generate different permutations by adding a random permutation generator (RPG) to exchange particles from different data paths. After local resampling, the particle exchange mechanism will quickly generate a new random permutation and transfer particles from the corresponding data path according to the index value in the permutation, so that the particles in each running chain can be exchanged algorithmically. For example Figure 3There are four channels entering RPG. The particles in these four channels are numbered 00, 01, 10, and 11. Before entering RPG, the number of each particle is consistent with the original channel; after entering RPG, RPG randomly generates a right-shifted number and randomly shifts the original numbers 00, 01, 10, and 11. For example, if the right-shifted number generated is 3, the original four channel numbers are changed from 00011011 to the current channel number 01100011, that is, the particles in the original 00 channel are now transferred to 01; the particles in the original 01 channel are transferred to the current channel 10; the particles in the original 10 channel are transferred to the current channel 00, etc. This is a random process. Experiments have shown that this method makes the particle weights of each channel evenly distributed after the exchange.

[0100] The execution order of step 3.8 and step 3.7 can be interchanged or performed simultaneously.

[0101] Step 3.9: Repeat steps 3.3 to 3.8 until the likelihood function value at time T is calculated.

[0102] Step 3.10: Select one particle from the particle set and use its state value from time 1 to time T as the model parameter θ of the current g-th iteration. g The sample status value of .

[0103] After the likelihood estimate for a set of parameters is calculated, the state value of a particle at time 1:T is selected from all particle states at all times as the sample state value of the parameter. This sample state value needs to be updated and iterated along with the parameters, representing the number of new COVID-19 cases per day. Particles are selected based on their weight at time T before resampling, with particles with larger weights being more likely to be selected.

[0104] Step 4: Use the MCMC proposal distribution function to obtain the candidate parameters θ for the g+1th iteration of the model g+1 ', and calculate the candidate parameter θ based on step 3 g+1 'The likelihood function value.

[0105] Step 5: Based on the model parameters θ g The likelihood function value and candidate parameters θ g+1 'Likelihood function value, calculate the candidate parameter θ g+1 'Acceptance rate α; then determine whether to accept the candidate parameter θ based on the acceptance rate α g+1 'As the optimal parameter θ for the g+1th iteration of the model g+1 If not accepted, the current model parameter θ g As the optimal parameter θ for the g+1th iteration of the model g+1 .

[0106] Step 6: When the model iteration times are completed, the state transition model based on the optimal parameters at this time is the constructed epidemic prediction model.

[0107] The pseudo code of the above process is shown in Table 1:

[0108]

[0109]

[0110] Example 2

[0111] This embodiment provides an epidemic prediction method. First, an epidemic prediction model is constructed using the method described in Example 1. Then, the true values of the numbers of susceptible persons, exposed persons, infected persons, and recovered persons at a certain time t0 are obtained, wherein the certain time t0 is the time closest to the current time and the numbers of susceptible persons, exposed persons, infected persons, and recovered persons are all known. Finally, based on the true values of the numbers of susceptible persons, exposed persons, infected persons, and recovered persons at time t0, the constructed epidemic prediction model is used to perform a rolling prediction of the numbers of susceptible persons, exposed persons, infected persons, and recovered persons at the time to be predicted.

[0112] The following is a complete example process for trend prediction of new coronavirus infection, refer to Figure 4 , including the following steps:

[0113] Step 1-1, data preparation:

[0114] Before the prediction, we counted the number of new COVID-19 infections per day from time period A to time period B based on official website data. This data was used as observation data for constructing pedigree information and calculating time series particle weights.

[0115] Step 1-2, model determination:

[0116] For epidemic prediction, it is necessary to determine the model and method used for the epidemic trend prediction. According to the steps of the method described in Example 1, the state transition equation and observation equation of the spatial state model (SSM) are determined using pedigree information, multi-chain PMCMC method and SEIR model.

[0117] Steps 1-3, prior settings:

[0118] During the implementation of this example, you need to set the initial values of the parameters R0, DE, DI, d, ρ, and τ. All parameter results in subsequent iterations must be predicted according to the algorithm. In addition, in this example implementation, you need to set the particle state at the initial moment, the total number of particles N to 1024, the number of observations to 1109, the multi-chain parallelism D to 4, and the number of MCMC iterations M to 10,000.

[0119] Steps 1-4, prediction process:

[0120] This step uses pedigree information and a multi-chain PMCMC algorithm, using the SEIR model as the state transition and daily COVID-19 infection data as the observation to predict epidemic trends. The particle weights W1 of the time series and the probability W2 of the pedigree within a given time interval are calculated simultaneously, and the likelihood estimates of the parameters are calculated based on W1 and W2. Once the prediction is complete, not only is the probability density distribution of each parameter provided, but also the predicted value of the daily epidemic increase data can be obtained based on the parameter sample values.

[0121] Steps 1-5, Error Assessment:

[0122] The root mean square error is calculated based on the predicted data and the actual data obtained in steps 1-4. The root mean square error is used to measure the difference between the model predicted value and the actual observed value. If the error is large, it means that the predicted data is inaccurate and needs to be adjusted and optimized.

[0123] Steps 1-6, adjustment and optimization:

[0124] The data is adjusted by adjusting the total number of particles, the number of MCMC iterations, the parameter proposal distribution step size, the number of parallel chains, the initial value of the parameters, etc. After each adjustment, steps 1-3 to 1-5 need to be rerun until the accuracy is high enough before running the next step.

[0125] Steps 1-7, model evaluation:

[0126] When the prediction accuracy obtained in the above steps is high, it is necessary to compare the prediction time of this method with that of the existing PMCMC algorithm to evaluate the acceleration effect of the present invention; in addition, the method used in this invention should be implemented on FPGA and GPU, and compared with CPU to measure the advantages of this method in resources and speed.

[0127] The above embodiments are preferred embodiments of the present application. Ordinary technicians in this field can also make various changes or improvements on this basis. Without departing from the overall concept of the present application, these changes or improvements should fall within the scope of protection required by the present application.

Claims

1. A method for constructing an epidemic prediction model based on pedigree information and multi-chain PMCMC, characterized in that: include: Step 1: Select the SEIR model as the state transition model for epidemic prediction and determine the model parameters to be estimated; Step 2: Initialize the model parameters to be estimated and use them as the model parameters θ0 of the g=0th iteration; Step 3, calculate the likelihood function value based on the current model parameters; Step 3.1, initialize the particle set, the state of each particle includes the number of susceptible, exposed, infected and recovered people; Step 3.2, divide the current particle set evenly into D parallel chains; Step 3.3: Based on the particle set and state transition model at time t-1, predict the state of each particle in each parallel chain at time t, and then calculate the predicted value of the number of newly infected people at time t; Step 3.4: Calculate the weight of each particle in each parallel chain at time t based on the predicted and actual number of newly infected people at time t. Step 3.5, calculate the probability of a common ancestor event occurring at time t based on the true number of susceptible and infected individuals in several subintervals at time t; Step 3.5 is performed after or simultaneously with Steps 3.3-3.4; Step 3.6, multiply the weight of each particle at time t by the probability of the common ancestor event occurring at time t to obtain the weighted weight of each particle at time t; Step 3.7, based on the weighted weights of each particle at time t, calculate the likelihood function value at time t; Step 3.8, resample the particles in each parallel chain; wherein, the execution order of step 3.8 and step 3.7 can be interchanged or performed simultaneously; Step 3.9: Repeat steps 3.3 to 3.8 to calculate the likelihood function value at time t+1 until the likelihood function value at time T is calculated; Step 3.10: Select one particle from the particle set and use its state value from time 1 to time T as the model parameter θ of the current g-th iteration. g The sample status value of Step 4: Use the MCMC proposal distribution function to obtain the candidate parameters θ for the g+1th iteration of the model g+1 ', and calculate the candidate parameter θ based on step 3 g+1 'likelihood function value; Step 5: Based on the model parameters θ g The likelihood function value and candidate parameters θ g+1 'Likelihood function value, calculate the candidate parameter θ g+1 'Acceptance rate α; then determine whether to accept the candidate parameter θ based on the acceptance rate α g+1 'As the optimal parameter θ for the g+1th iteration of the model g+1 If not accepted, the model parameter θ of the gth iteration g As the optimal parameter θ for the g+1th iteration of the model g+1 ; Step 6: Repeat steps 3 to 5 until the model iterations are completed. The state transition model based on the optimal parameters at this time is the constructed epidemic prediction model.

2. The method for constructing an epidemic prediction model based on pedigree information and multi-chain PMCMC according to claim 1, characterized in that: The SEIR model is selected as the state transition model for epidemic prediction, which is expressed as: I t+1 =I t +αE(t)-γI t R t+1 =R t +γI t Where S t 、E t , I t and R t Represent the number of susceptible, exposed, infected, and recovered individuals at time t; A is the population size, and A=S t +E t +I t +R t ; ξ is noise, which is a normal random variable with mean zero and variance 1; F is the constant noise amplitude; β is the transmission coefficient; α is the rate at which individuals enter the infectious period from the latent period; γ is the recovery probability; among them, β = R0 / γ, R0 is the basic reproduction number; γ = 1 / DI, DI represents the average infectious period; α = 1 / DE, DE represents the average latent period.

3. The method for constructing an epidemic prediction model based on pedigree information and multi-chain PMCMC according to claim 2, characterized in that: The observation error of epidemic prediction follows a normal distribution, so when y is observed at time t t The probability of a case is expressed as: P(y t |ΔI t )=N(y t |ρΔI t ,ρτΔI t ) Where ρ is the observation probability, τ is the scaling factor, and ΔI t For the number of new patients; The model parameters to be estimated include: basic reproduction number R0, average incubation period DE, average infectious period DI, observation probability ρ, scaling factor τ, and the length of time d from the first infection time to the actual discovery time of the epidemic virus.

4. The method for constructing an epidemic prediction model based on pedigree information and multi-chain PMCMC according to claim 3, characterized in that: The calculation formula for the predicted number of new infections at time t is: is_y t =ρ*I t Where est_y t represents the predicted value of the number of new epidemic cases at time t obtained according to the above observation equation.

5. The method for constructing an epidemic prediction model based on pedigree information and multi-chain PMCMC according to claim 1, characterized in that: Step 3.4 calculates the weight of each particle in each parallel chain at time t, expressed as: Where, is the weight of particle i at time t, normpdf() is the normal probability density function, real t represents the true value of the number of new epidemic infections at time t, est_y t It represents the predicted number of new epidemic infections at time t, and sqrt() represents square root operation.

6. The method for constructing an epidemic prediction model based on pedigree information and multi-chain PMCMC according to claim 2, characterized in that: Step 3.5 calculates the probability of a common ancestor event occurring at time t, expressed as: Where W2 represents the probability of a common ancestor event occurring at time t, k and m represent the number and index of common ancestor events occurring at time t, respectively. m represents the time when the mth common ancestor event occurs within time t, dt represents the time intervals divided within time t, and λ dt represents the common ancestor rate in the time interval dt under the SEIR model, I dt S represents the number of infected people predicted by the SEIR model within the time interval dt, dt It represents the number of susceptible individuals predicted by the SEIR model within the time interval dt.

7. The method for constructing an epidemic prediction model based on pedigree information and multi-chain PMCMC according to claim 1, characterized in that: The calculation formula for the likelihood function value in step 3.7 is: Where likelihood represents the likelihood function, ΣW s (X 1:N ) represents the sum of the weighted weights of all N particles in the particle set at time s.

8. The method for constructing an epidemic prediction model based on pedigree information and multi-chain PMCMC according to claim 1, characterized in that: Step 3.6 performs resampling in parallel through multiple chains. The higher the weight of a particle, the higher its probability of accepting resampling. After resampling, new particles with equal weights are obtained and the number is the same as before resampling. After resampling, a particle exchange mechanism is introduced between different chains, that is, particles in different chains are randomly exchanged.

9. A method for predicting an epidemic, characterized in that: include: Constructing an epidemic prediction model using the method described in any one of claims 1 to 8; Get the true value of the number of susceptible, exposed, infected and recovered people at a certain time t0; The time t0 is the time closest to the current time and the number of susceptible persons, exposed persons, infected persons, and recovered persons is known; According to the actual values of the number of susceptible people, exposed people, infected people and recovered people at time t0, the constructed epidemic prediction model is used to make a rolling prediction of the number of susceptible people, exposed people, infected people and recovered people at the time to be predicted.

Images