美文网首页数据科学与R语言R语言小技能
R语言——自定义函数求置信区间

R语言——自定义函数求置信区间

作者: 弓二竹 | 来源:发表于2017-12-10 19:04 被阅读4次

    【代码】:

    #求单正态均值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)
    }
    
    得到的函数

    选修课作业,自己写函数求单正态样本均值、方差置信区间,两个正态样本均值差、方差比的置信区间。
    求解时正态方差和均值默认为未知,函数具体样子可以参考题图。
    本来是想输出一段话,但是我不知道怎么换行,所以将就着看吧。
    我觉得我应该没有求错,如果求错了,麻烦告诉一声,谢啦!!

    相关文章

      网友评论

        本文标题:R语言——自定义函数求置信区间

        本文链接:https://www.haomeiwen.com/subject/drvvixtx.html