R語言——自定義函式求置信區間
阿新 • • 發佈:2019-01-09
選修課作業,自己寫函式求單正態樣本均值、方差置信區間,兩個正態樣本均值差、方差比的置信區間。求解時正態方差和均值預設為未知,函式具體樣子可以參考題圖。本來是想輸出一段話,但是我不知道怎麼換行,所以將就著看吧。我覺得我應該沒有求錯,如果求錯了,麻煩告訴一聲,謝啦!!#求單正態均值mu的置信區間 #引數依次為置信水平alpha,正態樣本x,已知總體方差(預設為未知) mu <- function(alpha,x,sigma=NA){ n <- length(x) meanx <- mean(x) if(is.na(sigma)){ t1 <- qt(1-alpha/2,n-1) t2 <- qt(1-alpha,n-1) mu11 <- meanx - t1*sqrt(sum((x-meanx)^2)/(n-1))/sqrt(n) mu12 <- meanx + t1*sqrt(sum((x-meanx)^2)/(n-1))/sqrt(n) mu21 <- meanx + t2*sqrt(sum((x-meanx)^2)/(n-1))/sqrt(n) mu22 <- meanx - t2*sqrt(sum((x-meanx)^2)/(n-1))/sqrt(n) } else{ u1 <- qnorm(1-alpha/2,0,1) u2 <- qnorm(1-alpha,0,1) mu11 <- meanx - u1*sigma/sqrt(n) mu12 <- meanx + u1*sigma/sqrt(n) mu21 <- meanx + u2*sigma/sqrt(n) mu22 <- meanx - u2*sigma/sqrt(n) } string1 <- paste('以1-',alpha,'為置信水平的mu雙側置信區間為:[',mu11,', ',mu12,']。',sep='') string2 <- paste('以1-',alpha,'為置信水平的mu單側置信區間上限為:',mu21,'。',sep='') string3 <- paste('以1-',alpha,'為置信水平的mu單側置信區間下限為:',mu22,'。',sep='') string <- data.frame(Confidence_Interval=c(string1,string2,string3)) return(string) } #求單正態方差sigma的置信區間 #引數依次為置信水平alpha,正態樣本x,已知總體均值(預設為未知) sigma <- function(alpha,x,mu=NA){ n <- length(x) if(is.na(mu)){ meanx <- mean(x) chisq11 <- qchisq(1-alpha/2,n-1) chisq12 <- qchisq(alpha/2,n-1) chisq21 <- qchisq(alpha,n-1) chisq22 <- qchisq(1-alpha,n-1) sigma11 <- sqrt(sum((x-meanx)^2)/chisq11) sigma12 <- sqrt(sum((x-meanx)^2)/chisq12) sigma21 <- sqrt(sum((x-meanx)^2)/chisq21) sigma22 <- sqrt(sum((x-meanx)^2)/chisq22) } else{ chisq11 <- qchisq(1-alpha/2,n) chisq12 <- qchisq(alpha/2,n) chisq21 <- qchisq(alpha,n) chisq22 <- qchisq(1-alpha,n) sigma11 <- sqrt(sum((x-mu)^2)/chisq11) sigma12 <- sqrt(sum((x-mu)^2)/chisq12) sigma21 <- sqrt(sum((x-mu)^2)/chisq21) sigma22 <- sqrt(sum((x-mu)^2)/chisq22) } string1 <- paste('以1-',alpha,'為置信水平的sigma雙側置信區間為:[',sigma11,', ',sigma12,']。',sep='') string2 <- paste('以1-',alpha,'為置信水平的sigma單側置信區間上限為:',sigma21,'。',sep='') string3 <- paste('以1-',alpha,'為置信水平的sigma單側置信區間下限為:',sigma22,'。',sep='') string <- data.frame(Confidence_Interval=c(string1,string2,string3)) return(string) } #求兩個正態均值差(mux-muy)的置信區間 #引數依次為置信水平alpha,正態樣本x,正態樣本y, #已知x總體方差sigmax(預設為未知),已知y總體方差sigmay(預設為未知) mux_muy <- function(alpha,x,y,sigmax=NA,sigmay=NA){ if(is.na(sigmax)|is.na(sigmay)){ meanx <- mean(x) meany <- mean(y) m <- length(x) n <- length(y) sx <- sqrt(sum((x-meanx)^2)/(m-1)) sy <- sqrt(sum((y-meany)^2)/(n-1)) sw <- sqrt((m-1)*sx^2/(m+n-2)+(n-1)*sy^2/(m+n-2)) mu11 <- (meanx-meany)+qt(1-alpha/2,m+n-2)*sw*sqrt(1/m+1/n) mu11 <- (meanx-meany)-qt(1-alpha/2,m+n-2)*sw*sqrt(1/m+1/n) } else{ meanx <- mean(x) meany <- mean(y) m <- length(x) n <- length(y) sx <- sqrt(sum((x-mux)^2)/m) sy <- sqrt(sum((y-muy)^2)/n) mu11 <- (meanx-meany)+qt(1-alpha/2,m+n)*sw*sqrt(1/m+1/n) mu11 <- (meanx-meany)-qt(1-alpha/2,m+n)*sw*sqrt(1/m+1/n) } string1 <- paste('以1-',alpha,'為置信水平的mux-muy雙側置信區間為:[',mu11,', ',mu12,']。',sep='') return(string1) } #求兩個正態標準差比sigmax/sigmay的置信區間 #引數依次為置信水平alpha,正態樣本x,正態樣本y, #已知x總體均值mux(預設為未知),已知y總體均值muy(預設為未知) sigmax_sigmay <- function(alpha,x,y,mux=NA,muy=NA){ alpha <- alpha mux <- mux muy <- muy if(is.na(mux)|is.na(muy)){ meanx <- mean(x) m <- length(x) meany <- mean(y) n <- length(y) F1 <- qf(1-alpha/2,m-1,n-1) F2 <- qf(alpha/2,m-1,n-1) sigma11 <- 1/F1*sum((x-meanx)^2)*(n-1)/sum((y-meany)^2)/(m-1) sigma12 <- 1/F2*sum((x-meanx)^2)*(n-1)/sum((y-meany)^2)/(m-1) } else{ m <- length(x) n <- length(y) F1 <- qf(1-alpha/2,m,n) F2 <- qf(alpha/2,m,n) sigma11 <- 1/F1*sum((x-mux)^2)*n/sum((y-muy)^2)/m sigma12 <- 1/F2*sum((x-mux)^2)*n/sum((y-muy)^2)/m } string1 <- paste('以1-',alpha,'為置信水平的sigmax-sigmay雙側置信區間為:[',sigma11,', ',sigma12,']。',sep='') return(string1) }