R语言—自定义函数求置信区间的操作
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
R语⾔—⾃定义函数求置信区间的操作看代码吧~
#求单正态均值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语⾔【估计单侧置信区间】
以上为个⼈经验,希望能给⼤家⼀个参考,也希望⼤家多多⽀持。
如有错误或未考虑完全的地⽅,望不吝赐教。