[過去ﾛｸﾞ] 臨床統計もおもしろいですよ、その2 (1002ﾚｽ)
上下前次1-新
抽出解除 必死ﾁｪｯｶｰ(本家) (べ) ﾚｽ栞 あぼーん
このｽﾚｯﾄﾞは過去ﾛｸﾞ倉庫に格納されています｡
次ｽﾚ検索 歴削→次ｽﾚ 栞削→次ｽﾚ 過去ﾛｸﾞﾒﾆｭｰ
388: 2018/11/30(金)10:07 ID:c4eruZjZ(1/8) AAS
rm(list=ls())
graphics.off()
par(mfrow=c(1,2))
a=360 ; b=1
R = function(t) ifelse(0<=t&t<=2*b,-a*t*(t-2*b),0)
N = integrate(R,0,2*b)$value ; (N=4*a*b^3/3)
A = function(t) ifelse(0<=t&t<=2*b,-a*t^3/3 +a*b*t^2,N)
curve(A(x),0,3,lwd=2,bty='l',xlab='t')
mu=100
n.win=2
省20
389: 2018/11/30(金)10:11 ID:c4eruZjZ(2/8) AAS
c2Wq <- function(c,a=360,b=1){ #-> Wq:平均待ち時間
# R(t): 到着率関数 -at(t-2b)
# c:サービス率
R = function(t) ifelse(0<=t&t<=2*b,-a*t*(t-2*b),0)
N=4*a*b^3/3
d = sqrt(b^2-c/a)
min(case1=9/4*a*d^4/N,case2=a/36*(b+d)^3/(b-d)*(4*b*d-b^2-d^2)/N)
}
c2Wq(300)
390: 2018/11/30(金)14:10 ID:c4eruZjZ(3/8) AAS
rm(list=ls())
graphics.off()
par(mfrow=c(2,1))
a=360 ; b=1 # R(t) at(t-2b) 到着率関数[0,2b]
R = function(t) ifelse(0<=t&t<=2*b,-a*t*(t-2*b),0)
curve(R(x),0,3,bty='l',xlab='t')
N = integrate(R,0,2*b)$value ; 4*a*b^3/3 # 総人数
A = function(t) ifelse(t<=2*b,-a*t^3/3 +a*b*t^2,N) # 流入関数=∫Rdt
# = integerate(function(t) R(t),0,t) 0<t<2b
curve(A(x),0,3,bty='l',xlab='t')
省20
391: 2018/11/30(金)14:12 ID:c4eruZjZ(4/8) AAS
Q.<- function(t){
if(t1<=t & t<= min(t4,2*b)) -a*(t-t1)^3/3 + a*d*(t-t1)^2}
# -a*(t-t1)^3/3 + a*d*(t-t1)^2 # 解析値[t1,min(t4,2b)]
tt=seq(t1,min(t4,2*b),len=1000)
lines(tt,sapply(tt,Q.),bty='l',type='l',xlab='t',ylab='Wq')
par(mfrow=c(1,1))
curve(A(x),0,3,bty='l',ylab='person',xlab='t',lwd=1,
main="到着率:at(t-2b) サービス率:200") # 累積入場者
curve(Q(x),0,3,col='navy',add=T,lwd=1,lty=3,type='h') # 待ち人数
tt=seq(0,3,len=1000)
省6
392: 2018/11/30(金)14:22 ID:c4eruZjZ(5/8) AAS
数式を追うだけだと身につかないからプログラムに入力して自分でグラフを書いてみると理解が捗る。
自分がどこができていないもよく分かる。必要な計算ができないとグラフが完成できないから。
プログラムしておくとあとで数値を変えて再利用できるのが( ・∀・)ｲｲ!!
画像ﾘﾝｸ[png]:i.imgur.com
394: 2018/11/30(金)23:12 ID:c4eruZjZ(6/8) AAS
#
source('tmp.tools.R') # 乱数発生にNeumann法
# 受付時間9:00-12:30,15:30-19:00
curve(10*(dgamma(x-9,2,1)+dgamma(x-16,8,5)),9,20,type='h') # 雛形
R <- function(x) ifelse((9<x&x<12.5)|(15.5<x&x<19),dgamma(x-9,2,1)+dgamma(x-16,8,5),0)
set.seed(123) ; data=vonNeumann2(function(x) R(x),9,19,Print=F)
N=100 # 来院患者数
n.win=1 # サービス窓口数
mu=8 # サービス率(１時間診察人数)
client=hist(data,freq=F,breaks=30,col='skyblue',main='',xlab='clock time')
省10
395: 2018/11/30(金)23:12 ID:c4eruZjZ(7/8) AAS
tt=seq(9,24,len=1000)
Rtt=R(tt)
plot(tt,Rtt,type='s',bty='l')
cumR=cumsum(Rtt)/sum(Rtt)*N # cumsumで累積来院数をgrid化
plot(tt,cumR,type='l',bty='l')
A <- function(t) cumR[which.max(t<=tt)] # 離散量を連続関数に
A=Vectorize(A)
curve(A(x),9,24,bty='l')
Q <- function(t){ # 時刻tでの待ち人数
if(t<t1) return(0)
省13
396: 2018/11/30(金)23:51 ID:c4eruZjZ(8/8) AAS
午前の受付９時から１２時３０分まで午後の受付１５時３０分から１９時までのクリニックに
図のような二峰性の分布で１００人が来院するとする。
画像ﾘﾝｸ[png]:i.imgur.com
> breaks
[1] 9.0 9.5 10.0 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5
[13] 15.0 15.5 16.0 16.5 17.0 17.5 18.0 18.5 19.0
> round(y)
[1] 5 8 10 8 6 5 3 0 0 0 0 0 0 0 0 7 19 17 8 3
医師は一人、診察時間は平均8分として待ち時間をグラフ化。
画像ﾘﾝｸ[png]:i.imgur.com
省5
上下前次1-新書関写板覧索設栞歴
ｽﾚ情報 赤ﾚｽ抽出 画像ﾚｽ抽出 歴の未読ｽﾚ AAｻﾑﾈｲﾙ
ぬこの手 ぬこTOP 2.148s*