臨床統計もおもしろいですよ、その2

[過去ﾛｸﾞ] 臨床統計もおもしろいですよ、その2 (1002ﾚｽ)
上下前次 1-新
通常表示 512ﾊﾞｲﾄ分割ﾚｽ栞

このｽﾚｯﾄﾞは過去ﾛｸﾞ倉庫に格納されています｡
次ｽﾚ検索歴削→次ｽﾚ栞削→次ｽﾚ過去ﾛｸﾞﾒﾆｭｰ

887: 卵の名無しさん [sage] 2020/01/31(金) 00:59:52 ID:24QiZJ0Y

f <- function(n,prec=1000){ # Σ 1/kを既約分数表示する
if(n==1){
cat(1,'\n')
invisible(1)
}else{
GCD <- function(a,b){　 # ユークリッドの互除法
r = a%%b               # a=bq+r ⇒ a%%b=b%%rで最大公約数表示
while(r!=0){a = b ; b = r ; r = a%%b}
b }
library(Rmpfr)
one = mpfr(1, prec) # 1(one)を1000桁精度に設定
nn = 1:n            # nn : 1 2 3 ... n
nume=numeric(n)     # 分子の容器
for(i in nn) nume[i] = prod(nn[-i])*one # nnからi番目の要素を除いて乗算し精度アップ
nume = mpfr2array(nume, dim = length(nume)) # mpfr2arrayで加算を可能にする
Nume = sum(nume)　# numeの総和を計算して分子に
Deno=factorialZ(n) # 分母 n! = factorial(n*one)
gcd = GCD(Nume,Deno) # Numerator/Denominator約分するため最大公約数を計算
res=list(nume = Nume/gcd,deno=Deno/gcd,ratio=as.numeric(Nume/Deno)) # 最大公約数で除算して
　# 分数表示 give.head=FALSEでheader除去,digits.dで桁数を指定
# capture.outputで変数に取り込み
nm = capture.output(str(res$nume, give.head=FALSE,digits.d = prec))
dn = capture.output(str(res$deno, give.head=FALSE,digits.d = prec))
cat(paste0(nm,'/',dn,'\n'))
invisible(res)
}}
for(i in 1:30) f(i)

http://egg.5ch.net/test/read.cgi/hosp/1540905566/887

895: 卵の名無しさん [sage] 2020/01/31(金) 11:25:54 ID:24QiZJ0Y

# simulation model using binominal random number
rm(list=ls())
reedfrost <- function(p, I0, S0, n, greenwood=FALSE) {
S <- St <- rep(S0, n)  # St : Suscepibles @ time t, S:
I <- It <- rep(I0, n)  # It : Infected @ time t
q <- 1-p               # probability of non-transmission

time <- 0
while (sum(It)>0) { # until no new transmission
if (greenwood)
It <- rbinom(n, St, ifelse(It>0, p, 0))
else
It <- rbinom(n, St, 1-(q^It)) # how many ppl newly trannsmitted among St
St <- St-It                     # how many ppl left non-infected
I <- rbind(I, It)
S <- rbind(S, St)
time <- time+1
}
Inames=NULL
for(i in 0:(nrow(I)-1)) Inames[i+1]=paste0('I',i)
rownames(I)=Inames
Snames=NULL
for(i in 0:(nrow(I)-1)) Snames[i+1]=paste0('S',i)
rownames(S)=Snames
list(I0=I0,S0=S0,p=p,n=n,I=I,S=S)
}
re=reedfrost(0.05,1,99,100)
nr=nrow(re$I)
plot(0:(nr-1),re$I[,1],bty='l',type='l',lwd=2,
ylim=c(0,max(re$I)),xlab='time',ylab='Infecteds')
for(i in 2:n) lines(0:(nr-1),re$I[,i],lwd=2,col=sample(colours()))

http://egg.5ch.net/test/read.cgi/hosp/1540905566/895

898: 卵の名無しさん [sage] 2020/01/31(金) 15:00:45 ID:24QiZJ0Y

# Compute Ro - Reproductive number.
Ro = beta_value / gamma_value

# Disease dynamics parameters.
parameter_list = c (beta = beta_value, gamma = gamma_value, delta = delta_value)

# Compute total population.
N = s + i + r + e

# Initial state values for the differential equations.
initial_values = c (S = s/N, E = e/N, I = i/N, R = r/N)

# Simulate the SEIR epidemic.
# ?lsoda # Solver for Ordinary Differential Equations (ODE), Switching Automatically Between Stiff and Non-stiff Methods
output = lsoda (initial_values, timepoints, seir_model, parameter_list)

# Plot dynamics of Susceptibles sub-population.
plot  (S ~ time, data = output, pch='S', type='b', col = 'blue',  bty='l', cex=0.75,
ylab='SEIR')
points(E ~ time, data = output, pch='E', type='b', col = 'yellow',  cex=0.75)
points(I ~ time, data = output, pch='I', type='b', col = 'red',   cex=0.75)
points(R ~ time, data = output, pch='R', type='b', col = 'green', cex=0.75)
}

SEIR()
args(SEIR)
SEIR(contact_rate=0.5,transmission_probability=0.05,infectious_period=365,latent_period=7,s=99,e=0,i=1,r=0,timepoints = 0:365*6)

http://egg.5ch.net/test/read.cgi/hosp/1540905566/898

904: 卵の名無しさん [sage] 2020/02/04(火) 16:35:57 ID:b64IHQrg

"https://this.kiji.is/597220008571339873?c=39550187727945729
中国湖北省武漢市からチャーター機で日本へ帰国した邦人の新型コロナウイルス感染率が高いと、
中国で驚きの声が上がっている。中国当局が発表した同市の感染者の割合に比べ「39倍も高い」というのだ。
現地は医療現場が混乱しているため、実際には発表よりかなり多くの感染者がいる可能性がある。
日本政府はチャーター機計3便を武漢市に派遣し、邦人565人が帰国した。厚生労働省によると、
チャーター機に乗っていた感染者は、症状のない人も含め計8人。感染率は1.416％だ。
一方、1月31日現在、武漢市の感染者数は3215人で、感染率は0.036％にとどまった。"
r1=8
r2=3215
n1=565
n2=round(3215/(0.036/100))
poisson.test(c(r1,r2),c(n1,n2))
prop.test(c(r1,r2),c(n1,n2))
library(BayesFactor)
mat=cbind(c(r1,r2),c(n1,n2)-c(r1,r2))
fisher.test(mat)
chisq.test(mat)$p.value
bf=BayesFactor::contingencyTableBF(mat,sampleType = 'poisson',posterior = TRUE,iteration=1e5)
rbf=round(bf)
(r1/n1)/(r2/n2)
frp <- function(x) (x[1]/(x[1]+x[2]))/ (x[3]/(x[3]+x[4]))
rrp=apply(rbf,1,frp)
hist(rrp)
BEST::plotPost(rrp)
(r1/(n1-r1))/(r2/(n2-r2))
fop <- function(x) (x[1]/x[2])/ (x[3]/x[4])
orp=apply(rbf,1,fop)
BEST::plotPost(orp)
quantile(orp,c(0.005,0.025,0.5,0.975,0.995))
library(exact2x2)
exact2x2::fisher.exact(mat)
fisher.test(mat)

http://egg.5ch.net/test/read.cgi/hosp/1540905566/904

上下前次 1-新書関写板覧索設栞歴

あと 98 ﾚｽあります
ｽﾚ情報赤ﾚｽ抽出画像ﾚｽ抽出歴の未読ｽﾚ AAｻﾑﾈｲﾙ

ぬこの手ぬこTOP 0.151s*