[過去ﾛｸﾞ] 臨床統計もおもしろいですよ、その１ [無断転載禁止]©2ch.net (747ﾚｽ)
上下前次1-新
このｽﾚｯﾄﾞは過去ﾛｸﾞ倉庫に格納されています｡
次ｽﾚ検索 歴削→次ｽﾚ 栞削→次ｽﾚ 過去ﾛｸﾞﾒﾆｭｰ
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
237: 2017/11/11(土)16:51 ID:AQmX3Wf1(3/8) AAS
#-----------------------------------------------------------------------------
# RUN THE CHAINS
parameters = c( "theta") # The parameters to be monitored
adaptSteps = 500 # Number of steps to adapt the samplers
burnInSteps = 500 # Number of steps to burn-in the chains
nChains = 4 # nChains should be 2 or more for diagnostics
thinSteps = 1
nIter = ceiling( ( numSavedSteps * thinSteps ) / nChains )
# Create, initialize, and adapt the model:
jagsModel = jags.model( "TEMPmodel.txt" , data=dataList , inits=initsList ,
省15
238: 2017/11/11(土)19:06 ID:AQmX3Wf1(4/8) AAS
library("rstan")
library("rstudioapi")
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
G1314model='
data{// Golgo13.14.stan
int N13; // 100
int N14; //10
int<lower=0,upper=1> G13[N13];
int<lower=1,upper=2> G14[N14];
省16
239: 2017/11/11(土)19:15 ID:AQmX3Wf1(5/8) AAS
N13=100
N14=10
G13=rep(1,N13)
G14=rep(1,N14)
data=list(N13=N13,N14=N14,G13=G13,G14=G14)
G1314.model=stan_model('G1314.stan')
saveRDS(G1314.model,file='G1314.model.rds')
# G1314.model=readRDS('G1314.model.rds')
fitG1314 <- sampling(G1314.model, data=data, seed=1234)
fitG1314
省10
240: 2017/11/11(土)19:25 ID:AQmX3Wf1(6/8) AAS
画像ﾘﾝｸ[png]:i.imgur.com
上下前次1-新書関写板覧索設栞歴
あと 507 ﾚｽあります
ｽﾚ情報 赤ﾚｽ抽出 画像ﾚｽ抽出 歴の未読ｽﾚ AAｻﾑﾈｲﾙ
ぬこの手 ぬこTOP 0.375s*