Last active
March 18, 2016 02:36
-
-
Save t-nissie/fe0cbab32e8472103652 to your computer and use it in GitHub Desktop.
エボラ出血熱の時刻t=0に感染した後、時刻t>0で発病する確率が潜伏期間の確率密度関数f(t)。μ=10、σ=6として水素原子の動径波動関数に似た関数で描いてみた。21日経っても5%がまだ発症しない。
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/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