Hatena::ブログ(Diary)

驚異のアニヲタ社会復帰への道

Prima Project

2016-12-07

MikuHatsune2016-12-07

球面調和関数の回転

球面調和関数RでやったPythonでやったりしていたのだが、球面調和関数展開して係数をもっているときに、物体が回転してしまうとルジャンドル関数が回転してしまうので、係数も回転に応じて適当に変換しないといけない。

これはたいてい、回転後のパラメータがカンマつきで表すと

f(¥theta^{’},¥phi{’})=¥sum_{l=0}^¥infty ¥sum_{m{’}=-l}^l c_l^{m^{’}} Y_l^{{’}}(¥theta{’},¥phi{’})

という、回転後の係数c_l^{m^{’}} で考えないといけない。ちなみに回転は、オイラー(¥alpha, ¥beta, ¥gamma) で、zyz 軸がどの教科書や論文をみてもそうかいてある。

回転自体はこれがわかりやすい。

 

回転しても、もとの係数c_l^m から計算できるように偉い人が考えていて、

c_l^{m^{’}}=¥sum_{m=-l}^l D_{l}^{m^{’},m} c_l^m

D_l^{m^{’},m}=e^{im¥gamma}d_l^{m^{’},m}(¥beta)e^{im¥alpha}

d_l^{m^{’},m}(¥beta)=¥sum_{t=max(0,m-m^{’})}^{min(l+m,l-m^{’})} (-1)^t ¥frac{¥sqrt{(l+m)!(l-m)!(l+m^{’})!(l-m^{’})!}}{(l+m-t)!(l-m^{’}-t)!(t+m^{’}-m)!t!}(¥cos¥frac{¥beta}{2})^{2l+m-m^{’}-2t}(¥sin¥frac{¥beta}{2})^{2t+m^{’}-m}

とすごいことになっている。

R でがんばってみたが、計算は遅いし結果もおかしいしで、Python に頼る。pyshtools はもともとfortran でshtools というものがあるっぽいのだが、それのPythonモジュールである。を参考にやってみる。

 

import pyshtools as p
import numpy as np
import math
import itertools

pi = math.pi

# 球面調和関数係数をおさめる箱を適当につくる
x = p.SHCoeffs.from_random([10, 1, 0])
x.get_coeffs()

lmax = 2
x = p.SHCoeffs.from_zeros(lmax)
ls = [0, 1, 1, 1, 2, 2, 2, 2, 2]     # l の番地
ms = [0, -1, 0, 1, -2, -1, 0, 1, 2]  # m の番地
value = [10, 1, 2, 0, 1, 5, 3, 4, 0] # (l, m) の値

# sphericalharmonics coefficients
for (l, m, v) in zip(ls, ms, value):
  x.set_coeffs(v, l, m)

x.get_coeffs() でnumpy array float 型のリストを取得できる。l=2 のデータを作ってみると、長さ2、中身はl+1 * l+1 次元の行列になっている。

array([[[ 10.,   0.,   0.],
        [  2.,   0.,   0.],
        [  3.,   4.,   0.]],

       [[  0.,   0.,   0.],
        [  0.,   1.,   0.],
        [  0.,   5.,   1.]]])

リストの最初はm¥geq 0 の項で、2つめはm<0 の項である。ひとつめのリストには、l=0,1,2 の要素がはいっていて、各l のリストにはm=0,1,2 の要素が入っている。

m<0 のほうは少し変わっていて、各l のなかにm=0,-1,-2 の順番で要素がはいっている。要素をいじるにはset_coeffs() で値、l, m の番地を指定すれば任意にいじることができる。

係数を取り出したかったら

def shc_extract(coef, lmax):
  c = [coef.get_coeffs()[0][0][0]]
  for i in range(1, lmax+1):
    for k in [1, 0]:
      if k == 0:
        for j in range(k, i+1):
          c += [coef.get_coeffs()[k][i][j]]
      else:
        for j in reversed(range(k, i+1)):
          c += [coef.get_coeffs()[k][i][j]]
  return c

shc_extract(x_rot, lmax)

f:id:MikuHatsune:20161207144806p:image

 

球面調和関数は球面上のポテンシャルの(数学的には厳密な)級数展開なので、球面上のなにかを比較したかったら、係数を比較してもいいことになっている。もとの形を球面調和関数で記述してあるとして、適当にぐるぐる回転してみて、もっとも近いもしくは遠い回転を探してみる。

たぶんこれも、何かいい方法があるのだろうが、愚直に回転させまくってみる。10度ずつの回転をzyz に対して総当りする。総当りにはitertools が使い勝手がよく、めっちゃ速い。ちなみに上でzip を使って一括iterate しているが、これもパクっている。

 

div = 10                  # 何度ずつ回転するか
degs = range(1, 360, div) # 回転角度
rot_res = []
for i, j, k in itertools.product(degs, degs, degs):
  irad, jrad, krad = i*pi/180, j*pi/180, k*pi/180 # 回転角度をradian で記録しておく
  x_rot = x.rotate(irad, jrad, krad, degrees=False)
  rot_res += [[irad, jrad, krad, shc_extract(x_rot, lmax)]]


# R じゃないといろいろできない情弱なのでファイルに書き出す
w0 = open("rotate_shc.txt", "w")
for z in rot_res:
  txt = " ".join(map(str, z[:3])) + " " + " ".join(map(str, z[3])) + "\n"
  w0.write(txt)

w0.close()

 

実際は単位球面上に模様のように色が付いているが、その値に応じて半径を決めると、球面調和関数で画像をググったときにでてくるようなチョウチョみたいなアレができる。

差が最小の時の係数の回転である。ほとんど重なっているように見えるが、これでも(321, 1, 41)度回転しているらしい。

f:id:MikuHatsune:20161207144807p:image

 

差が最大の時の係数の回転である。このときは(51, 161, 131)度回転している。

f:id:MikuHatsune:20161207144808p:image

 

回転自体は、どれくらい細かくグリッドサーチするかと、d_l^{m,m^{’}}O(l^3) のオーダーで計算量がかかるらしいので、d_l^{m,m^{’}}再帰的に求める方法もあるらしい。

Journal of Computational Physics 228(16):5621-5627 September 2009

 

d_l^{m,m^{’}} を計算したかったらdjpi2 という関数が使える。これを使うと、最終的な回転はrotate 関数でできるが、

alpha = pi/8
beta = pi/9
gamma = pi/10
ang = [pi/1.2, pi/1.8, pi/3]
dj = p.shtools.djpi2(lmax)

x_rot = p.shtools.SHRotateRealCoef(x.get_coeffs(), np.array(ang), dj)
# 同じ
x_rot = x.rotate(alpha, beta, gamma, degrees=False)

とやっても同様な結果が得られる。

2016-12-06

MikuHatsune2016-12-06

フィッシャーの正確確率検定とカイ自乗検定と尤度比検定

遺伝統計学の基礎―Rによる遺伝因子解析・遺伝子機能解析―

遺伝統計学の基礎―Rによる遺伝因子解析・遺伝子機能解析―

3つの検定法の比較というところで、2*2 の分割表に対して正確確率検定、カイ自乗検定、尤度比検定の差を考えている。

分割表といえば、もっとも単純なものは

¥begin{array}{cc|c}a&b&e¥¥c&d&f¥¥¥hline g&h&n¥end{array}

というもので、周辺度数のe, f,, g, h, n を固定して、各項のa, b, c, d がどんな具合か、を調べる。

 

Fisher の正確確率検定は、前にもやったけど周辺度数が決まっているときに、考えられる分割表のパターンを列挙して、その生起確率を求める。生起確率は

P=¥frac{e!f!g!h!}{n!a!b!c!d!}

となる。階乗が含まれるので計算機的には対数をとる。周辺度数が大きくなると取りうる分割表のパターンが増えるので、計算量としては厳しい。しかし、すべての場合の数を網羅するので、正確と言われる。

Rによるフィッシャーの直接確率検定

 

カイ自乗検定は、分割表に対して検定をすれば、2群が独立なのかそれともなんらかの関係をもって各項が分布しているのかを検定するし、1次元ベクトルに使えば、例えばメンデルアサガオの色的なあれで言えばどれくらい想定している確率に従って分布しているかの適合度検定とかに使われる。これはどちらも、周辺度数(か想定確率)によって各項の期待値というものが計算できて、そこからどれくらいずれているかを定量化してp値にしている。ずれは自乗値で計算されていて、そこから分割表の行と列数に応じた自由度が決まって、カイ自乗分布から統計量→p値への変換がなされる。

各項の値をn_{ij}, 各項の期待値をe_{ij} とすると、統計量は

¥chi ^2=¥sum_{i,j}¥frac{(n_{ij}-e_{ij})^2}{e_{ij}}

となる。e_{ij} は周辺度数から決まり、予測値なので整数でなく実数をとるので¥chi ^2は連続値となる。

周辺度数が小さいときにはe_{ij} がかなりおおざっぱな範囲しか担当しないため、分子もそれ相応にいじらないと、値がずれる。しかも、¥chi^2 が大きくなる方向にずれるため、結果としてp値が過小評価される傾向にある。このため、周辺度数が小さいとか、各項に5以下の要素がある場合は、補正しましょうとかなんとか言われる。

正確確率検定のように、すべての場合の数は計算せず、周辺度数から決まる各ij 項分の計算をすればいいので、高速である。ただし、周辺度数が小さい場合は誤差が大きい。しかし、周辺度数がそこそこ大きくなれば理論値に近づくので、よい。でも、本当に理論値が欲しい時はアレ。

 

尤度比検定は、比較したいふたつの仮説の尤度比をとって、それがあまりにも偏っていたらその値が大きくなるので、そうなった場合は基準となる仮説を棄却しましょう、という考え。

例えば分割表では、帰無仮説H_0 は「どちらの行も同じ確率で生じたデータ」、対立仮説H_1 は「ふたつの行のデータを生じた確率は異なる」というのが一般的である。

H_0H_1 の尤度比L はそれぞれ

L_{H_0}=¥frac{(a+c)!}{a!c!}¥frac{(b+d)!}{b!d!}(¥frac{a}{a+b})^a (¥frac{b}{a+b})^b (¥frac{c}{c+d})^c (¥frac{d}{c+d})^d

L_{H_1}=¥frac{(a+c)!}{a!c!}¥frac{(b+d)!}{b!d!}(¥frac{a+c}{n})^{a+c} (¥frac{b+d}{n})^{b+d}

ふたつの大小を考えるときに、商をとると階乗が打ち消しあって捗るから、

LR=¥frac{L_{H_1}}{L_{H_0}}=¥frac{n^n a^a b^b c^c d^d}{(a+b)^{a+b} (c+d)^{c+d} (a+c)^{a+c} (b+d)^{b+d}}=¥prod_{i,j}(¥frac{n_{ij}}{e_{ij}})^{n_{ij}}

掛け算は対数で足し算になること、対数の掛け算は乗数になることを利用すれば、2倍の対数が自乗っぽくなって、カイ自乗分布に帰着できる。

帰無仮説では、片方の群の確率を決めれば、周辺度数が決まっているのでもう片方も同時に決まるので、自由度は1。一方で、対立仮説では、片方の群を決めても、ふたつの群は別々だろうという想定のもとで検定をしようとしているので、もう片方の確率は決まらない。よってこれも独自に確率を決めなければならず、自由度は2。尤度比検定をしようとするときに、比をとると、自由度が1から2に上がるので、尤度比検定をしようとするときの自由度は2-1=1 である。なので自由度1 のカイ自乗分布を採用する。

尤度比検定は、尤度が定義できて、比較したいふたつの仮説があれば、理論上はなんにでも応用可能である。理論的な部分を含むので、連続的ではあるが、サンプルサイズが小さい時のズレの問題がある。

 

3つの手法を本的にまとめたものがこちら。

検定方法計算量小さいサンプルサイズ対称性漸近近似連続・離散
正確確率保守的・正確非対称正確離散
ピアソン正確確率検定との乖離が大きい対称漸近連続
尤度比正確確率検定との乖離が大きい非対称漸近連続

 

試してみる。横軸はある2*2分割表を入力したときに、自由度1 である項を増減させることができるが、その増減をすべての場合にやり尽くした時の取りうる分割表について、y軸をp値としてプロットしている。

周辺度数が小さい分割表でやってみると、なぜかカイ自乗分布(補正あり)がもっとも大きいp値が返ってきてしまった。横軸の破線はp=0.05 を示しているが、例えば、3つ(にカイ自乗の補正ありなしか、Exact パッケージのexact.test の結果を入れている)の手法でp=0.05 をまたいだりまたがなかったりする瞬間があるわけだが、このときに、尤度比検定をすればもっともp=0.05 を下回りやすいからする、というのはやってはいけないことである。というのは、p-hacking だからである。

であれば、解析を始めるまえから、「尤度比検定を使います」と宣言して、尤度比検定しか使わないのはどうか。これは、個人的にはギリギリ許してあげようかと思う。けれども、「尤度比検定はそもそも小さなp値がでやすい。それすなわち論文になりやすい!」というような邪な考えならば、絶対に許さない絶対ニダ。

統計的な判断は、どんだけ突っ込まれようともスキのない結果を提示すべきだと思っているから、こういうときは保守的(どう頑張っても大きめなp値が出やすい検定なんだけど、それでやってもp<0.05 なんですよ、これなら文句ないでしょ?)な方法を選択しておくのが無難。

f:id:MikuHatsune:20161206163014p:image

 

周辺度数が大きな分割表では、もはやp=0.05 のようなチンケな領域においてほとんど手法間に差はない。しかし、分割表がばらついてきて(¥Delta が0 もしくは大きいところは、入力した2*2 分割表をさらに偏らせた方向に行っているので、ものすごく小さなp値が出やすい)、小さなp値が出ている。この領域では、尤度比検定がものすごい小さなp値を出す傾向にあるようである。しかし、-60 乗のp値ってそんなの誰も気にしない領域だと思うので、まあいいんじゃないかとは思うけど。

f:id:MikuHatsune:20161206163013p:image

library(Exact)
mat <- matrix(c(10, 20, 30, 40), 2) # 大きい分割表
mat <- matrix(c(3, 6, 9, 12), 2)    # 小さい分割表

# 期待値
makeExptable <- function(mat){
  m1 <- margin.table(mat, 1)
  m2 <- margin.table(mat, 2)
  N  <- sum(mat)
  etable <- m1 %*% t(m2) / N
  return(etable)
}

# 尤度比検定
LR.test <- function(mat){
  etable <- makeExptable(mat)
  chi2 <- 2 * sum(log(mat / etable) * mat)
  df <- prod(dim(mat) - 1)
  p <- pchisq(chi2, df, lower.tail=FALSE)
  return(list(stat=chi2, p.value=p))
}

# すべてのとりうる分割表
tabs <- function(mat){
  idx <- matrix(c(1, -1, -1, 1), 2)
  i <- -mat[1,1]:sum(mat)
  bs <- mapply(function(z) mat + idx*z, i, SIMPLIFY=FALSE)
  bs[ apply(sapply(bs, ">=", 0), 2, all) ]
 }

# 3つの検定法をすべてのとりうる分割表で試す
p_three <- function(mat){
  b0 <- tabs(mat)
  pv <- mapply(function(z) c(fisher=fisher.test(z)$p.value,
                             chisqF=chisq.test(z, correct=FALSE)$p.value,
                             chisqT=chisq.test(z, correct=TRUE)$p.value,
                             LR=LR.test(z)$p.value
                             ), b0)
  return(list(tables=b0, p.value=t(pv)))
}
x <- p_three(mat)

matplot(x$p.value, lty=1, type="l", log="y", xlab=expression(Delta), ylab="p value", lwd=3)
abline(h=0.05, lty=3, lwd=3)
legend("topright", legend=colnames(x$p.value), lty=1, col=seq(ncol(x$p.value)), bg="white", lwd=3)
title("Large sample size")

2016-11-27

Rmd でreveal.js のhtml スライドプレゼンテーション

Japan.R 2016 が終わりました。参加された方々はお疲れ様でした。

発表者不在というなんとも謎なLT をしたやつがいたらしいですが、Rmd でRpresentation を作って自動プレゼンするやり方について。

自動プレゼン自体はppt にも実装されています。リハーサルモードで行い、スライドめくり時間を記憶で保存すれば、勝手にスライド送りが発動します。

同時に録音ができるのでマイクで吹き込めば勝手にビデオになりますし、後から音声ファイルを追加しても可能です。

 

Rmd の出力形式をいじれば、html でスライドプレゼンテーションができるようになります。そのなかでも縦にスライド送りができるreveal.js はそこそこいいです。

Rmd のヘッダー部分にこう書いておくと個人的にはいい感じの設定になっている。

---
title: "hoge title"
subtitle: <span style="font-size:32pt">hoge subtitle</span>
author: <img src=icon96.gif style="border:none"> hoge author
date: 20161127
output:
  revealjs::revealjs_presentation:
   #self_contained: false
    transition: none 
    highlight: pygments
    #default fade slide convex concave zoom none
    reveal_options:
      #autoSlide: 5000
      width: 960
      height: 720
      control: true
      progress: true
      slideNumber: true
      previewLinks: true
      history: false
      overview: true
      center: false
      left: true
      mouseWheel: false
    fig_width: 7
    fig_height: 7
    fig_caption: true
    #mathjax: local
    #reveal_plugins: ["notes"] #search
    css: "for-revealjs.css"
# render("20161127.Rmd"); browseURL("20161127.html")
# http://rmarkdown.rstudio.com/revealjs_presentation_format.html
---
```{r eval=TRUE, echo=FALSE, warning=FALSE, message=FALSE, comment=""}
knitr::opts_chunk$set(echo=TRUE, eval=TRUE, warning=FALSE, message=FALSE, comment="", fig.height=7, fig.width=7, out.height=400, out.width=400)
knitr::knit_hooks$set(rgl=hook_webgl)
```

 

for-reveal.css として次のファイルを設定しておくとhtml に流し込まれる。

/* for {revealjs} */
/* https://github.com/hakimel/reveal.js/blob/master/css/reveal.css */
/* slide setting */
.reveal .progress {
  height: 15px;
  margin-bottom: 5px;
}

/* slide font setting */

.reveal h1,
.reveal h2,
.reveal h3,
.reveal h4,
.reveal h5,
.reveal h6 {
  text-transform: none;
}

.reveal .slide {
  text-align:left;
  font-family: segoe UI, Meiryo;
  vertical-align:top;
  /*background-color:aliceblue*/
}

/* slide-title settintg*/

.reveal section.slide>h1 {
  color: black;
  text-align:center;
  vertical-align:middle;
}

.level1>h1 {
  margin-top: 25%;
}
span.header-section-number {
  margin-right: 0.5em;
}
.reveal li {
  margin-top:0.2em;
}
.reveal ul ul li {
  font-size:0.85em;
  margin-top:0;
}

/* image setting */
.reveal img {
  /*display:block;*/
  /*text-align:center;*/
  /*max-height: 500px;*/
  top: 0;
  bottom: 0;
  left 0;
  right 0;
  margin: auto;
  border:none;
  max-width: auto;
  max-height: auto;
}

/* table setting */
.reveal table{
  border-style: solid;
  background-color:rgba(255, 255, 255, 0.8);
  width:auto;
  height:auto;
}
.reveal table th {
  border-style: solid;
  width:auto;
  height:auto;
}
.reveal table td {
  border-style: solid;
  width:auto;
  height:auto;
}

/* column setting */
.reveal .slides section .no-column {
  width: 96%;
}
.reveal .slides section .column1 {
  float: left;
  width: 48%;
}

.reveal .slides section .column2 {
  float: right;
  width: 48%;
}

/* for htmlwidgets to iframe */
.reveal .slides iframe {
	height: 450px;
	width: 100%;
	margin-left:auto;
	margin-right:auto;
}

/* slide-header setting */
#slide-header {
  text-align: right;
  background-color: #fff;
  opacity: 0.7;
}

#slide-header img {
  max-height: 50px;
  margin: 3px 10px;
  vertical-align: middle;
  opacity: 1.0;
}

#slide-header p {
  font-size: 18px;
  padding: 5px 20px 5px;
}

/* botsu */
/*
.reveal .slides section .column {
   position: fixed;
   width: 48%;
   top: 3em;
   bottom: 0;
}
.reveal .slides section .column1 {
   left:  0;
}
.reveal .slides section .column2 {
   right: 0;
}
*/

/*
.reveal .slides section .column img {
   max-width: 190%;
   max-height: 190%;
   height: auto; 
}
*/

/* other */
/*
.section .reveal .state-background {
  background: #002b36;
}
*

 

自動プレゼンのために音声ファイルを作る。推しの声優に原稿を読んでもらえるのが一番だが、声優とのコネもスタジオもなければPC マイクにちまちま吹き込む。

このときはubuntuこちらを参考に

arecord -t wav -f dat -d 3 -q | lame -b 128 -m s - out.mp3

としたあとに読み上げを開始する。

 

だらだらと吹き込んだあとはこちらを参考に、いい感じに切り出す。

ffmpeg -i input.mp4 -ss 10 -t 20 output.mp4

 

細切れのファイルがいやならまとめる

concatinate するファイル数をconcat=n= で指定する。

ffmpeg -i input1.mp4 -i input2.mp4 -i input3.mp4 -filter_complex "concat=n=3:v=1:a=1" output.mp4

 

さて、プレゼンするhtml スライドと、音声ファイルができたとしよう。

スライド送り時間は、スライドタイトルを指定した # のあとに

{data-autoslide=milliseconds}

でミリ秒を指定する。

音声をバックグランドで流したければ

{data-background-video=hoge.mp3}

で指定する。video なのでmp4 でもよい。

ここで、背景画像も指定したかったので

{data-background=hoge.png data-gackground-video=hoge.mp3 data-autoslide=1000}

と指定すると、hoge.mp3hoge.png も反映されない。

この回避のために頑張ってhoge.png静止画で表示されるmp4 を作ってもよかったけど、面倒だったのでやめた。

 

あとはmarkdwon のかんたん記法もいいけど、もうちょっと凝りたかったらやはりhtml 記法が必要になる。たとえば

改行の

</br>

文字色、サイズ、背景の

<span style="color:brack;background-color:rgba(155,155,155,0.5)">hoge text</span>

 

いいと思います(他力本願

2016-11-21

MikuHatsune2016-11-21

antimicrobial cycling

読んだ。

Proc Natl Acad Sci U S A. 2004 Sep 7;101(36):13285-90.

 

antimicrobial cycling という概念がある。細菌感染症に対して抗生物質を使うのだが、うまく治療していても耐性菌は生じることもあるし、適当に使ってしまえばさらに耐性菌のリスクがあがる。

ある抗生物質を使うことが耐性菌を生み出してしまうならば、ある一定期間抗生物質を使ったら、次の期間はごっそりと採用している抗生物質を入れ替えてしまえば、いままでいた耐性菌はいなくなるのではないか(代わりに新たに採用した抗生物質の耐性菌は出てくるけど)、という考えである。

数理モデルを使ったこの論文では、cycling による耐性菌の増加は防げない、ということになっている。

cycling そのものが効果があるかどうかはPLoS Pathog. 2014 Jun 26;10(6):e1004225. などがある。"adjustable cycling" なんて言っちゃって、結局いい感じに調整したものがいいと言っている。

cycling に対する方法として、mixing と言っている。mixing というのは、抗生物質がいくつかあって、そのうちのひとつを使って治療されるのがいいだろう、と言ってる。

 

数理モデルとして、病院内(院内感染)と病院外(市中肺炎)のふたつの領域があって、病院内に4つのコンポーネントがある。

X:uncolonized. 感染していない人。

S:colonized. 感染した人。

R:resistant. 耐性菌。drug1 とdrug2 を考えており、同時耐性は考えていない(R_1+R_2 で表される。論文中では積分)。

 

抗生物質が2種類あって、4つのコンポーネントを考えるときの各コンポーネントの流出入(微分方程式)は以下のようになる。

¥Delta S=(m-S)¥mu - (¥tau_1+¥tau_2+¥gamma)S + ¥beta S X + ¥sigma ¥beta (c_1 R_1+c_2 R_2)S

¥Delta R_1=(m_1-R_1)¥mu - (¥tau_2+¥gamma)R_1 + ¥beta(1-c_1)R_1 X + ¥sigma ¥beta (c_1 S+(c_1 - c_2) R_2)R_1

¥Delta R_2=(m_2-R_2)¥mu - (¥tau_1+¥gamma)R_2 + ¥beta(1-c_2)R_2 X + ¥sigma ¥beta (c_2 S+(c_2 - c_1) R_1)R_2

¥Delta X=(1-m-m_1-m_2-X)¥mu + (¥tau_1+¥tau_2+¥gamma)S + (¥tau_2 + ¥gamma)R_1 + (¥tau_1 + ¥gamma)R_2 - ¥beta X(S+(1-c_1)R_1+(1-c_2) R_2)

ここで、パラメータ

m:各コンポーネントの流出入係数

¥mu:患者のターンオーバー率。入院日数でいうと平均して¥frac{1}{¥mu} 日いると考える。

¥gamma:治療されないcolonized の人が¥frac{1}{¥gamma} 日はcolonized のまま。

c_i:fitness cost といって、抗生物質がないときに細菌が勝手に耐性菌になる頻度。論文中ではc=0 と書いてあるので両方0 にしてみる。

¥beta:uncolonized な人がcolonized になる比率

¥tau_i抗生物質i が使われる率。論文では¥tau_1+¥tau_2=0.5 としている。なので、drug 1 が使われるときは¥tau_2=0 にする。抗生物質を使えば、耐性でない限り、確実に細菌は死ぬとする。

¥sigma:supercolonization という、colonized された人が病院内で二次的に別の細菌でcolonized されたような感じの係数。スクリプトではdelta と空目してd になっている。

¥alpha:医者が抗生物質を処方するコンプライアンス、というのは、ある期間では抗生物質1のみ、と言っていても、医者の裁量(わがまま?)で抗生物質2を処方するかもしれない。¥alpha=1 だと100% 順守するけど、¥alpha=0 だと半々の確率でしか抗生物質1 (期間1 のときに)をきちんと処方しない。これはつまり¥alpha=0¥to 0.5, ¥alpha=1¥to 1写像し直す必要がある。

90日間隔にすると論文のFig 2 が再現できる。ある程度時間が経つと、抗菌薬のどちらかに耐性であるR_1+R_2 は0.65 くらいでプラトーに達する。

f:id:MikuHatsune:20161121194341p:image

 

2週間くらいでcycle すれば耐性割合がプラトーに達する前に耐性菌割合を減らすことができる。しかし、day=0 のときのR_1+R_2(=0.2+0.2) よりも基本的に多い割合で耐性菌が存在する。

f:id:MikuHatsune:20161121194342p:image

 

mixing にすれば、R_1+R_2 はほぼ0.4 のまま推移する。

f:id:MikuHatsune:20161121194340p:image

 

2日単位でcycle すれば、R_1+R_2 は0.4 よりもわずかに減る。しかし、2日単位でcycle するのはどう考えても現実的ではない。

f:id:MikuHatsune:20161121194343p:image

 

mixing がなぜcycle より耐性菌割合を減らせるかは、Fig 5 にあるが新しい抗生剤に出会う頻度が増えるからだとか、あとはdiscussion のところにいろいろ書いてある。

意識低い系臨床医であれば、とりあえずスルバシリンで治療を始めて、困ったらメロペンバンコアネメトロでいいんじゃない? (適当

m <- 0.7
m1 <- 0.05
m2 <- 0.05
b <- 1 # beta
c1 <- 0
c2 <- 0
g <- 0.03 # gamma
tau <- 0.5 # sum of taus
t1 <- 0.2
t2 <- 0.3
mu <- 0.1
d <- 0.25 # sigma
a <- 0.8  # alpha
switch_cycle <- 90

niter <- 350
X <- S <- R1 <- R2 <- rep(0, niter)
X[1] <- 0.35
S[1] <- 0.25
R1[1] <- 0.2
R2[1] <- 0.2

for(i in 1:(niter-1)){
  drug <- gtools::odd(i %/% switch_cycle) + 1 # drug used in time period
  t12 <- tau/2 * (1 + a*list(c(1, -1), c(-1, 1))[[drug]])
  t1 <- t12[1]; t2 <- t12[2]
  # minute deviation
  dS  <- (m-S[i])*mu - (t1+t2+g)*S[i] + b*S[i]*X[i] + d*b*(c1*R1[i]+c2*R2[i])*S[i]
  dR1 <- (m1-R1[i])*mu - (t2+g)*R1[i] + b*(1-c1)*R1[i]*X[i] - d*b*(c1*S[i]+(c1-c2)*R2[i])*R1[i]
  dR2 <- (m2-R2[i])*mu - (t1+g)*R2[i] + b*(1-c2)*R2[i]*X[i] - d*b*(c2*S[i]+(c2-c1)*R1[i])*R2[i]
  dX  <- (1-m-m1-m2-X[i])*mu + (t1+t2+g)*S[i] + (t2+g)*R1[i] + (t1+g)*R2[i] - b*X[i]*(S[i]+(1-c1)*R1[i]+(1-c2)*R2[i])
  # next step
  S[i+1]  <- S[i] + dS 
  R1[i+1] <- R1[i] + dR1
  R2[i+1] <- R2[i] + dR2
  X[i+1]  <- X[i] + dX
}

Y <- cbind(X, S, R1, R2, R1+R2)
colnames(Y) <- c("X", "S", "R1", "R2", "R1+R2")
yl <- c(0, 1)
matplot(Y, type="l", lty=1, lwd=3, ylim=yl, xlab="Day", ylab="Proportion", main=paste0("Antimicrobial cycling: ", switch_cycle, " days"))
abline(v=seq(niter %/% switch_cycle)*switch_cycle, lty=3, lwd=3)
legend("topright", legend=colnames(Y), col=seq(nrow(Y)), lty=1, lwd=3, bg="white")

2016-11-19

29cm vs 28.3cm

こんなページを見つけた。

理系大学院生の普段使いリュックの中身を見せます! | 時短師匠

メンズ大学院生の生活を支えるカバンの中身を細かく紹介! - Abstract Life

入院生活者のかばんの中身らしい。

 

私のかばんには、いつでも解析できるようにレッツノートと読み物用にxperia z3 tablet が絶対必要である。

レッツノートにはubuntu が入って窓OSデュアルブートできる環境が研究用(SZ5)、パワポワードだけ使う用(SX2)、格下げでもはや家で転がっているだけ用(SX1)の3台使っている。3台ともubuntu 環境がそろっている。

xperia z3 には論文電子書籍がはいっている。

パナソニック CF-SZ5ADCVS Lets note SZ5 シルバー

パナソニック CF-SZ5ADCVS Lets note SZ5 シルバー

ソニー Xperia Z3 Tablet Compact SGP612 ブラック

ソニー Xperia Z3 Tablet Compact SGP612 ブラック

 

移動は9割自転車、ときたま電車なので両手のあくリュックが必須。

ということでリュックが欲しいのだが、個人的な条件として、

・リュック単体が軽い

・雨に強い。欲を言えば完全防水、できれば耐水

・ペットボトルホルダーが欲しい

・出張でも遊び使いでもどちらでもいける

・1泊2日くらいならこれでいきたい(コミケ

・充電コード(iphone, microUSB)類をどこかに入れたい

・値段はそこそこ

という条件を満たした神のリュックがこちら。

 

これがどれくらい神かというと

・リュック単体で700g とあるが、マジで軽い

・リュック記事自体はそこそこ耐水。底の部分に防水カバーが入っていてそれを被せられる。が、これは上の部分がどうしても濡れるのでアレ

・両側にペットボトルホルダーあり。ゴミ袋詰めている

・黒を基調としていて、いつのまにかグレーのカラーバリエーションが出てた

・がんばれば3泊くらいいける(レッツノートが軽いというのがあるけど)

・サブポケット充実

・2年くらい前で5000円、いまだに値下がり中というアレ

という神っぷりである。

 

が、最近、レッツノートがSX 系からSZ 系になったときに、大きさが小さくなったので、ちょっと出かけるときにリュックというよりかボディバッグ的なアレで出かけたい気分になった。しかし、ボディバッグはどの店で見てもレッツノートは入らなさそう。

というわけで神リュックを売っているsumdex 神が小型かばんを売り出していた。

公称ではこのかばんが高さ29cm、レッツノートが28.3cm なので頑張れば入るッ…

 

結論から言うとチャックが締まりませんでした。

 

これは横幅がデカすぎで15インチPC を期待していそうなのでアレ。