sdwfrost

sdwfrost

sdwfrost

#LMS change Detector 
function [state,h]=ChangeDetectionLMSShiftMultipleSensor1(x,WinSlow,WinFast,WinDesired,muBase)
%Slow Response to change, accurate in steady state

[m,n]=size(x);
if m>n
   x=x';
end
[NumSensors,n]=size(x);
shift=300;

sdwfrost

# Ordinary differential equation model using Petri.jl
Simon Frost (@sdwfrost), 2020-07-08

## Introduction


This implementation considers the SIR model as a Petri net, based on the example of the [documentation of `Petri.jl`](https://mehalter.github.io/Petri.jl/stable/usage/), which is then used to generate an ordinary differential equation model.

## Libraries

sdwfrost

# Libraries
using ModelingToolkit
using OrdinaryDiffEq

# Define size of the system
nstates = 3 # S,I,R
nrates = 2 # infection, recovery
# Model specific parameters
@parameters β c γ S₀ I₀ R₀

sdwfrost

# Libraries
using ModelingToolkit
using OrdinaryDiffEq
using LinearAlgebra: Diagonal
import Base: sqrt

# Utility function

function sqrt(x::Diagonal)
    Base.sqrt.(x)

sdwfrost

---
title: "SIR in phydynR"
author: "Simon Frost"
date: '`r Sys.Date()`'
output: github_document
---

```{r setup, include=FALSE}
library(simecol)
library(phydynR)

sdwfrost

Branching process model of 2019nCoV spread

This is an attempt to reproduce Figure 1 from Imai et al. 2020, using code from Riou and Althaus (code here). I have included a Jupyter notebook (which includes the output), as well as RMarkdown and plain R files, generated using Jupytext.

Details missing from Imai et al. that I had to guess:

1. The initial time (I took it to be 2019-12-01).
2. The distribution of generation times (I took it to be gamma with standard deviation of 3.8).
3. The summary curve is generated using the median of all cases at a given time (similarly with the 95 percentile range). The language in the figure legend is ambiguous, and could mean generating a median trajectory.

sdwfrost

#!/usr/bin/env run-cargo-script
//!
//! ```cargo
//! [dependencies]
//! time = "0.1.25"
//! float_extras = "0.1.1"
//! xorshift = "0.1"
//! ```

#![allow(non_snake_case)]

sdwfrost

// zig build-exe sir.zig --release-fast --strip
const std = @import("std");

inline fn randbn(n: i64, p: f64, rng: *std.rand.Xoroshiro128) i64 {
    const log_q: f64 = std.math.ln(1.000000 - p);
    var x: i64 = i64(0);
    var sum: f64 = 0.000000;
    while (1 != 0) {
        sum += std.math.ln(rng.random.float(f64))/@intToFloat(f64, n - x);
        if (sum < log_q) break;