[過去ﾛｸﾞ] 臨床統計もおもしろいですよ、その１ [無断転載禁止]©2ch.net (747ﾚｽ)
上下前次1-新
このｽﾚｯﾄﾞは過去ﾛｸﾞ倉庫に格納されています｡
次ｽﾚ検索 歴削→次ｽﾚ 栞削→次ｽﾚ 過去ﾛｸﾞﾒﾆｭｰ
256: 2017/11/15(水)14:20 ID:RUmagzE5(1) AAS
In Bayesian data analysis, evidence is the marginal likelihood （Integrate P(D|Θ)P(Θ)dΘ） which MCMC cannot yield.
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
259: 2017/11/17(金)16:06 ID:nii6SWM6(3/3) AAS
Golgo Script will be reduced as follows:
N shoots with z hits
Nz <- function(N,z,...){
curve(dbeta(x,1+z,1+N-z),xlab='Probability of Hit',
ylab=paste('Probabilty of',z,'Hits out of ',N,'Shoots'),...)
hdi=HDInterval::hdi(qbeta,shape1=1+z,shape2=1+N-z)
print(hdi,digits=4)
}
260: 2017/11/18(土)10:52 ID:V/eIhsOZ(1) AAS
modelString =
'model {
z ~ dbin(0.5,N)
N ~ dpois(lambda)
p=z/N
}
'
writeLines(modelString , con="TEMPmodel.txt" )
f <- function (lambda){
dataList=list(lambda=lambda)
省10
261: 2017/11/18(土)16:59 ID:FiDbN6uZ(1) AAS
graphics.off()
p=binom.test(7,24,alt='less')$p.value
1-(1-p)^2 # 0.06289339
z1=7
N1=24
N2=24
plot(0:N1/N1,dbinom(0:N1,N1,0.5),type='h',lwd=5,col='skyblue',ann=FALSE,ylim=c(0,0.25))
points(z1/N1,0,pch='+',cex=2)
binom.test(z1,N1,alt='less')$p.value
0:N1/N1 < z1/N1
省14
262: 2017/11/20(月)19:17 ID:xJug4kDO(1/4) AAS
#
pdf2hdi <- function(pdf,xMIN=0,xMAX=1,cred=0.95,Print=FALSE){
nxx=1001
xx=seq(xMIN,xMAX,length=nxx)
xx=xx[-nxx]
xx=xx[-1]
xmin=xx[1]
xmax=xx[nxx-2]
AUC=integrate(pdf,xmin,xmax)$value
cdf <- function(x) integrate(pdf,xmin,x)$value/AUC
省12
263: 2017/11/20(月)20:38 ID:xJug4kDO(2/4) AAS
pdf2hdi <- function(pdf,xMIN=0,xMAX=1,cred=0.95,Print=TRUE){
nxx=1001
xx=seq(xMIN,xMAX,length=nxx)
xx=xx[-nxx]
xx=xx[-1]
xmin=xx[1]
xmax=xx[nxx-2]
AUC=integrate(pdf,xmin,xmax)$value
PDF=function(x)pdf(x)/AUC
cdf <- function(x) integrate(PDF,xmin,x)$value
省16
264: 2017/11/20(月)21:23 ID:xJug4kDO(3/4) AAS
pdf2hdi <- function(pdf,xMIN=0,xMAX=1,cred=0.95,Print=TRUE){
nxx=1001
xx=seq(xMIN,xMAX,length=nxx)
xx=xx[-nxx]
xx=xx[-1]
xmin=xx[1]
xmax=xx[nxx-2]
AUC=integrate(pdf,xmin,xmax)$value
PDF=function(x)pdf(x)/AUC
cdf <- function(x) integrate(PDF,xmin,x)$value
省16
265: 2017/11/20(月)22:09 ID:xJug4kDO(4/4) AAS
According to this posting, 2chｽﾚ:doctor as many as 15 freshmen flunk.
Let's assume there are 120 students in one class and the number of flunker is distributed as poisson distribution.
In what range will students are expected to flunk next year? Calculate the number with 99% confidence interval.
pdf2hdi <- function(pdf,xMIN=0,xMAX=1,cred=0.95,Print=TRUE,FUN=FALSE){
nxx=1001
xx=seq(xMIN,xMAX,length=nxx)
xx=xx[-nxx]
xx=xx[-1]
xmin=xx[1]
xmax=xx[nxx-2]
省20
266: 2017/11/22(水)10:27 ID:WHgLvp/3(1/3) AAS
ab2mv<-function(a,b){
m<-a/(a+b)
v<-m*(1-m)/(a+b+1)
mv<-c(m,v)
return(mv)
}
mv2ab<-function(m,v){
a=(-m^3+m^2-m*v)/v
b=(m^3-2*m^2+m*v+m-v)/v
ab<-c(a,b)
省9
267: 2017/11/22(水)12:17 ID:WHgLvp/3(2/3) AAS
Scale <- function(x) (x-mean(x))/sqrt(sum((x-mean(x))^2)/(length(x)-1))
m=50;s=10
rn=rnorm(1000,m,s)
mean(rn)
sd(rn)
y=Scale(rn)*s+m
mean(y)
sd(y)
268: 2017/11/22(水)19:57 ID:WHgLvp/3(3/3) AAS
動画ﾘﾝｸ[YouTube]
269: 2017/11/23(木)15:22 ID:Ue5tZuwc(1/2) AAS
par(mfrow=c(2,2))
theta=0.5
NN=1000
N=1:NN
flip=numeric()
for(i in N){
flip=append(flip,rbinom(1,1,theta))
}
z=cumsum(flip)
z_N=z/N
省14
270: 2017/11/23(木)15:22 ID:Ue5tZuwc(2/2) AAS
hdi=NULL
for(i in N){
y=flip[1:i]
s=rep(1,i)
data=data.frame(y,s)
Ntotal=i
Nsubj=1
dataList=list(y=y,s=s,Ntotal=Ntotal,Nsubj=Nsubj)
js=genMCMC2(data)
hdi=rbind(hdi,HDIofMCMC(as.matrix(js)))
省7
271: 2017/11/25(土)12:10 ID:3kOACkIe(1) AAS
Metropolis Algo
source('DBDA2E-utilities.R')
Bern <- function(x) rbinom(1,1,x)
Metro <- function(
SD=0.02,
.z=14,
.N=20,
a=1,
b=1,
k=50000,
省18
272: 2017/11/26(日)19:08 ID:11LwtYgG(1/3) AAS
# JAGS for proportion
graphics.off()
rm(list=ls())
zi=100; Ni=100; zj=10; Nj=10
(y=c(rep(1,zi),rep(0,Ni-zi),rep(1,zj),rep(0,Nj-zj)))
(s=as.numeric(factor(c(rep('D',Ni),rep('U',Nj)))))
myData=data.frame(y=y,s=s)
Ntotal = length(y)
Nsubj = length(unique(s))
dataList = list(
省18
273: 2017/11/26(日)19:09 ID:11LwtYgG(2/3) AAS
writeLines( modelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList)
#jags.model(file, data, inits,n.chains = 1, n.adapt=1000, quiet=FALSE)
update( jagsModel)
codaSamples = coda.samples( jagsModel , variable="theta", n.iter=10000 )
#coda.samples(model, variable.names, n.iter, thin = 1, na.rm=TRUE, ...)
summary(codaSamples)
plot(codaSamples,col=sample(colors()))
mcmcMat=as.matrix(codaSamples)
par(mfrow=c(2,2))
省7
274: 2017/11/26(日)19:09 ID:11LwtYgG(3/3) AAS
source("Kruschke_tools.R") # for genMCMC2
(myData=data.frame(y,s))
mcmcCoda = genMCMC2( data=myData , numSavedSteps=10000,a=1,b=1)
# mcmcCoda = genMCMC( data=myData , numSavedSteps=10000)
# genMCMC
# Display diagnostics of chain, for specified parameter:
source('Jags-Ydich-XnomSsubj-MbernBeta.R')
diagMCMC( mcmcCoda , parName="theta[1]" )
diagMCMC( mcmcCoda , parName="theta[2]" )
# Display numerical summary statistics of chain:
省14
275: 2017/11/28(火)18:49 ID:QFCizNPN(1/3) AAS
logistic <- function(x,gain,threshold){
1/(1 + exp(-gain*(x-threshold)))
}
b2gt <- function(b0,b1,b2){
gain=sqrt(b1^2+b2^2)
threshold=-b0/gain
return(gain=gain,threshold=threshold)
}
上下前次1-新書関写板覧索設栞歴
あと 472 ﾚｽあります
ｽﾚ情報 赤ﾚｽ抽出 画像ﾚｽ抽出 歴の未読ｽﾚ AAｻﾑﾈｲﾙ
ぬこの手 ぬこTOP 0.427s*