# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # Skutecne pokryti ruznych intervalu spolehlivosti pokryti <- function(n=20, alpha=0.05){ # n ... rozsah vyberu # alpha ... hladina testu n.p0 <- 100; p0 <- (1:n.p0)/(n.p0+1) vysl <- list(); vysl$wald <- numeric(n.p0) vysl$wilson <- numeric(n.p0) vysl$logit <- numeric(n.p0); vysl$presny <- numeric(n.p0); delka <- list(); delka$wald <- numeric(n.p0) delka$wilson <- numeric(n.p0) delka$logit <- numeric(n.p0); delka$presny <- numeric(n.p0); # Mozne odhady Xn <- 0:n; phat <- Xn/n; kvant <- qnorm(1-alpha/2); IS.wald <- cbind(phat - qnorm(1-alpha/2)*sqrt(phat*(1-phat)/n), phat + qnorm(1-alpha/2)*sqrt(phat*(1-phat)/n)); IS.wilson <- 1/(1 + kvant^2/n)*cbind(phat + kvant^2/(2*n) - kvant * sqrt(phat*(1 - phat)/n+kvant^2/(4*n^2)), phat + kvant^2/(2*n) + kvant * sqrt(phat*(1 - phat)/n+kvant^2/(4*n^2))); # Logit thetan <- log(phat/(1-phat)); Dn <- sqrt(n/(Xn * (n - Xn))) IS.log.sance <- cbind(thetan - Dn*qnorm(1-alpha/2), thetan + Dn*qnorm(1-alpha/2)); IS.logit <- exp(IS.log.sance)/(1+exp(IS.log.sance)) IS.logit[1,] <- c(0,0); IS.logit[n+1,] <- c(1,1); # Presny IS qL <- qf(alpha/2, df1=2*Xn, df2=2*(n-Xn+1)); qL[is.na(qL)] <- 0; qU <- qf(1-alpha/2, df1=2*(Xn+1), df2=2*(n-Xn)); qU[is.na(qU)] <- 1; IS.presny <- cbind(Xn*qL/(Xn*qL + n - Xn + 1), (Xn+1)*qU/((Xn+1)*qU + n - Xn)) # # Je pokryto? fce <- function(x, p) as.logical(x[,1] <= p)*as.logical(x[,2] >= p); fce_delka <- function(x) x[,2]-x[,1] for(i in 1:n.p0){ psti <- dbinom(0:n, size=n, prob=p0[i]); vysl$wald[i] <- sum(fce(IS.wald, p0[i])*psti) vysl$wilson[i] <- sum(fce(IS.wilson, p0[i])*psti) vysl$logit[i] <- sum(fce(IS.logit, p0[i])*psti) vysl$presny[i] <- sum(fce(IS.presny, p0[i])*psti) delka$wald[i] <- sum(fce_delka(IS.wald)*psti) delka$wilson[i] <- sum(fce_delka(IS.wilson)*psti) delka$logit[i] <- sum(fce_delka(IS.logit)*psti) delka$presny[i] <- sum(fce_delka(IS.presny)*psti) } # print(vysl); plot(p0, vysl$wald, main = paste("Skutecne pokryti pro n =", n), type="l", ylim=range(vysl), xlab="Hodnota parametru p", ylab="Pravdepodobnost pokryti", lwd=2); lines(p0, vysl$wilson, col="red", lwd=2); lines(p0, vysl$logit, col="green", lwd=2); lines(p0, vysl$presny, col="blue", lwd=2); abline(h=0.95, lty="dotted"); legend(0.5, min(vysl$wald, vysl$logit), lty=rep("solid", 4), lwd=2, col=c("black", "red", "green", "blue"), legend=c("Wald", "Wilson", "Logit", "Presny"), xjust=0.5, yjust=0) delka <- round(matrix(unlist(delka), nrow=n.p0, ncol=4),2); colnames(delka) <- c("Wald", "Wilson", "Logit", "Presny"); rownames(delka) <- round(p0,2) print("Prumerna delka intervalu spolehlivosti:") print(delka[seq(1,n.p0, length=10),]) # invisible(vysl); }