ゴーヤのために支柱を立てて、ネットを張った。
あと、枝豆の苗とパプリカを植えて、またまたトマトの苗を植えた。娘が、大根を植えると言ってゆずらないので二十日大根の種を買ってまいたのと、つるなしインゲンもついでにまいてみた。
なかなか畑らしくなった。
macbookにgraphviz,pygraphviz,networkxを入れた
コンパートメントモデルを書くのにいちいちInkspaceを立ち上げるのもめんどうなのでGraphvizにでも頼るかと思ったらmacbookには入ってなかった。networkxもついでに入れといた。
sudo port install graphviz
sudo easy_install-2.6 pygraphviz
sudo easy_install-2.6 networkx
使う
>>> import pygraphviz as pgv
>>> d={'1': {'2': None}, '2': {'1': None, '3': None}, '3': {'2': None}}
>>> A=pgv.AGraph(d)
>>> A.layout(prog='dot')
>>> A.draw('pyg.png')
PK1-3(Digoxinの血漿中濃度推移)
設問1
フィッティングのための初期値を求める。
t <- c(2.1, 4.1, 6.1, 8.1, 10.1, 14.1, 18.1, 22.1, 30.1, 45.2, 60.2,
120.3, 180.5, 241.8, 470, 948, 1423, 2187, 2831, 4305)
conc <- c(20.5, 17.5, 14.5, 12.5, 13.0, 12.0, 11.0, 9.1, 9.6, 5.6, 4.9,
3.2, 2.0, 1.8, 0.90, 0.85, 0.70, 0.45, 0.43, 0.39)
logc <- log(conc)
lm(logc[15:20] ~ t[15:20])
Call:
lm(formula = logc[15:20] ~ t[15:20])
Coefficients:
(Intercept) t[15:20]
-0.0399743 -0.0002441
exp(-0.0399)
[1] 0.9608855
beta <- function(t){exp(-0.0399743-t)}
ac <- conc[1:14] - beta(conc[1:14])
lm(ac ~ t[1:14])
Call:
lm(formula = ac ~ t[1:14])
Coefficients:
(Intercept) t[1:14]
13.19420 -0.06273
設問2
2,3exp式でフィッティングし、あてはめの良さをくらべる。
e2 <- nls(conc ~ a*exp(-alpha*t)+b*exp(-beta*t),
start=list(a=16,alpha=0.028,b=2,beta=0.000637))
e3 <- nls(conc ~ a*exp(-alpha*t)+b*exp(-beta*t)+c*exp(-gamma*t),
start=list(a=16,alpha=0.13,b=8.4,beta=0.0114,c=0.96,gamma=0.000239))
あてはめのよさは赤池の情報量基準で。Rの場合はAIC関数でよい。
AIC(e2) # [1] 56.79388
AIC(e3) # [1] 44.93942
というわけで、exp3式のほうがあてはまりがよい。
設問3
skip
ファーマコキネティクス―演習による理解macbookにGAMESS
DMPKブームがきた！とか思っていたのだけど、突然量子化学計算に引き戻された。
macbookにGAMESSを入れる。とりあえず64bitのを入れた。
溶媒効果を取り入れた計算がしたかったので、exampleを見てたらPCMのサンプルがあったので、これをながしてみる(exam31)
$contrl scftyp=rhf runtyp=optimize nzvar=12 $end
$system timlim=2 mwords=2 $end
$pcm solvnt=water $end
$basis gbasis=n31 ngauss=6 ndfunc=1 $end
$guess guess=huckel $end
$zmat izmat(1)=1,1,2, 1,2,3, 1,3,4, 1,3,5, 1,3,6,
2,1,2,3, 2,2,3,4, 2,2,3,5, 2,2,3,6,
3,1,2,3,4, 3,1,2,3,5, 3,1,2,3,6 $end
$statpt opttol=1d-5 $end
$data
Methanol in PCM water...starting at gas phase geom
Cs
H 1.0 -1.0616171503 0.8036449245 0.0000000000
O 8.0 -0.6870131482 -0.0653470836 0.0000000000
C 6.0 0.7093551399 0.0291827007 0.0000000000
H 1.0 1.0836641283 0.5408321444 0.8835398105
H 1.0 1.0975386849 -0.9797829903 0.0000000000
$end
実行結果(抜粋)
FREE ENERGY IN SOLVENT = -72190.31 KCAL/MOL
INTERNAL ENERGY IN SOLVENT = -72185.33 KCAL/MOL
DELTA INTERNAL ENERGY = 0.00 KCAL/MOL
ELECTROSTATIC INTERACTION = -4.97 KCAL/MOL
PIEROTTI CAVITATION ENERGY = 0.00 KCAL/MOL
DISPERSION FREE ENERGY = 0.00 KCAL/MOL
REPULSION FREE ENERGY = 0.00 KCAL/MOL
TOTAL INTERACTION = -4.97 KCAL/MOL
TOTAL FREE ENERGY IN SOLVENT = -72190.31 KCAL/MOL
さすがに一瞬で終わる。実用を考えるとPerlMolとかopenbabelのオブジェクトからインプットを吐き出して計算させてその結果をオブジェクトにして返すモジュールがあれば素晴らしいのだが。
すぐできる 量子化学計算ビギナーズマニュアル (KS化学専門書)量子化学専門のヒトは長いインプットから長いアウトプットを出して、それをなめまわす傾向はあるよな(偏見？)
普通のプログラムだとデバッグオプションとかverbose的な扱いの冗長ログが標準というか、、、、
丹念に見るって言うのかなぁ。「必要に応じて覗く」じゃなくて「常に完全に把握する」みたいな感じか。
22052009 chemoinformatics
茄子の三本仕立て
本に茄子は花が咲いたら上から二つを残してわき芽をかいて三本仕立てにすべしとあったので、今朝早起きして早速かきにいった。
加茂茄子だけどな。
現代にいきづく京の伝統野菜 古都の食文化を担って20052009 life
いまさらながら「金持ち父さん」
文章が読みにくい。U隊長は苦行だといってた。
金持ち父さん貧乏父さん投資のすすめ。リターンを生むようなものをこつこつ貯める。そんな感じ。 いずれにせよ、手を動かさなくてもある程度一定の見返りが得られるような 仕組みは必要だろう。
それを、自給という道に求めれば半農半Xかな。
半農半Xという生き方 (ソニー・マガジンズ新書)レアなキャリアを持てば、会社に対してもある程度通しやすいものがあるし、しょうもない仕事を嫌々やる必要もないし。
結局、何かへ極度に依存することを避け、選択肢を広げておくために色々工夫をしろということかな。
19052009 life
「お金の学校」を読んだ
TOEICのアンケートに答えたら図書カードがもらえたので。
勝間和代のお金の学校―サブプライムに負けない金融リテラシーamazonでの評価はあまりよろしくないようだが、僕は結構面白かった。考えるための材料が対談のあちこちに散らばっていたので。
1章の竹中平蔵さんと3章の太田忠さんは普通に面白かったが、4章の「金融から未来を変える」はかなりよかった。
4章読み始めた時はなんでこんなしょうもないのをここにもってくるのかな？と不思議だったが、読み進めていくにしたがい、あーそういう風に考えるのかと思った。ただ一方で、その目指している方向を実現するには、もっと小規模なというか規模を追求しないスタンスの会社でないとサステイナブルな社会というのはなかなか難しいのではないかなぁと思った。
副題に金融リテラシーとか書いてあるけど、これってどっちかというと生き方指南だよなぁと思ったが、よくよく考えてみるにこれはトートロジーではないかと妙に納得した。
ので、酒飲んで寝よーっと。
PK1-2(コンパートメントモデル)
微分方程式を解くためにodesolveを利用する。初期値の与えかたとか、ぐるぐる回していくイメージをつかむのに微妙に時間がかかったが、覚えてしまうと簡単な方程式を解かせたりシミュレーションするのには楽になる。
設問1
吸収過程を考えないコンパートメントモデル
library(odesolve)
params <- c(kel = 0.037)
times <- c(0,(1:40))
dydt <- function(t,y,p){
kel <- p['kel']
y <- -kel*y
list(c(y))
}
y <- lsoda(c(y = 8.6),times,dydt,params)
plot(y,ylab="Conc",xlab="time (hour)")
設問2&3
今度は吸収過程を考える。設問3は吸収速度定数を1/10にした時にどういう挙動を示すか。
params <- c(kel = 0.037,ka = 0.1, D = 5000, V = 580)
times <- c(0,(1:120))
dydt <- function(t,y,p){
kel <- p['kel']
ka <- p['ka']
D <- p['D']
V <- p['V']
y <- ka * (D/V) * exp(-ka*t) -kel*y
list(c(y))
}
y <- lsoda(c(y = 0),times,dydt,params)
params2 <- c(kel = 0.037,ka = 0.01, D = 5000, V = 580)
y2 <- lsoda(c(y = 0),times,dydt,params2)
plot(y,ylim=c(0,6),ylab="Conc",xlab="time(hour)")
par(new=T)
plot(y2,ylim=c(0,6),axes=F,xlab="",ylab="")
今回バイオアバイラビリティ(F)を1とおいているが、これを変化させた時の挙動も気になるのでやってみた。設問の式にFを組み込んでF=0.1とした時の結果
dydt2 <- function(t,y,p){
kel <- p['kel']
ka <- p['ka']
D <- p['D']
V <- p['V']
F <- p['F']
y <- ka * (D/V) * F * exp(-ka*t) -kel*y
list(c(y))
}
params3 <- c(kel = 0.037,ka = 0.1, D = 5000, V = 580, F = 0.1)
y3 <- lsoda(c(y = 0),times,dydt2,params3)
plot(y2,ylim=c(0,3),axes=F,xlab="time(hour)",ylab="Conc")
par(new=T)
plot(y3,ylim=c(0,3),axes=F,xlab="",ylab="")
あーなるほど、そうだよなぁと一人で妙に納得したのであった。
