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

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

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

884: 卵の名無しさん [sage] 2020/01/30(木) 21:21:54 ID:m/EdOA9B

Rはデフォルトでは不定長さ整数を扱えないので工夫がいる。
その点はpythonやHaskellの方が有利。

"
H(n) = Σ[k=1,2,...,n] 1/k
とする。H(n)を既約分数で表したときの分子の整数をf(n)と表す。
(1）lim[n→∞] H(n) を求めよ、答えのみで良い。
(2）n=1,2,...に対して、f(n)に現れる1桁の整数を全て求めよ
"
rm(list=ls())

H <- function(n,prec=1000){ # Σ 1/kを既約分数表示する
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 = factorial(n)*one # 分母のn!を精度アップ
gcd = GCD(Nume,Deno) # Numerator/Denominator約分するため最大公約数を計算
res=list(nume = Nume/gcd,deno=Deno/gcd,ratio=as.numeric(Nume/Deno)) # 最大公約数で除算して
cat(as.numeric(res$nume),'/',as.numeric(res$deno),'\n') # 数値化して分数表示
invisible(res)
}

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

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

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

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

ぬこの手ぬこTOP 0.007s