MATLAB课业代写 Final Project代写 model代写 ODE代写

ENG3 Final Project

MATLAB课业代写 Modeling the Covid-19 pandemic The SIR (susceptible-infected-removed) epidemiological model consists of a simplified mathematical model

Modeling the Covid-19 pandemic

The SIR (susceptible-infected-removed) epidemiological model consists of a simplified mathematical model for the spread of infectious diseases. The idea of the SIR model is simple: consider a population of N individuals. We then separate the population into three groups: susceptible, infected and removed (deceased). The sizes of these subpopulations at time t are denoted by S(t), I(t), and R(t), respectively. Mathematically, this model can be written as a system of three coupled, non-linear ordinary differential equations (ODEs)

where we assume that the disease transmission rate β = β(t) and the mortality rate γ = γ(t) are real, time-dependent functions. Given that the ODEs are nonlinear and coupled, finding an analytical solution for Eqs. (1)-(3) is impractical. Nonetheless, one can search for solutions numerically, using, for example, ODE solvers such as the ones available in Matlab. We will use this classic epidemiological model to model the spread of Covid-19 in Santa Barbara County.

Part 1: Exploring the data (10 points) MATLAB课业代写

The first step of tackling a problem such as this is getting to know the data! Therefore, your first task consists of loading the data file covid_data.mat provided on GauchoSpace and plotting its contents.

The file consists of a dataset for COVID-19 cases and deaths for the 90 days following March 15 in a fictional Santa Barbara county. It has 2 columns: the first contains the cumulative number of infected individuals in the population at a certain date, and the second contains the cumulative number of deaths recorded on the same date.

(a)

Create a script named explore_data.m that loads the file covid_data.mat provided.

(b)

Still in your explore_data.m script, create a figure with two plots: i) number of infected individuals vs time and ii) deaths vs time. Make sure that you:

• Change the line thickness to 3 for each line.

• Use legends “Infected Individuals”, “Deaths” to identify each curve.

• Name the x and y axes as “Time (days)” and “Number of individuals”, respectively. You can use legend(’location’,’northwest’) to move your legend box to the top left corner, so it doesn’t cover your plot. MATLAB课业代写

• Give your plot the following title: “Covid-19 data”.

• Turn the grid on using the command grid on.

• Increase the font size of legends, axis labels and title using the command set(gca,’FontSize’,15).

(c)

As you can see from the figure plotted in (b), the length scales used to plot make it difficult to discern changes in the number of deaths over time. Let’s try a different approach. Create another figure with the same two plots as in (b), but this time using a log scale on the y axis.

• Follow all the same instructions for this figure as those in (b), but use the semilogy command to plot instead of the plot command.

At the end, your plots should look like the ones displayed in Fig. 1. Upload explore_data.m to GauchoSpace.

Part 2: Estimating the parameters from the data (20 points) MATLAB课业代写

Now that you have familiarized yourself with the data, let us come back to the SIR model. Our goal is to use it to predict the growth of the number of infected individuals as well as the number of deaths in the near future. But the problem is that we don’t know the values for the functions β(t), γ(t) in Eqs. (1)-(3). Despite that, observing the data from places where the Covid-19 epidemic got under control via quarantine measures can give us some hints about these parameters. For example, it is reasonable to consider that both the infection and the mortality rates decreased as more and more people started taking precautionary measures, e.g., stay-at-home shelters, self-quarantine, crowd control in public transportation systems and the usage of masks to avoid the spread of the virus.

Based on this information, let us postulate that both the infection and the mortality rate decreased exponentially after the quarantine measures were implemented, i.e.,

γ(t) = γ₀ e^−0.0199t,   (4)

β(t) = (7.50 × 10⁻⁷) e^−0.05088t,   (5)

where γ₀ is a constant to be determined. Your task is to estimate γ₀ using the data you plotted in Part 1.

There exist many methods to estimate parameters (like γ₀ from a given set of data points. For simplicity, you are going to use a Euclidean norm approach to minimize the error in the total number of deaths R(t). The idea is simple: solve Eqs. (1)-(3) for various γ₀ values within a range, and for each of these values, compute the error given by the norm of the difference between the model results R(t) and the data, i.e.,

where N days is the number of days covered by the data (i.e., the length of the covid_data matrix).

(a)

Write a function named sir_model.m that takes as input t (time) and X (a vector of S, I, R), and returns the right-hand side of the ODE system given by Eqs. (1)-(3).

(b)

Write a script named sir_main.m that solves the ODE system for various values of γ₀ (use a for loop for that).

• Use the ode45() function to solve the ODE system. For the initial conditions, consider that the susceptible population of Santa Barbara County was 500,000, and assume that the infection started with 1 person infected and zero deaths. Use a time span of 90 days.

• Use the following range of values for γ₀ : gamma0 = (1:10:1000)./S0, where S0 is the initial susceptible population

• Compute the error using Eq. (6) for each value of γ₀ tested, and store these values.

(c)

Still in your script sir_main.m, plot γ0 vs err (using the plot() function). Make sure that:

• Your plots look like Fig. 2.

• Give your plot the following title: “Error analysis (SIR model vs data)”.

• Change the line thickness to 3.

• Name the x and y axis as “gamma0” and “Error”, respectively.

• Turn the grid on using the command grid on.

• Increase the font size of legends, axis labels and title using the command set(gca,’FontSize’,15).

(d)

Still in your script sir_main.m, plot the number of deaths R(t) vs. time for the best value of γ₀ found in part (c). In addition, overlap the curves for R with the number of deaths in the covid_data.mat data provided on GauchoSpace, so that you can compare simulations with data. Make sure that:

• Your plots look like Fig. 3.

• Change the line thickness to 3 for each line.

• Give your plot the following title: “SIR model, Estimated Gamma0 = X”, where X must be replaced by the respective value of β₀ found in (c).

• Use legends “Deaths (SIR)”, “Deaths (data)” in your plots.

• Name the x and y axis as “Time (days since March 15)” and “Number of individuals”, respectively. You can use legend(’location’,’best’) to move your legend box to the best available spot, so it doesn’t cover your plots.

• Turn the grid on using the command grid on.

• Increase the font size of legends, axis labels and title using the command set(gca,’FontSize’,15).

Note: The best γ₀ value obtained in (c) should be close to γ₀ ≈ 0.0012. If you have not been able to accomplish the optimization, please proceed with γ₀ = 0.0012 for Part 3 of this project.

Upload sir_model.m and sir_main.m to GauchoSpace.

Part 3: Stochastic simulation of SIR model (20 points) MATLAB课业代写

Let’s take a look back at your plots from Part 2-(c). At this point you probably noticed that we can find a γ0 that fits the data reasonably well, but no matter how much you change it, we will never find the perfect fit. One of the reasons for that involves the fact that traditional deterministic (ODE) models like the one given by Eqs. (1)-(3) cannot accurately predict the early stages of the infection, since they rely on the hypothesis of a well-mixed population and interactions of millions of individuals. In contrast, for small populations, random events are often relevant, and therefore situations like jumps in the number of cases are prone to occur; deterministic models, on the other hand, disregard these events. To solve these problems, we will use stochastic models.

Notice that the current number of infected individuals in Santa Barbara County is relatively small compared to the total population. Thus, the deterministic model given by Eqs. (1)-(3) might lead to incorrect predictions. Your goal for this assignment is to solve a stochastic SIR model using a method called the Gillespie algorithm.

The Gillespie algorithm belongs to a class of methods called dynamic Monte Carlo methods. As you probably remember from homeworks 6 and 7, Monte Carlo simulations allow us to estimate probabilities of events by testing multiple trials (also called realizations). Since each realization involves random numbers, the results of Monte Carlo simulations are never exactly the same. Similarly, each realization (simulation) via the Gillespie algorithm will generate unique curves for S(t), I(t), R(t).

The stochastic SIR model considers the following chains of events (also called Markov chains)

where γ(t), β(t) are the same functions given by Eqs. (4)-(5). The chain given by Eq. (7) describes the process of infection: each time a susceptible individual meets an infected one, there’s a propensity (which is proportional to β(t)) that dictates how often the susceptible individual becomes infected. Similarly, the chain in Eq. (8) describes the process of infected individuals succumbing to the virus and dying, with a propensity proportional to γ(t). The propensity function α_i  gives the probability that an event in the chain i will occur in the time interval (t, t + τ ):

α_i = number of individuals interacting at chain ”i” × rate,  (9)

which, for the chains given by Eq. (7)-(8), read

α₁ = S(t)I(t)β(t),    (10)

α₂ = I(t)γ(t). (11)

For more details about the implementation of this algorithm, see p. 13 of this final project assignment: Recipe for Stochastic SIR using Gillespie Algorithm.

(a)

Using the recipe provided on p. 13, write a function named ssir_model.m, that takes as input tspan (time span) and X0 = [* ; * ; *] (column vector of initial conditions for S, I, R) and returns one realization (i.e., column vectors containing S(t), I(t), R(t), resulting from the algorithm provided in p. 9) of the Gillespie algorithm, and a time vector tVec. Use the functions γ(t), β(t) given by Eqs. (4)-(5), and use the value of γ₀ that you found best fit the data in Part 2-(c), i.e., γ₀ = 0.0012 .

(b)

Write a script named ssir_main.m to compute Nr = 5 realizations of the stochastic SIR model by calling your function ssir_model 5 times (use a for loop for that). For the initial conditions, assume that the susceptible population of Santa Barbara County was 500,000, and assume that the infection started with 1 person infected and zero deaths. Use a time span of 90 days.

(c) MATLAB课业代写

Still in your script ssir_main.m, plot the number of deceased individuals R(t) vs. time, for each of the 5 realizations computed in Part 3-(b). In addition to that, overlap the curves for R with the covid_data.mat data provided on GauchoSpace, so that you can compare simulations with data. Make sure that:

• Your plots, legends and title look like Fig. 4. Note: since the method is

stochastic, your curves will be different than the ones shown in Fig. 4.

• Instead of using the usual plot(), use the command semilogy() to plot in semi-log scale. This will make your two curves look more distinguish able.

• Use a linewidth of 3 for the sSIR model curves and a markersize of 15 for the covid death data.

• Give your plot the following title: “sSIR model”.

• Use legends ’Deaths (sSIR)’, ’Deaths (data)’ in your plots. Note: Since all the realizations have the same legends, use ’HandleVisibility’,’off’ in your semilogy() commands for realizations 2 to Nr, to avoid legend repetition.

• Name the x and y axis as “Time (days since March 15)” and “Number of individuals”, respectively. Again, you can use legend(’location’,’best’) to move your legend box to the best available spot, so it doesn’t cover your plots.

• Turn the grid on using the command grid on.

(d)

In a similar manner as you did for your script ssir_main.m, create another script, named ssir_hist.m, to compute Nr = 2000 realizations of the stochastic SIR model. At the end of each realization i, store the last value of the total number of deaths, i.e., deaths(i) = R(end). This value represents a prediction of the number of deaths in Santa Barbara County after tspan days (i.e., 90 days). With these values, create a histogram of the array deaths using the command histogram(deaths). Your histogram should look similar to Fig. 5.

For the plot: Use the options ’FaceColor’ and ’EdgeColor’ in the histogram command to create the histogram plotted in blue with black edges.

Upload ssir_model.m, ssir_main.m, ssir_hist.m to GauchoSpace.

Recipe for Stochastic SIR using the Gillespie Algorithm MATLAB课业代写

While t < tspan , do the following:

1. Generate two random numbers, r₁, r₂, uniformly distributed in the interval (0, 1).

2. Compute the propensity function α for each chain. Since we have two chains, the propensities are given by

α₁ = S(t)I(t)β(t),   (12)

α₂ = I(t)γ(t), (13)

where γ(t), β(t) are the functions given by Eqs. (4)-(5). Use the value of β₀ that you found best fit the data in Part 2-(c).

3. Compute the total propensity α₀ = α₁ + α₂.

6. Make sure that the system always includes at least one infected person, i.e., check if I(t + τ ) > 0. If I(t + τ ) becomes 0, include another infected person at time t + τ by setting I(t + τ ) = 1.

7. Save your new time as a point in the time vector: tVec(i) = t + tau (so you can plot S(t), I(t), R(t) vs t later).

更多代写：计算机新西兰网课代上代修  GMAT代考  英国统计学代考  英语essay代写北美  新西兰音乐学论文代写 英国网课代写

合作平台：essay代写 论文代写 写手招聘 英国留学生代写