R語言—自定義函數(shù)求置信區(qū)間的操作
更新時間:2021年04月20日 08:34:43 作者:弓二竹
這篇文章主要介紹了R語言—自定義函數(shù)求置信區(qū)間的操作,具有很好的參考價值,希望對大家有所幫助。一起跟隨小編過來看看吧
看代碼吧~
#求單正態(tài)均值mu的置信區(qū)間
#參數(shù)依次為置信水平alpha,正態(tài)樣本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雙側置信區(qū)間為:[',mu11,', ',mu12,']。',sep='')
string2 <- paste('以1-',alpha,'為置信水平的mu單側置信區(qū)間上限為:',mu21,'。',sep='')
string3 <- paste('以1-',alpha,'為置信水平的mu單側置信區(qū)間下限為:',mu22,'。',sep='')
string <- data.frame(Confidence_Interval=c(string1,string2,string3))
return(string)
}
#求單正態(tài)方差sigma的置信區(qū)間
#參數(shù)依次為置信水平alpha,正態(tài)樣本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雙側置信區(qū)間為:[',sigma11,', ',sigma12,']。',sep='')
string2 <- paste('以1-',alpha,'為置信水平的sigma單側置信區(qū)間上限為:',sigma21,'。',sep='')
string3 <- paste('以1-',alpha,'為置信水平的sigma單側置信區(qū)間下限為:',sigma22,'。',sep='')
string <- data.frame(Confidence_Interval=c(string1,string2,string3))
return(string)
}
#求兩個正態(tài)均值差(mux-muy)的置信區(qū)間
#參數(shù)依次為置信水平alpha,正態(tài)樣本x,正態(tài)樣本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雙側置信區(qū)間為:[',mu11,', ',mu12,']。',sep='')
return(string1)
}
#求兩個正態(tài)標準差比sigmax/sigmay的置信區(qū)間
#參數(shù)依次為置信水平alpha,正態(tài)樣本x,正態(tài)樣本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雙側置信區(qū)間為:[',sigma11,', ',sigma12,']。',sep='')
return(string1)
}
選修課作業(yè),自己寫函數(shù)求單正態(tài)樣本均值、方差置信區(qū)間,兩個正態(tài)樣本均值差、方差比的置信區(qū)間。
求解時正態(tài)方差和均值默認為未知,函數(shù)具體樣子可以參考題圖。
本來是想輸出一段話,但是我不知道怎么換行,所以將就著看吧。
補充:R語言【估計單側置信區(qū)間】

以上為個人經驗,希望能給大家一個參考,也希望大家多多支持腳本之家。如有錯誤或未考慮完全的地方,望不吝賜教。
相關文章
R語言中dnorm,pnorm,qnorm和rnorm的區(qū)別淺析
正在學習R語言統(tǒng)計學的小伙伴,可能會被各種專有名詞所困擾,下面這篇文章主要給大家介紹了關于R語言中dnorm,pnorm,qnorm和rnorm區(qū)別的相關資料,需要的朋友可以參考下2022-12-12

