######################################################### ## ## ## noncentral t-distribution ## ## ## ######################################################### # 本程序用于求非中心t分布非中心参数在自由度已知情况下的最大似然估计及相应的对数似然函数值 # 注意对数似然函数值为包含了非中心参数那一部分的对数似然函数的最大值 ############################################## ## ## ## > source("noncentT.txt") ## ## ## ## > noncentT("Exam.txt") ## ## ## ############################################## # 输入变量data表示数据文件 # 第一列为自由度,其它各列为样本数据 # 参考文献: # 汪鹏,周基元.非中心t分布非中心参数的最大似然估计.2012.(已投稿). noncentT <- function(data) { # 读取数据文件 data <- as.vector(read.table(data, as.is = T)) data0 <- as.numeric(data) n.sample <- length(data0) - 1 n <- data0[1] # degree of freedom x <- data0[-1] # 样本 ############################################## ## ## ## 函数g, 可以在下面修改 ## ## ## ############################################## # g是一个关于非中心参数的单调函数,其中自由度是参变量,非中心参数δ与待估参数k之间满足δ=g(n,k) g <- function(n, k) { ## 默认g(n,k) = k,即直接估计非中心参数δ,然而可以在此处做修改 k } #################################### ## ## ## iv函数用于计算满足精度的初始值 ## ## ## #################################### iv <- function(u){ value1 <- 0; value2 <- sqrt(pi/2) for (i in 1:13) { value1<-value1+(-1)^(i-1)*u^(2*i-2)/(2^(i-1)*factorial(i-1)) value2<-value2+(-1)^(i-1)*u^(2*i-1)/(2^(i-1)*(2*i-1)*factorial(i-1)) } value1<-value1+(-1)^14*u^(26)/(2^14*factorial(14)) list(value1=value1,value2=value2) } ############### number1函数说明 ################### ## ## ## 用途:计算运用迭代公式(1)误差增长倍数的上限 ## ## ## ################################################## number1 <- function(n,u){ if (u<-1||u>0) r<-1 if (n==1) r<-0 else{ a <- as.vector(n) a[1] <- u+iv(u)$value1/iv(u)$value2 b <- 1 for (i in 1:(n-1)){ a[i+1] <- u+i/a[i] if (i%%2==0) b<-b*abs(i*(i-1)/((u*a[i-1]+i-1)*(u^2+i-abs(u)*(sqrt(i)-u)))) } r <-2.7* b*1e-14 } return(r) } ################## number2函数说明 ############### ## ## ## 用途:计算运用迭代公式(2)所需的迭代次数m ## ## ## ################################################## number2<-function(n,u){ i <- n+2 repeat { c <- 1 a <- as.vector(i-n+1) l1 <- max(i/(sqrt(i+1)+abs(u)),sqrt(i-1)-abs(u)/2) l2 <- min(-i/u,sqrt(i)) a[i] <- (l1+l2)/2 for (j in i:(n+1)){ a[j-1] <- (j-1)/(a[j]-u) if ((j-n)%%2==0) { c <- c*(j-1)/(j-1+u^2-u*a[j]) } } d <- c*0.5*(l2-l1)/l1 if (d<=1e-10) break else { if (n<=100) i <- i+2*floor(n/2) else i <- i+max(2,2*floor(n/12)) } } return(i) } ######### f函数说明 ############### ## ## ## 用途:计算迭代公式 ## ## 参数说明: ## ## ## ## k:待估参数 ## ## x:样本值 ## ## n:自由度 ## ################################### f <- function(n,k,x){ u <- g(n,k)*x/sqrt(x^2+n) value<-log(sqrt(2 * pi) )+pnorm(u,log.p=TRUE) if (u>=0||(number1(n,u)<=1e-10)) { a <- as.vector(n) a[1] <- u + exp(-u^2/2)/(sqrt(2*pi)*pnorm(u)) if(n == 1) value<- log(a[1]) + value else { b <- log(a[1]) for (i in 1:(n - 1)) { a[i + 1] <- u + i/a[i] b <- b + log(a[i + 1]) } value <- b + value } } else { r <- number2(n,u) a <- as.vector(r) l1 <- max(0,sqrt(r-1)-abs(u)/2) l2 <- min(-r/u,sqrt(r)) a[r] <- (l1+l2)/2 c <- 0 for (i in r:2) { a[i-1] <- (i-1)/(a[i]-u) if (i<=n+1) c<-c+log(a[i-1]) } value <- c + value } return(value) } ############## 计算最大似然估计及对数似然函数最大值 ################## ## 参数说明: ## ## ## ## freedom:似然函数 ## ## estimator:最大似然估计 ## ## lmaxvalue:对数似然函数最大值 ## ## ## #################################################################### freedom<-function(k){ logL <- 0 for (i in (1:n.sample)){ logL <- logL+n*g(n,k)^2/(2*(x[i]^2+n))-f(n,k,x[i]) } return(logL) } upper<-mean(x)+1.96*sd(x) lower<-mean(x)-1.96*sd(x) estimator <- optimize(freedom,c(lower,upper))$minimum maxvalue <- -optimize(freedom,c(lower,upper))$objective list(df = n, MLE = estimator, max_log_likelihood = maxvalue) }