[過去ﾛｸﾞ] 臨床統計もおもしろいですよ、その１ [無断転載禁止]©2ch.net (747ﾚｽ)
上下前次1-新
このｽﾚｯﾄﾞは過去ﾛｸﾞ倉庫に格納されています｡
次ｽﾚ検索 歴削→次ｽﾚ 栞削→次ｽﾚ 過去ﾛｸﾞﾒﾆｭｰ
217: 2017/11/05(日)20:18 ID:G4ZpDCNG(4/8) AAS
Bayesian Estimation Supersedes the t-test (BEST) - online
外部ﾘﾝｸ:www.sumsar.net
トレースプロットがみられるのがうれしい。
218: 2017/11/05(日)20:22 ID:G4ZpDCNG(5/8) AAS
新版の予約受付中か。
外部ﾘﾝｸ:www.amazon.co.jp
ベイズ統計がどれくらい組み込まれたをみてから買うかな。
219: 2017/11/05(日)20:28 ID:G4ZpDCNG(6/8) AAS
# A君の彼女は女子大生、B君の彼女は女子高生。
# Y1女子大生n1=100人とY2女子高生n2=100人の胸囲を測定して
# 前者が平均82 , 標準偏差3
# 後者が平均81 , 標準偏差3
# 有意差はあるか？
T.test=function(n1,n2,m1,m2,sd1,sd2){
SE12=sqrt((1/n1+1/n2)*((n1-1)*sd1^2+(n2-1)*sd2^2)/((n1-1)+(n2-1)))
T=(m1-m2)/SE12
2*pt(abs(T),n1-1+n2-1,lower.tail = FALSE)
}
省17
220: 2017/11/05(日)20:33 ID:G4ZpDCNG(7/8) AAS
ProData <- makeData(mu1=82, sd1=3, mu2=81, sd2=3, nPerGrp=100,
pcntOut=10, sdOutMult=2.0, rnd.seed=NULL,showPlot=TRUE)
proMCMC <- BESTmcmc(proData$y1, proData$y2, numSavedSteps=2000)
N1plan <- N2plan <- 100
powerPro <- BESTpower(proMCMC, N1=N1plan, N2=N2plan,
ROPEm=c(-2,2), ROPEsd=c(-0,0), ROPEeff=c(-0,0),
maxHDIWm=5.0, maxHDIWsd=2.0, maxHDIWeff=1.0, nRep=1000)
powerPro
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
222: 2017/11/06(月)21:07 ID:1uZMaFLW(1/2) AAS
##
library(BayesFactor)
tBF <- function(x,y){
t=t.test(x,y,var.equal = TRUE)$statistic
ttest.tstat(t,length(x),length(y),simple=TRUE)
}
NHST2 <- function(N){
x=rnorm(3)
y=rnorm(3)
while(length(x) < N & tBF(x,y) <= 1){
省18
223: 2017/11/06(月)21:51 ID:1uZMaFLW(2/2) AAS
>>216
Sequential Samplingに修正
library(BEST)
BA <- function(x,y){
mc=BEST::BESTmcmc(x,y,num=5000)
re=summary(mc,ROPEm=c(-0.1,0.1))
return(re[['muDiff','%InROPE']])
}
x=rnorm(3)
y=rnorm(3)
省8
224: 2017/11/07(火)00:05 ID:2r/5Mjhu(1/4) AAS
>>201
n=1;r=1
f=function(p)choose(n,r)*p^r*(1-p)^(n-r)
AUC=integrate(f,0,1)$value
integrate(function(x) x*f(x)/AUC,0,1)$value
225: 2017/11/07(火)00:49 ID:2r/5Mjhu(2/4) AAS
library(BayesFactor)
tBF2 <- function(x,y,...){
tt=t.test(x,y,var.equal = TRUE)
bf=ttest.tstat(tt$statistic,length(x),length(y),simple=TRUE,...)
p=tt$p.value
return(c(bf=1/bf,p=p))
}
N=10
BF=numeric()
p.value=numeric()
省12
226: 2017/11/07(火)20:40 ID:2r/5Mjhu(3/4) AAS
# crooked coin or dice
crooked <- function(n,r,H0=0.5,d=0.1,cl=0.95,Print=TRUE,...){
hdi=HDInterval::hdi(qbeta,cl,shape1=1+r,shape2=1+n-r)
if(Print){
ROPEl=H0*(1-d)
ROPEu=H0*(1+d)
curve(dbeta(x,1+r,1+n-r),lwd=1,xlab='p',ylab='density',...) # posterior
segments(ROPEl,0,ROPEu,lwd=5,col='navy')
segments(hdi[1],0,hdi[1],dbeta(hdi[1],1+r,1+n-r),col='blue')
segments(hdi[2],0,hdi[2],dbeta(hdi[2],1+r,1+n-r),col='blue')
省8
227: 2017/11/07(火)21:15 ID:2r/5Mjhu(4/4) AAS
# crooked coin or dice
crooked <- function(n,r,H0=0.5,d=0.1,credMass=0.95,Print=TRUE,...){
hdi=HDInterval::hdi(qbeta,credMass,shape1=1+r,shape2=1+n-r)
if(Print){
ROPEl=H0*(1-d)
ROPEu=H0*(1+d)
curve(dbeta(x,1+r,1+n-r),lwd=1,xlab='p',ylab='density',...) # posterior
segments(ROPEl,0,ROPEu,lwd=5,col='navy')
segments(hdi[1],0,hdi[1],dbeta(hdi[1],1+r,1+n-r),col='blue')
segments(hdi[2],0,hdi[2],dbeta(hdi[2],1+r,1+n-r),col='blue')
省6
228: 2017/11/08(水)19:29 ID:vjtLqicz(1/4) AAS
data{// coin5d.stan
int N;
int<lower=0,upper=1> Y[N];
real b; // border
}
parameters{
real<lower=0,upper=1> p;
}
transformed parameters{
real ura = step(b-p); // ifelse(b-p>0,1,0)
省5
229: 2017/11/08(水)19:36 ID:vjtLqicz(2/4) AAS
N=10
r=5
Y=c(rep(1,r),rep(0,N-r))
b=0.6
data <- list(Y=Y, N=N,b=b)
model.coin5d <- stan_model('coin5d.stan')
fit.coin5 <- sampling(model.coin5d, data=data, seed=123)
print(fit.coin5,digits=4)
Inference for Stan model: coin5d.
4 chains, each with iter=2000; warmup=1000; thin=1;
省5
230: 2017/11/08(水)20:57 ID:vjtLqicz(3/4) AAS
# NHST Has 100% False Alarm Rate in Sequential Testing
# NHST : Null Hypothesis Significance Testing
# Under NHST, sequential testing of data generated from the null
# hypothesis will eventually lead to a false alarm. With infinite
# patience, there is 100% probability of falsely rejecting the null.
# This is known as “sampling to reach a foregone conclusion” (e.g.,
# Anscombe, 1954). To illustrate this phenomenon,
# a computer simulation generated random values from a normal distribution
# with mean zero and standard deviation one, assigning each sequential
# value alternately to one or the other of two groups, and at each
省4
231: 2017/11/08(水)20:57 ID:vjtLqicz(4/4) AAS
FUN = wilcox.test
FUN = perm::permTS
FUN = lawstat::brunner.munzel.test
NHST <- function(N){
x=rnorm(3)
y=rnorm(3)
while(length(x) < N & FUN(x,y)$p.value >= 0.05){
x=append(x,rnorm(1))
y=append(y,rnorm(1))
}
省15
232: 2017/11/11(土)09:46 ID:EghxJpRr(1/2) AAS
source('外部ﾘﾝｸ:github.com
source('外部ﾘﾝｸ:github.com
a=1;b=1
N=21
z=21
pp=seq(0,1,by=0.001)
likelihood=pp^z*(1-pp)^(N-z)
prior=dbeta(pp,a,b)
posterior = BernGrid(pp,prior/sum(prior),c(rep(1,z),rep(0,N-z)), plotType="Bars" ,
showCentTend="Mean" , showHDI=TRUE , showpD=TRUE )
233: 2017/11/11(土)10:09 ID:EghxJpRr(2/2) AAS
# core concept
par(mfrow=c(2,2))
z=1 ; N=1
pp=seq(0,1,by=0.001)
prior=pmin(pp,1-pp) # a isosceles triangular distribution
plot(pp,prior,type='h',lwd=5,col='lightblue')
likelihood=pp^z*(1-pp)^(N-z)
plot(pp,likelihood,type='h',col='lightblue')
(pD=sum(prior*likelihood))
posterior=prior*likelihood/pD
省1
234: 2017/11/11(土)11:35 ID:ARZ4w8db(1) AAS
require(rjags)
z=21
N=21
y=c(rep(1,z),rep(0,N-z))
data=list(N=N,y=y)
modelString = "
model {
for (i in 1:N) {
y[i] ~ dbern(theta)
}
省11
235: 2017/11/11(土)16:50 ID:AQmX3Wf1(1/8) AAS
genMCMC <- # prior : theta ~ Beta(a,b)
function( data , numSavedSteps=50000 , saveName=NULL, a=1,b=1) {
require(rjags)
y = data$y
s = as.numeric(data$s) # converts character to consecutive integer levels
# Do some checking that data make sense:
if ( any( y!=0 & y!=1 ) ) { stop("All y values must be 0 or 1.") }
Ntotal = length(y)
Nsubj = length(unique(s))
# Specify the data in a list, for later shipment to JAGS:
省20
236: 2017/11/11(土)16:50 ID:AQmX3Wf1(2/8) AAS
#-----------------------------------------------------------------------------
# INTIALIZE THE CHAINS.
# Initial values of MCMC chains based on data:
# Option 1: Use single initial value for all chains:
# thetaInit = rep(0,Nsubj)
# for ( sIdx in 1:Nsubj ) { # for each subject
# includeRows = ( s == sIdx ) # identify rows of this subject
# yThisSubj = y[includeRows] # extract data of this subject
# thetaInit[sIdx] = sum(yThisSubj)/length(yThisSubj) # proportion
# }
省13
上下前次1-新書関写板覧索設栞歴
あと 511 ﾚｽあります
ｽﾚ情報 赤ﾚｽ抽出 画像ﾚｽ抽出 歴の未読ｽﾚ AAｻﾑﾈｲﾙ
ぬこの手 ぬこTOP 0.269s*