[過去ログ] 臨床統計もおもしろいですよ、その1 [無断転載禁止]©2ch.net (747レス)
前次1-
抽出解除 必死チェッカー(本家) (べ) レス栞 あぼーん

このスレッドは過去ログ倉庫に格納されています。
次スレ検索 歴削→次スレ 栞削→次スレ 過去ログメニュー
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
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
前次1-
スレ情報 赤レス抽出 画像レス抽出 歴の未読スレ AAサムネイル

ぬこの手 ぬこTOP 0.765s*