[過去ﾛｸﾞ] 臨床統計もおもしろいですよ、その１ [無断転載禁止]©2ch.net (747ﾚｽ)
上下前次1-新
抽出解除 ﾚｽ栞
このｽﾚｯﾄﾞは過去ﾛｸﾞ倉庫に格納されています｡
次ｽﾚ検索 歴削→次ｽﾚ 栞削→次ｽﾚ 過去ﾛｸﾞﾒﾆｭｰ
23: 2017/05/06(土)17:10 ID:NK3Wh0hW(4/4) AAS
MCMCでシミュレーションしてみた
Inference for Stan model: boobs.
4 chains, each with iter=15000; warmup=5000; thin=1;
post-warmup draws per chain=10000, total post-warmup draws=40000.
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
d 0.990 0.024 4.251 -7.352 0.974 9.308 31862 1
JD 0.175 0.002 0.380 0.000 0.000 1.000 28355 1
JK 0.078 0.002 0.268 0.000 0.000 1.000 26103 1
dは胸囲の差
分布をグラフにすると
省4
57(1): 2017/06/18(日)16:18 ID:nxMdZkja(1/2) AAS
# NN(100)個中当たりSS個、N(5)個サンプルしてS個当たりからSSを推測する。
foo <-function(SS,NN=100,N=5){
YY=c(rep(1,SS),rep(0,NN-SS))
return(sum(sample(YY,N)))
}
ss=0:100
S=0
k=10^4
SS=NULL
for(i in 1:k){
省18
58: 2017/06/18(日)16:40 ID:nxMdZkja(2/2) AAS
当たり確率を連続量のpとしてstanを走らせると似たような結果が得られた。
> print(fit, pars='p',probs=c(.025,.5,.975),digits=4)
Inference for Stan model: coin5.
4 chains, each with iter=100000; warmup=50000; thin=1;
post-warmup draws per chain=50000, total post-warmup draws=200000.
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
p 0.1429 0.0005 0.123 0.0042 0.1097 0.4569 68121 1
画像ﾘﾝｸ[png]:i.imgur.com
99: innuendo ◆kCkk5BVA12 2017/09/27(水)16:37 ID:By0OihfH(3/9) AAS
巨乳女子大100人と桃尻女子大100人のどちらが誘いに応じやすいか検討してみる。
巨乳女子大では10人に声をかけて2人が誘いに応じた。
桃尻女子大では無作為順に声をかけていったら5人めが誘いに応じた。
どちらが誘いに応じやすいといえるか？
画像ﾘﾝｸ[png]:i.imgur.com
> lottery.mci(100,5)
mean mod
28.14286 21.00000
2.5% 97.5%
4 63
省6
102: innuendo ◆kCkk5BVA12 2017/09/27(水)17:59 ID:By0OihfH(6/9) AAS
画像ﾘﾝｸ[png]:i.imgur.com
平均値と信頼区間のどちらを選択するかだな。
135(1): innuendo ◆kCkk5BVA12 2017/10/09(月)15:13 ID:N8CcYMr4(1) AAS
>>127
画像ﾘﾝｸ[png]:imagizer.imageshack.com
136: 2017/10/09(月)15:39 ID:Wked2yci(4/6) AAS
>>135
画像ﾘﾝｸ[png]:i.imgur.com
141: 2017/10/10(火)21:23 ID:VF4JdnZT(1) AAS
>>129
宝くじAは1000本中300本が当たりで当たりは 100万円の賞金、
宝くじBは1000本中 10本が当たりで当たりは1000万円の賞金、
とAの期待値＝Bの期待値の3倍のとき
画像ﾘﾝｸ[png]:i.imgur.com
142: 2017/10/11(水)06:48 ID:8dgUXTG5(1/5) AAS
>>139
続きのコードが書き込めないので画像でアップ
画像ﾘﾝｸ[jpg]:imagizer.imageshack.com
155: 2017/10/15(日)07:15 ID:hnCh54BK(1) AAS
外部ﾘﾝｸ:en.wikipedia.org
を参考にオイラー関数をグラフにしてみた。
画像ﾘﾝｸ[png]:i.imgur.com
phi <- function(n){
nn=0:(n-1)
f=function(x,y) (x*y)%%n
names(nn)=paste0('n',nn)
z=outer(nn,nn,f)
i=which(gmp::gcd(1:(n-1),n)==1)
length(i)
省12
178: 2017/10/21(土)20:19 ID:TcIT7+c6(4/12) AAS
画像ﾘﾝｸ[png]:i.imgur.com
196: 2017/10/25(水)19:24 ID:nTDNW+TM(4/6) AAS
>>195
n=10 の時のPDFと９５％CI下限値0.7615958
画像ﾘﾝｸ[png]:i.imgur.com
197: 2017/10/25(水)19:24 ID:nTDNW+TM(5/6) AAS
>>195
n=10 の時のPDFと９５％CI下限値0.7615958
画像ﾘﾝｸ[png]:i.imgur.com
221: 2017/11/05(日)20:51 ID:G4ZpDCNG(8/8) AAS
Jikkan <- function(x){
re=summary(BESTout,ROPEm=c(-x,x))
re[['muDiff','%InROPE']]
}
xx=seq(0,2.5,by=0.01)
InRope=sapply(xx,Jikkan)
plot(xx,InRope,type='l',lwd=2, xlab='Difference',ylab='% Practically Equivalent')
abline(h=95, lty=3)
(Diff=uniroot(function(x,u0=95) Jikkan(x)-u0, c(1.5,2))$root)
plot(BESTout,ROPE=c(-Diff,Diff))
省1
240: 2017/11/11(土)19:25 ID:AQmX3Wf1(6/8) AAS
画像ﾘﾝｸ[png]:i.imgur.com
241: 2017/11/11(土)20:30 ID:AQmX3Wf1(7/8) AAS
Golgo13 vs Golgo14
Let's MCMC with JAGS.
I set the ROPE ( Range of Practical Equivalence ) as -0.01 to 0.01
> smryMCMC( mcmcCoda2 , compVal=NULL ,ropeDiff = c(-0.01,0.01))
Mean Median Mode HDIlow HDIhigh CompVal PcntGtCompVal
Diff 0.07384483 0.05247579 0.01267166 -0.02478474 0.2345431 0 90.18
PcntLtROPE PcntInROPE PcntGtROPE
Diff 3.07 15.66 81.27
It'll be illustrated as follows.
画像ﾘﾝｸ[png]:i.imgur.com
省18
242: 2017/11/11(土)23:28 ID:AQmX3Wf1(8/8) AAS
0/5 vs 1/20
画像ﾘﾝｸ[png]:i.imgur.com
243: 2017/11/12(日)16:15 ID:OKdyWAPj(1) AAS
posterior ∝ likelihood * prior
画像ﾘﾝｸ[png]:i.imgur.com
257: 2017/11/17(金)12:56 ID:nii6SWM6(1/3) AAS
# randam numbers by inverse culmutive density function
RandICDF <- function(ICDF,PDF,...){
U=runif(10000)
rand=ICDF(U,...)
hist(rand,freq=FALSE,breaks=30,
col=sample(colors(),2),main='')
curve(PDF(x,...),add=TRUE,lwd=2)
invisible(rand)
}
par(mfrow=c(2,2))
省5
258: 2017/11/17(金)13:24 ID:nii6SWM6(2/3) AAS
# randam numbers following PDF by Jhon von Neuman's method
vonNeumann <- function(PDF,xmin=0,xmax=1){
N=10000
ymax=max(PDF(seq(xmin,xmax,length=N+1)))
Ux=runif(N,xmin,xmax)
Uy=runif(N,0,ymax)
Rand=Ux[which(Uy<=PDF(Ux))]
hist(Rand,xlim=c(xmin,xmax),freq=FALSE,breaks=30,col=sample(colors(),2),main='')
curve(PDF,add=TRUE,lwd=2)
invisible(Rand)
省9
上下前次1-新書関写板覧索設栞歴
ｽﾚ情報 赤ﾚｽ抽出 画像ﾚｽ抽出 歴の未読ｽﾚ AAｻﾑﾈｲﾙ
ぬこの手 ぬこTOP 0.597s*