Skip to content

Instantly share code, notes, and snippets.

@t-nissie
Last active March 18, 2016 02:36
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save t-nissie/fe0cbab32e8472103652 to your computer and use it in GitHub Desktop.
Save t-nissie/fe0cbab32e8472103652 to your computer and use it in GitHub Desktop.
エボラ出血熱の時刻t=0に感染した後、時刻t>0で発病する確率が潜伏期間の確率密度関数f(t)。μ=10、σ=6として水素原子の動径波動関数に似た関数で描いてみた。21日経っても5%がまだ発症しない。
#!/usr/bin/env gnuplot
# ebola.gp: Probability density function f(t) of incubation period
# Time-stamp: <2016-03-18 11:34:36 takeshi>
# Author: Takeshi Nishimatsu
# Ref1: WHO Ebola Response Team: N. Engl. J. Med. 37, 1481-1495 (2014), doi:10.1056/NEJMoa1411100 .
# Ref2: Hiroshi Nishiura: Emerging Themes in Epidemiology 2007, 4:2, doi:10.1186/1742-7622-4-2 .
# Gamma distribution: en: http://en.wikipedia.org/wiki/Gamma_distribution
# ja: http://ja.wikipedia.org/wiki/%E3%82%AC%E3%83%B3%E3%83%9E%E5%88%86%E5%B8%83
# Log-normal distribution: en: http://en.wikipedia.org/wiki/Log-normal_distribution
# ja: http://ja.wikipedia.org/wiki/%E5%AF%BE%E6%95%B0%E6%AD%A3%E8%A6%8F%E5%88%86%E5%B8%83
# Gist: https://gist.github.com/t-nissie/fe0cbab32e8472103652
# Clone: git clone git@gist.github.com:fe0cbab32e8472103652.git ebola
##
mu = 10.0
sigma2 = 36.0
normal(x) = exp(-(x-mu)**2/2/sigma2)/sqrt(sigma2)/sqrt(2*pi)
ti=sprintf("{/Symbol m}=%5.2f, {/Symbol s}^2=%5.2f",mu,sigma2)
set title ti
set sample 600
a = mu**2/sigma2-1
b = mu/sigma2
c = gamma(a+1)/b**(a+1)
print 'a= ',a, ' b= ',b, ' c= ', c
f(x) = x**a*exp(-b*x)/c
F(x) = igamma(a+1,b*x)
print 100.0*(1.0-F(21.0))
l_sigma2 = log(sigma2/mu**2+1)
l_mu = log(mu)-l_sigma2/2
print 'l_mu= ', l_mu, ' l_sigma2= ', l_sigma2
log_normal(x) = exp(-(log(x)-l_mu)**2/2/l_sigma2)/sqrt(l_sigma2)/sqrt(2*pi)/x
set terminal postscript eps enhanced color "Times-Roman" 28
set output "ebola.eps"
set grid
set size 1.0,2.0
set multiplot
#set key at 4.9,0.9 spacing 1.3
set xrange [-5:25]
set format y '%.2f'
set origin 0.0,1.0
set size 1.0,1.0
set yrange [0:0.1]
set ytics 0.02
plot f(x) title '{/Times-Italic f}({/Times-Italic t}) = {/Times-Italic t^a e^{-bt} / c}' \
with lines lt 1 lw 3,\
normal(x) title 'normal dist.' with lines lt 3 lw 1,\
log_normal(x) title 'log-normal dist.' with lines lt 4 lw 1
set origin 0.0,0.0
set size 1.0,1.0
set tmargin 0
set yrange [0:1.0]
set ytics 0.1
set key at 12,0.88
set notitle
set xlabel "{/Times-Italic t} [days after infection]"
plot F(x) title '{/Times-Italic F}({/Times-Italic t}) = {/Symbol \362}@^{/Times-Italic t}_0 {/Times-Italic f}({/Times-Italic t}{/Symbol \242}){/Times-Italic dt{/Symbol \242}}' \
with lines lt 2 lw 3
#Local variables:
# compile-command: "gnuplot ebola.gp"
#End:
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment