首页 > 编程语言 > R语言—自定义函数求置信区间的操作
2021
05-19

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

看代码吧~

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
#求单正态均值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语言【估计单侧置信区间】

在这里插入图片描述

以上为个人经验，希望能给大家一个参考，也希望大家多多支持自学编程网。如有错误或未考虑完全的地方，望不吝赐教。

编程技巧