CDOoDocuments.StdDocumentDescDocuments.DocumentDescContainers.ViewDescViews.ViewDescStores.StoreDescqDocuments.ModelDescContainers.ModelDescModels.ModelDescStores.ElemDesc TextViews.StdViewDescTextViews.ViewDescaTextModels.StdModelDescTextModels.ModelDesckcTextModels.AttributesDesc1$Courier New  1$Courier New < CourierPS1 CourierPSi   #------------------------------------------------------------------------------------------------------------------------------ Version 1 #------------------------------------------------------------------------------------------------------------------------------ model{ # model's likelihood for (i in 1:n){ y[i]~dnorm( mu, tau ) y.rep[i]~dnorm( mu, tau ) } # # log-normal prior mu ~ dnorm( 0.0, 0.001) tau ~ dgamma( 0.001, 0.001 ) s<-1/sqrt(tau) # # predictive cummulative frequencies for each y_i for (i in 1:n){ for (k in 1:n){ pred.lower.yi[i,k] <- step( y[i] - y.rep[k] ) } F.pred[i] <- sum( pred.lower.yi[i,1:n] )/n } #-------------------------------------------------------------------------- # residuals for (i in 1:n){ # ordered replicated/predicted data ordered.y.rep[i]<- ranked( y.rep[1:n], i ) ordered.y[i] <- ranked( y[1:n], i ) rs[i] <- (y[i]-mu)*sqrt(tau) pr[i] <- step(y.rep[i]-y[i]) pro[i]<- step(ordered.y.rep[i]-ordered.y[i]) } #-------------------------------------------------------------------------- # gof's #-------------------------------------------------------------------------- e<-0.0000001 # precision parameter # zcut = cut points for N(0,1) # calculation of rescaled cut points # ycut = y' (cut points for y) for (i in 1:ncut){ ycut[i] <- mu + zcut[i]*s } # # calculation of observed cumulative frequencies (for actual and replicated data) for (i in 1:ncut){ for (j in 1:n){ bin.freq.obs[i,j] <- step( ycut[i]-ordered.y[j] ) bin.freq.rep[i,j] <- step( ycut[i]-ordered.y.rep[j] ) } F.obs[i] <- sum( bin.freq.obs[i,1:n] )/n F.rep[i] <- sum( bin.freq.rep[i,1:n] )/n F.exp[i] <- phi( zcut[i] ) } # # calculation of frequencies values within interval y_k', y_{k-1}' f.exp[1] <- F.exp[1] f.obs[1] <- F.obs[1] f.rep[1] <- F.rep[1] for (i in 2:ncut ){ f.obs[i] <- F.obs[i]-F.obs[i-1] f.rep[i] <- F.rep[i]-F.rep[i-1] f.exp[i] <- F.exp[i]-F.exp[i-1] } f.obs[ncut+1] <- 1 - F.obs[ncut] f.rep[ncut+1] <- 1 - F.rep[ncut] f.exp[ncut+1] <- 1 - F.exp[ncut] # for (i in 1:(ncut+1)){ # setting zero expected frequencies equal to e f.exp2[i] <- f.exp[i] + e*equals(f.exp[i],0) chisq.f.obs.vec[i] <- pow(f.obs[i]-f.exp2[i],2)/(f.exp2[i]*(1-f.exp2[i])) chisq.f.rep.vec[i] <- pow(f.rep[i]-f.exp2[i],2)/(f.exp2[i]*(1-f.exp2[i])) } # chisq values chisq.f.obs <- sum( chisq.f.obs.vec[] ) chisq.f.rep <- sum( chisq.f.rep.vec[] ) # chisq p-value chisq.f.p <- step( chisq.f.rep - chisq.f.obs ) # 1st version of chi^2_F for (i in 1:(ncut)){ F.exp2[i] <- F.exp[i] + e*equals(F.exp[i],0) chisq.F.obs.vec[i] <- pow(F.obs[i]-F.exp2[i],2)/(F.exp2[i]*(1-F.exp2[i])) chisq.F.rep.vec[i] <- pow(F.rep[i]-F.exp2[i],2)/(F.exp2[i]*(1-F.exp2[i])) } # chisq values chisq.F.obs <- sum( chisq.F.obs.vec[] ) chisq.F.rep <- sum( chisq.F.rep.vec[] ) # chisq p-value chisq.F.p <- step( chisq.F.rep - chisq.F.obs ) # 2nd version of chi^2_F for (i in 1:n){ F.obs2[i] <- rank( ordered.y[], i )/n F.rep2[i] <- rank( ordered.y.rep[], i )/n F.exp.obs2[i] <- phi( (ordered.y[i]-mu)*sqrt(tau) ) F.exp.rep2[i] <- phi( (ordered.y.rep[i]-mu)*sqrt(tau) ) chisq.F2.obs.vec[i] <- pow(F.obs2[i]-F.exp.obs2[i],2)/(F.exp.obs2[i]*(1-F.exp.obs2[i])) chisq.F2.rep.vec[i] <- pow(F.rep2[i]-F.exp.rep2[i],2)/(F.exp.rep2[i]*(1-F.exp.rep2[i])) } # chisq values chisq.F2.obs <- sum( chisq.F2.obs.vec[] ) chisq.F2.rep <- sum( chisq.F2.rep.vec[] ) # chisq p-value chisq.F2.p <- step( chisq.F2.rep - chisq.F2.obs ) # KS statistic # for (i in 1:n){ F.diff.obs[i] <- abs( F.obs2[i] - F.exp.obs2[i] ) F.diff.rep[i] <- abs( F.rep2[i] - F.exp.rep2[i] ) } ks.obs <- ranked( F.diff.obs[], n ) ks.rep <- ranked( F.diff.rep[], n ) # chisq p-value ks.p <- step( ks.rep - ks.obs ) } INITS list( mu=0.0, tau=1 ) CUTPOINTS list( ncut=19, zcut=c(-1.64, -1.28, -1.04, -0.84, -0.67, -0.52, -0.39, -0.25, -0.13, 0, 0.13, 0.25, 0.39, 0.52, 0.67, 0.84, 1.04, 1.28, 1.64) ) list( ncut=9, zcut=c(-1.28, -0.84, -0.52, -0.25, 0, 0.25, 0.52, 0.84, 1.28) ) list( ncut=39, zcut=c(-1.96, -1.64, -1.44, -1.28, -1.15, -1.04, -0.93, -0.84, -0.76, -0.67, -0.6, -0.52, -0.45, -0.39, -0.32, -0.25, -0.19, -0.13, -0.06, 0, 0.06, 0.13, 0.19, 0.25, 0.32, 0.39, 0.45, 0.52, 0.6, 0.67, 0.76, 0.84, 0.93, 1.04, 1.15, 1.28, 1.44, 1.64, 1.96 ) ) DATA1 (LIST) list( n=19, y=c(0.51, 0.1, -2.53, 1, 0.65, -0.95, 2.76, 1.33, 0.25, 1.48, -0.25, -0.45, 2.11, -0.76, -1.51, -0.35, 0.18, 1.35, -0.18) ) DATA2 (LIST) list( n=20, y=c(0.51, 0.1, -2.53, 1, 0.65, -0.95, 2.76, 1.33, 0.25, 1.48, -0.25, -0.45, 2.11, -0.76, -1.51, -0.35, 0.18, 1.35, -0.18, 5) ) DATA3 (LIST) list( n=20, y=c(0.51, 0.1, -2.53, 1, 0.65, -0.95, 2.76, 1.33, 0.25, 1.48, -0.25, -0.45, 2.11, -0.76, -1.51, -0.35, 0.18, 1.35, -0.18, 10) ) DATA4 (from log(gamma(2,1)) list(n=19, y=c(1.43, 1.16, 1.64, 1.06, -1.76, -0.29, -1.37, -0.72, -0.94, 1.2, 1.41, 0.74, 0.77, 0.71, 1.83, 1.58, 0.18, 0.18, 1.43)) DATA5 (from gamma(2,1) ) list(n=100, y=c(2.02, 0, 0.34, 0.19, 0.03, 0.02, 0.01, 0.21, 0.04, 0, 0, 0.79, 0, 0.04, 0.01, 0.04, 0, 0.22, 0.03, 0.04, 0, 0.11, 0, 0.1, 0.31, 0.05, 0.13, 0, 0.67, 0.06, 0.08, 0.26, 0, 0, 0, 0, 0.05, 0.98, 0.06, 0, 0, 0, 0, 1.07, 0.08, 0, 0, 2.11, 0, 0.01, 0, 0, 0.01, 0, 0.01, 0.48, 0.07, 0, 0, 0.22, 0.06, 0, 0, 0, 0, 0.04, 0.49, 1.36, 0.01, 0, 0, 0, 0, 0, 0, 0, 0, 0.06, 0.01, 0, 0.66, 0.24, 0.02, 0.02, 0.05, 0, 0.15, 0.07, 2.34, 0, 0.02, 0.12, 0.38, 0, 0.02, 0.01, 0, 0.79, 0.03, 0.21)) TextControllers.StdCtrlDescTextControllers.ControllerDescContainers.ControllerDescControllers.ControllerDesc TextRulers.StdRulerDescTextRulers.RulerDescTextRulers.StdStyleDescTextRulers.StyleDescZTextRulers.AttributesDesc$ ZGo * ,[ @Documents.ControllerDesc t]s ' `h*