四种方法结果的比较

最新推荐文章于 2023-12-23 04:48:33 发布
最新推荐文章于 2023-12-23 04:48:33 发布
阅读量759 收藏 1
点赞数 1
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/Cet777/article/details/104086190
版权

以下是关于整个模拟程序的解释与程序：
加载所需要的包，生成所需要的数据：目标分布服从威布尔分布，截断分布服从指数分布。

install.packages("rootSolve")
install.packages("Rsolnp")
library(Rsolnp)
library(rootSolve)
#################################################################
#####  define a new log function  log_n ###
log_n <- function(z,n){
  an=n^2
  ind=(z<1/an)
  y1 =   log(1/an)-1.5+2*an*z-(an*z)^2/2
  y1[!ind]=0
  z[ind] = 1
  y2=log(z)
  y=y1+y2 
  return(y)
}
################################################################################
####   Function used to calculate    log( gamma(N+1)/gamma(N-n+1) )  
digam<-function(n_big, n_small)
{ out=0
if(n_small>0 &  n_big>=n_small)  out=sum(log(n_big+1-c(1:n_small)))    
out
}

###########################################################
################  Function used to generate data
gen_obs <- function(wshape,wscale,erate,n_big){
  obs.x <- rweibull(n_big,shape = wshape,scale = wscale)
  obs.y <- rexp(n_big,rate = erate)   
  ind =(obs.x > obs.y)    
  obs <- data.frame(x=obs.x[ind],y=obs.y[ind])
  return(obs)
}

1，直接极大似然函数

max.like函数用于求解全似然函数的最优解：

#########################################################
######      full log-likelihood function
###
###   Function max.like is to calculate the estiamtor of [N,alpha,theta],N is n_big
###           para= [N,alpha,theta] in max.like , corresponding estimate is par_hat
###
###   Before the function log_likelihood , define a data frame-obs
###
##################################################################
max.like<-function(obs){
  n = nrow(obs)
  x = obs[, 1]
  y = obs[, 2]
  
  log_likelihood <- function(para){ #negative likelihood
    fn <- function(lambda){
      ter_log <- log_n(dexp(y, para[3]), n)  -  log_n(1+lambda*(pexp(x,para[3])-para[2]),n)
      t <- sum(ter_log)+(para[1]-n)*log(1-para[2])+digam(para[1],n)
      return(t)
    }
    z <- (-1)*optimize(fn, interval=c(-1000,1000), tol=0.0001)$objective
    return(z)
  }
  par_hat <- nlminb( c(n,0.5, 1), log_likelihood, lower=c(n, 0.0001, 0.001), upper=c(100*n, 0.999, 100))$par
  return(par_hat)
}

max.like2函数用于当N取真值N0时，求 α \alpha α和 θ \theta θ的极大似然估计，进而计算N的似然比函数，进一步构建全似然方法的覆盖率。这里，在模拟的时候每次N设定不同的值时，这里的N0都要调整为相同的值，否则结果偏差很大。

###################################################################################################
########  calculate full likelihood confidence interval for parameter  par
###
###   max.like2 is to calculate the estiamtor of [alpha,theta] when given true n_big
###             para = [alpha,theta] in max.like2,  corresponding estimate is par_Nhat
###
###   con_int is to calculate (1-a) Interval lower bound , para_star is a parameter vector
###           para_star = [par_hat,par_Nhat,N0,1-a] 
###
###
####################################################################
max.like2<-function(obs){
  n = nrow(obs)
  x = obs[, 1]
  y = obs[, 2]
  
  log_likelihood.N0 <- function(para){ #negative likelihood
    # N take true value:N0
    N0=300
    fn <- function(lambda){
      ter_log <- log_n(dexp(y, para[2]), n)  -  log_n(1+lambda*(pexp(x,para[2])-para[1]),n)
      t <- sum(ter_log)+digam(N0,n)+(N0-n)*log(1-para[1])
      return(t)
    }
    z <- (-1)*optimize(fn, interval=c(-500,500), tol=0.0001)$objective
    return(z)
  }
  par_Nhat <- nlminb(c(0.5, 1),log_likelihood.N0,lower=c(0.0001, 0.001), upper=c(0.999, 100))$par
  return(par_Nhat)
}

计算全似然方法的覆盖率：

###########################################################################################
#####   Construct (1-a) confidence interval lower bound  for full likelihood 
###
con_int <- function(obs,para_star,N_true){
  n = nrow(obs)
  x = obs[, 1]
  y = obs[, 2]
  
  log_likelihood <- function(para){ ##negative likelihood
    fn <- function(lambda){
      ter_log <- log_n(dexp(y, para[3]), n)  -  log_n(1+lambda*(pexp(x,para[3])-para[2]),n)
      t <- sum(ter_log)+(para[1]-n)*log(1-para[2])+digam(para[1],n)
      return(t)
    }
    z <- optimize(fn, interval=c(-500,500), tol=0.001)$objective
    return(z)
  }
  
  clog_likelihood <- function(N){  # N is scalar
    fn <- function(lambda){
      ter_log <- log_n(dexp(y,para_star[5]),n) - log_n(1+lambda*(pexp(x,para_star[5])-para_star[4]),n)
      t <- sum(ter_log)+(N-n)*log(1-para_star[4])+digam(N,n)
      return(t)
    }
    z <- optimize(fn, interval=c(-500,500), tol=0.001)$objective
    return(z)
  }
  R2N <- 2*(log_likelihood(para = para_star[1:3]) - clog_likelihood(N=N_true))
  
  ifelse(R2N <= qchisq(tail(para_star,1),df=1),1,0)
}
#####################################################################################
#########   calculate Coverage percentage for full likeligood
##   INPUT:
##         int : Vector 0 and 1 to determine if the true value is in the confidence interval
##               1 represents the truth value in the confidence interval
##         K :  Number of simulations
##
##   OUTPUT:
##         P_f : Coverage percentage for full likrlihood
##
######################################################
p_f <- function(K,int){
  sum(int)/K
}

条件似然（Wang）的极大似然估计：

######################################################################
##########     Wang: condition likelihood  
###
###   Function max.clike is to calculate the estiamtor of [N,alpha,theta],N is n_big
###           para = [N,alpha,theta] in max.clike , corresponding estimate is par_tilde
###
###   sigma : the variance of the asymptotic distribution
###
###
###
###############################################################
max.clike<-function(obs){ 
  n = nrow(obs)
  x = obs[, 1]
  y = obs[, 2]
  
  llc <- function(theta){ 
    ter_y = log_n(dexp(y,theta),n)-log_n(pexp(x,theta),n)
    -sum(ter_y)
  }
  
  theta_tilde = nlminb(1, llc, lower=c(0.001), upper=c(100))$par  
  N_tilde = sum(1/pexp(x, theta_tilde))
  alpha_tilde = n/N_tilde
  par_tilde <- c(N_tilde, alpha_tilde, theta_tilde)
  #names(par_tilde) <- c("N_tilde", "alpha_tilde", "theta_tilde")
  return(par_tilde)
}

计算条件似然方法的覆盖率：

max.clike<-function(obs){ 
  n = nrow(obs)
  x = obs[, 1]
  y = obs[, 2]
  
  llc <- function(theta){ 
    ter_y = log_n(dexp(y,theta),n)-log_n(pexp(x,theta),n)
    -sum(ter_y)
  }
  
  theta_tilde = nlminb(1, llc, lower=c(0.001), upper=c(100))$par  
  N_tilde = sum(1/pexp(x, theta_tilde))
  alpha_tilde = n/N_tilde
  par_tilde <- c(N_tilde, alpha_tilde, theta_tilde)
  #names(par_tilde) <- c("N_tilde", "alpha_tilde", "theta_tilde")
  return(par_tilde)
}
#########################################################################
########################  Estimates of sigma #######
sigma_cest<- function(obs,theta_tilde,N_tilde){
  n = nrow(obs)
  x = obs[, 1]
  y = obs[, 2]
  #estimate of density function of truncated variable
  g_tilde <- function(t){
    ifelse(t >= 0,theta_tilde*exp(-theta_tilde*t),0)
  }
  #derivative of g about theta_tilde
  g1_tilde <- function(t){ 
    ifelse(t >= 0,(1+theta_tilde^2)*exp(-theta_tilde*t),0)
  }
  #estimate of distribution function 
  G_tilde <- function(t){
    ifelse(t>=0,1-exp(-theta_tilde*t),0)
  }
  #derivative of G about theta
  G1_tilde <- function(t){
    ifelse(t>=0,t*exp(-theta_tilde*t),0)  
  }
  g <- function(z) {
    g1_tilde(z)^2/g_tilde(z)  
  }
  integ <- NULL
  for (i  in 1:n) {
    xi <- x[i]
    ter <- integrate(g,0,xi)$value
    integ <- c(integ,ter)
  }
  beta_tilde <- (1/N_tilde)*sum((G_tilde(x)^(-2)))
  V23_tilde  <- (1/N_tilde)*sum(G1_tilde(x)/(G_tilde(x)^2))
  V33_tilde  <- (1/N_tilde)*(sum((G1_tilde(x)^2)/(G_tilde(x)^2))-sum(integ/G_tilde(x)))
  return(beta_tilde-1-V23_tilde*(1/V33_tilde)*V23_tilde)
}
#####################################################################################
#########   calculate Coverage percentage for condition likeligood
##   INPUT:
##         K :  Number of simulations
##         sigma_tilde :   the condition estimate of sigma,is a vector in p_c
##         N_tilde :  the condition estimate of N,is a vector in p_c
##         N_true  :  true n_big
##          a      :  Significance level
##
##   OUTPUT:
##         P_c :   Two type Coverage percentage for condition likelihood
##
######################################################
p_c  <- function(K,sigma_tilde,N_tilde,N_true,a){
  I2 <- ifelse(abs(N_tilde^(-1/2)*(N_tilde-N_true))<=sqrt(sigma_tilde)*qnorm(1-a/2),1,0)
  I3 <- ifelse(abs(N_tilde^(1/2)*log(N_tilde/N_true))<=sqrt(sigma_tilde)*qnorm(1-a/2),1,0)
  list(p2=sum(I2)/K,p3=sum(I3)/K)
}

2，EM算法

#######################################################################################################
#########################     EM Algorithm  for full likelihood  #####################################
#############################################################################################
####   k1 : The number of  unique values of  observed x 
####   par =[p1,p2,...,pk1,theta]: (k1+1) dim vector,initial value of parameter
####   ep=1e-5:precision
####   it_max=1e4:The number of maximum iteration
##################################################################################
##########################  Iterative process of EM algorithm   ###################
EM_it <- function(obs){
  
  n = nrow(obs)
  x = obs[,1]
  y = obs[,2]
  tilde_x = unique(x)
  k1 = length(tilde_x)
  p_0 = 1/k1
  theta_0 = sum(y)/n
  alpha_0 <- sum(pexp(tilde_x,theta_0)*p_0)
  par_0 = c(rep(p_0,k1),alpha_0,theta_0)
  
  #indicator <- function(x,t){ifelse(x==t,1,0)}
  # expoential density function
  g <- function(x,theta){
    theta*exp(-1*theta*x)
  }
  my.em <- function(obs,par=par_0,ep=0.001,it_max=50000){
    n = nrow(obs)
    x = obs[,1]
    y = obs[,2]
    tilde_x = unique(x)
    k1 = length(tilde_x)
    
    it_k  <- 0
    while(it_k <= it_max){
      
      par0 <- par
      
      ### iteration
      a <- NULL
      for (j in 1:k1) {
        aj <- sum(ifelse(x==tilde_x[j],1,0)) +(n/par0[k1+1])*par0[j]*(1-pexp(tilde_x[j],par0[k1+2]))
        a <- c(a,aj)
      }
      fn <- function(theta){
        b1 = sum(log_n(dexp(y,theta),n))
        
        log_dexp <- function(x){
          (log(theta)-theta*x)*par0[k1+2]*exp(-1*par0[k1+2]*x)
        }
        
        intb <- NULL
        for (j in 1:k1) {
          intbj = integrate(f = log_dexp, tilde_x[j],Inf)$value
          intb = c(intb,intbj)
        }
        b2 = (n/par0[k1+1])*sum(par0[1:k1]*intb)
        
        return(b1+b2)
      }
      
      p_1 <- a/sum(a)
      theta_1 <- optimize(f=fn,interval = c(0.0001,100),tol = 0.0001,maximum = TRUE)$maximum
      alpha_1 <- sum(pexp(x,theta_1)*p_1)
      par_1  <- c(p_1,alpha_1,theta_1)
      par <- par_1
      
      
      ## Whether to stop iteration
      norm <- sqrt(t(par-par0)%*%(par-par0))
      if(norm < ep)
      {
        return(par)
        break
      }
      else
      {it_k <- it_k+1}
    }     # while
  }#my.em
  par    <- my.em(obs,par=par_0,ep=0.001,it_max=50000)
  p_em   <- par[1:k1]
  alp_em <- par[k1+1]
  the_em <- par[k1+2]
  N_em   <- n/alp_em
  
  par_em <- data.frame(N_em=N_em,alp_em=alp_em,the_em)
  
  #return(list(p_em=p_hat,the_alp_em= the_alp_hat))
  return(par_em)
}#EM.it
## EM算法无法计算的覆盖率

3，解方程的方法

#########################################################################
##################### the method of solving the equation ################
See <- function(obs){
  n <- nrow(obs)
  x <- obs[,1]
  y <- obs[,2]
  ## expoential density function
  g <- function(x,theta){
    ifelse(x>=0,theta/exp(theta*x),0)
  }
  ## The first derivative of theta (g) 
  g1 <- function(x,theta){
    ifelse(x>=0,(1-theta*x)/exp(theta*x),0)
  }
  ## expoential distribution function
  G <-function(x,theta){
    ifelse(x>=0,1-1/exp(theta*x),0)
  }
   ## the first derivative of theta (G)
  G1 <- function(x,theta){
    ifelse(x>=0,x/exp(theta*x),0)
  }
  
  ### parameter is a vector par=(N,alpha,theta,lambda)
  model <- function(par){
    c(F1 = sum(1/((par[1]-n)+(1:n)))+log_n(1-par[2],n),
      F2 = (par[1]-n)/(1-par[2])- par[4]*sum(1/(1+par[4]*(G(x,par[3])-par[2]))),
      F3 = sum(g1(y,par[3])/g(y,par[3])) - par[4]*sum(G1(x,par[3])/(1+par[4]*(G(x,par[3])-par[2]))),
      F4 = sum((G(x,par[3])-par[2])/(1+par[4]*(G(x,par[3])-par[2])))
    )
  }
  par0 <- c(n+10,0.5,1,1) # the initial value of par，affect the simulating results
  par_See <- multiroot(f = model,start = par0,rtol = 1e-8,atol = 1e-8, ctol = 1e-8)$root
  return(par_See)
}

4.1，构建截面似然（内部使用解方程的方法：multiroot函数）

###########################################################################
###################  the fourth method ###################################
profile_max.like1 <- function(obs){
  n <- nrow(obs)
  x <- obs[,1]
  y <- obs[,2]
  ## expoential density function
  g <- function(x,theta){
    ifelse(x>=0,theta/exp(theta*x),0)
  }
  ## The first derivative of theta (g) 
  g1 <- function(x,theta){
    ifelse(x>=0,(1-theta*x)/exp(theta*x),0)
  }
  ## expoential distribution function
  G <-function(x,theta){
    ifelse(x>=0,1-1/exp(theta*x),0)
  }
  ## the first derivative of theta (G)
  G1 <- function(x,theta){
    ifelse(x>=0,x/exp(theta*x),0)
  }
  
  ### the profile likelihood of N by solving the equation about (alpha,theta,lambda)
  f <- function(N){
    ## par is a vector, par = (alpha,theta,lambda)
    model0 <- function(par){
      c(
        F1 = (N-n)/(1-par[1])- par[3]*sum(1/(1+par[3]*(G(x,par[2])-par[1]))),
        F2 = sum(g1(y,par[2])/g(y,par[2])) - par[3]*sum(G1(x,par[2])/(1+par[3]*(G(x,par[2])-par[1]))),
        F3 = sum((G(x,par[2])-par[1])/(1+par[3]*(G(x,par[2])-par[1])))  
      )
    }
    par0 <- c(0.5,1,1) # the initial value of par
    atl <- multiroot(f = model0,start = par0, rtol = 1e-8,atol = 1e-8, ctol = 1e-8)$root # the root of (alpha,theta,lambda)
    sum(log_n((N-n)+(1:n),n))+(N-n)*log_n(1-atl[1],n)+sum(log_n(g(y,atl[2]),n))-sum(log_n(1+atl[3]*(G(x,atl[2])-atl[1]),n))
  }
  N_est <- optimize(f, interval=c(n,400), tol=0.001,maximum = TRUE)$objective # the estimate of N
  
  ## par is a vector, par = (alpha,theta,lambda)
  model <- function(par){
    c(
      F1 = (N_est-n)/(1-par[1])- par[3]*sum(1/(1+par[3]*(G(x,par[2])-par[1]))),
      F2 = sum(g1(y,par[2])/g(y,par[2])) - par[3]*sum(G1(x,par[2])/(1+par[3]*(G(x,par[2])-par[1]))),
      F3 = sum((G(x,par[2])-par[1])/(1+par[3]*(G(x,par[2])-par[1])))  
    )
  }
  par0 <- c(0.5,1,1) # the initial value of par
  atl_est <- multiroot(f = model,start = par0, rtol = 1e-8,atol = 1e-8, ctol = 1e-8)$root 
  ## profile_par is a vector,profile_par.est=[N_est,alpha_est,theta_est,lambda_est]
  profile_par.est <- c(N_est,atl_est) 
  return(profile_par.est)
}

4.2，构建截面似然（内部使用极大化似然的方法：nlminb函数）

############## Using the function nlminb to solve equations ######################
############################################################################
profile_max.like2 <- function(obs){
  n <- nrow(obs)
  x <- obs[,1]
  y <- obs[,2]
  ## expoential density function
  g <- function(x,theta){
    ifelse(x>=0,theta/exp(theta*x),0)
  }
  ## expoential distribution function
  G <-function(x,theta){
    ifelse(x>=0,1-1/exp(theta*x),0)
  }
  ############################################
  ###
  ##  f is the profile likelihood of N  
  ##  par is a vector, par = (alpha,theta,lambda)
  ####################################################
  f <- function(N) { 
    fun <- function(par){
      
      (-1)*(sum(log_n((N-n)+(1:n),n))+(N-n)*log_n(1-par[1],n)+sum(log_n(g(y,par[2]),n))-sum(log_n(1+par[3]*(G(x,par[2])-par[1]),n)))
      
    }
    
    atl <- nlminb(start = c(0.5,5,5),objective = fun,lower = c(0.001,0.001,-500),upper = c(0.9999,100,500))$par
    sum(log_n((N-n)+(1:n),n))+(N-n)*log_n(1-atl[1],n)+sum(log_n(g(y,atl[2]),n))-sum(log_n(1+atl[3]*(G(x,atl[2])-atl[1]),n))
  }
  N_est <- optimize(f, interval=c(n,500), tol=0.001,maximum = TRUE)$objective 
  ### Replace N with N_est
  fun <- function(par){
    (-1)*(sum(log_n((N_est-n)+(1:n),n))+(N_est-n)*log_n(1-par[1],n)+sum(log_n(g(y,par[2]),n))-sum(log_n(1+par[3]*(G(x,par[2])-par[1]),n)))
  }
  atl_est <- nlminb(start = c(0.5,5,5),objective = fun,lower = c(0.001,0.001,-500),upper = c(0.9999,100,500))$par
  ## par_full.est is a vector: par_full.est=[N_est,alpha_est,theta_est,lambda_est]
  par_full.est <- c(N_est,atl_est)
  return(par_full.est)
}

5，四种模拟结果

5.1模拟参数设定

#################################################################################
####             Simulation in one scenario
#################################################################################
N_true <- 200   
K <- 2000      #repeat times
a <- 0.05       #significance level
wshape <- 5
wscale <- 2
erate  <- 2   # true theta
theta_true <- erate
comf <- function(x){
  pexp(x,rate = erate)*dweibull(x,shape = wshape,scale = wscale)
}
alpha_true <- integrate(f=comf,lower = 0,upper = Inf)$value;alpha_true  # true alpha
par_true <- rep(c(N_true,alpha_true,theta_true),2);par_true

此处 的N取值为200，因此前面max.like2函数中N0取值应更改为200. 为考虑不同的入样概率，选择以下几种情形的参数设置( α \alpha α是入样概率)：

参数wshapewscaleerate（ θ \theta θ） α \alpha α
设定15220.963
设定21220.8
设定31210.667
设定41110.5

5.2模拟结果以及时间

所有的K=2000,N=200.

1，直接极大似然方法
#######################  simulating Result ################
### 方法1
par_est =NULL     
N_hat   = NULL
N_tilde = NULL     
sigma_tilde = NULL    #条件似然方法的sigma的估计
int <- NULL
t0 <- proc.time()
for (i in 1:K) {
  obs <-  gen_obs(wshape, wscale, erate, N_true)  
  par_est <- rbind(par_est, c(max.like(obs), max.clike(obs)))
  ter_para_star <- c(max.like(obs), max.like2(obs),N_true,1-a)
  N_hat   <- rbind(N_hat,max.like(obs)[1])
  N_tilde <-  rbind(N_tilde,max.clike(obs)[1])
  sigma_tilde <- rbind(sigma_tilde,sigma_cest(obs = obs,theta_tilde = max.clike(obs)[3],N_tilde = max.clike(obs)[1]))
  ter_int <- con_int(obs = obs,para_star = ter_para_star,N_true = N_true)
  int <- rbind(int,ter_int)
  print(i) 
}
t <- proc.time()-t0
print(paste0("直接极大似然方法执行时间：",t[3][[1]],"秒"))

par_true <- rep(c(N_true,alpha_true,theta_true),2);par_true
m    <-  apply(par_est, 2, mean);m
sd   <-  apply(par_est, 2, sd);sd
bias <-  abs(apply(par_est,2,mean) - par_true);bias
percentage_f <- p_f(K,int = int);percentage_f 
percentage_c <- p_c(K,sigma_tilde,N_tilde,N_true,a);percentage_c

模拟结果表

模拟结果表中依次为参数估计，参数估计标准差和偏差。

参数估计N_f α \alpha α_f θ \theta θ_fpercentage_fN_c α \alpha α_c θ \theta θ_cpercentage_c2percentage_c3运行时间
设定1199.3616 3.6174 0.63840.9635 0.0108 0.00042.0320 0.1888 0.032094.05%200.1870 3.6645 0.18700.9620 0.0111 0.00112.0065 0.1869 0.006588.15%87.9%1416秒
设定2198.2373 17.6154 1.76270.8093 0.0533 0.00932.0228 0.2109 0.022891%200.2741 19.3338 0.27410.8037 0.0563 0.00372.0092 0.2097 0.009285.8%86.1%1424.87秒
设定3198.9731 32.7968 1.02690.6793 0.0749 0.01271.0237 0.1430 0.023791.7%201.5994 34.7736 1.59940.6730 0.0764 0.00631.0139 0.1420 0.013983.05%83.4%1182.45秒
设定4201.2625 45.9229 1.26250.5156 0.0918 0.01561.0285 0.2315 0.028593.75%205.3044 47.6765 5.30440.5074 0.0920 0.00741.0104 0.2299 0.010477.35%76.3%918.92秒
2，EM方法
par_em <- NULL
t0 <- proc.time()
for (i in 1:K) {
  obs <-  gen_obs(wshape, wscale, erate, N_true)  
  ter_em <- EM_it(obs)
  par_em <- rbind(par_em,ter_em)
  print(i) 
}
t <- proc.time()-t0
print(paste0("EM方法执行时间：",t[3][[1]],"秒"))
par_true <- rep(c(N_true,alpha_true,theta_true),2);par_true
m_em    <-  apply(par_em, 2, mean);m_em
sd_em   <-  apply(par_em, 2, sd);sd_em
bias_em <-  abs(apply(par_em,2,mean) - par_true[1:3]);bias_em

EM方法模拟结果表：
结果依次为参数估计、参数估计标准差和偏差。

参数估计N_em α \alpha α_em θ \theta θ_em运行时间
设定1200.2807 3.5829 0.28070.9623 0.0108 0.00082.0119 0.1832 0.01193809.43秒
设定2198.9456 14.6940 1.05440.8070 0.0492 0.00702.0127 0.2060 0.01273711.39秒
设定3199.8875 25.5943 0.11250.6768 0.0688 0.01011.0105 0.1425 0.01053803.96秒
设定4202.7574 44.7711 2.75740.5103 0.0916 0.01031.0161 0.2407 0.01617032.74秒
3，解方程方法
par_See <- NULL # par_See = [N_est,alpha_est,theta_est,lambda_est]
t0 <- proc.time()
for (i in 1:K) {
  obs <-  gen_obs(wshape, wscale, erate, N_true)
  par_See <- rbind(par_See,See(obs))
  print(i)
}
t <- proc.time()-t0
print(paste0("解方程方法执行时间：",t[3][[1]],"秒"))
par_true <- rep(c(N_true,alpha_true,theta_true),2);par_true
m_See    <-  apply(par_See[,1:3], 2, mean);m_See
sd_See   <-  apply(par_See[,1:3], 2, sd);sd_See
bias_See <-  abs(apply(par_See[,1:3],2,mean) - par_true[1:3]);bias_See

程序运行中出现警告，

4，构建截面似然方程方法
# 内部解方程方法
par_profile1 <- NULL # par_profile1 =[N_est,alpha_est,theta_est,lambda_est]
t0 <- proc.time()
for (i in 1:K) {
  obs <-  gen_obs(wshape, wscale, erate, N_true)
  par_profile1 <- rbind(par_profile1,profile_max.like1(obs))
  print(i)
}
t <- proc.time()-t0
print(paste0("构建截面似然(内部解方程)方法执行时间：",t[3][[1]],"秒"))
par_true <- rep(c(N_true,alpha_true,theta_true),2);par_true
m_profile1    <-  apply(par_profile1[,1:3], 2, mean);m_profile1
sd_profile1   <-  apply(par_profile1[,1:3], 2, sd);sd_profile1
bias_profile1 <-  abs(apply(par_profile1[,1:3],2,mean) - par_true[1:3]);bias_profile1

# 内部直接极大似然(nlminb)
par_profile2 <- NULL # par_profile2 =[N_est,alpha_est,theta_est,lambda_est]
t0 <- proc.time()
for (i in 1:K) {
  obs <-  gen_obs(wshape, wscale, erate, N_true)
  par_profile2 <- rbind(par_profile2,profile_max.like2(obs))
  print(i)
}
t <- proc.time()-t0
print(paste0("构建截面似然(内部解方程)方法执行时间：",t[3][[1]],"秒"))
par_true <- rep(c(N_true,alpha_true,theta_true),2);par_true
m_profile2    <-  apply(par_profile2[,1:3], 2, mean);m_profile2
sd_profile2   <-  apply(par_profile2[,1:3], 2, sd);sd_profile2
bias_profile2 <-  abs(apply(par_profile2[,1:3],2,mean) - par_true[1:3]);bias_profile2

最后两种求解方法速度特别快，但是用到的R函数求解出来的结果都有一些问题，具体问题见文档。

Cet777
关注
基于朴素贝叶斯分类器的 20-news-group分类及结果对比(Python3)
weixin_39831037的博客
10-18 2687
之前看了很多CSDN文章，很多都是根据stack overflow 或者一些英文网站的照搬。导致我看了一整天最后一点收获都没有。 这个作业也借鉴了很多外文网站的帮助 但是是基于自己理解写的。环境是python3(海外留学 This Task is going to produce a text classification for 20-newgroups corpus. It will be i...
Multiboot 规范
uframer(焦冶)的专栏
05-12 9642
Multiboot规范Multiboot规范本文定义了Multiboot规范——提议中的引导过程标准。本文是此规范的0.6.93版。 概览 术语 规范 示例 历史 索引节点: 概览, 前进: 术语, 后退: 顶层, 向上: 顶层 Multiboot规范简介本章描述了一些关于Multiboot规范的粗略的信息。注意，这不是规范本身的一部分。 动机 架构 操作系统 引导源 引导时配置 为操作系统带来的
《统计建模与R软件》第四章笔记
最新发布
qq_69383479的博客
12-23 1435
参数估计是统计建模的关键步骤之一，它涉及根据样本数据推断总体参数的过程。在统计学中，参数通常用于描述总体的特征，如均值、方差等。通过参数估计，我们可以利用样本信息对这些未知参数进行推断，从而对总体进行更深入的了解。思想：当我们面对一个统计问题时，通常我们不能观察到整个总体的所有数据，而只能通过取一部分样本来进行研究。为了从这个样本中了解总体的性质，我们引入了一种思想，即使用样本的一些数字特征（矩）来估计总体相应的特征。在这个过程中，我们关注的是总体的矩，而这些矩与总体的参数有密切的关系，从而允许我们得出对总
GSL中的多维求根
yxmlhc的专栏
10-03 389
多维求根 本章介绍多维求根（用n个未知数中的n个方程求解非线性系统）的函数。 本库为各种迭代求解器和收敛测试提供了低级组件。用户可以将它们组合起来以实现所需的解决方案，同时可以完全访问迭代的中间步骤。每类方法都使用相同的框架，因此您可以在运行时在求解器之间切换，而无需重新编译程序。求解器的每个实例都会跟踪自己的状态，从而允许将求解器用于多线程程序中。求解器是基于原始的Fortran库MINPACK。 头文件gsl_multiroot.h包含多维求根函数和相关声明的原型。 38.1 概述...
【数据分析 R语言实战】学习笔记 第六章 参数估计与R实现（上）
05-22 4588
6.1点估计及R实现 6.1.1矩估计 R中的解方程函数: 函数及所在包：功能 uniroot()@stats：求解一元（非线性)方程 multiroot()@rootSolve：给定n个(非线性)方程，求解n个根 uniroot.all()@rootSolve：在一个区问内求解一个方程的多个根 BBsolve()@BB：使用Barzilai-Borwein步长求解非线性方程组 u...
R语言学习笔记（四）参数估计
zorchp
07-04 1万+
总结一下数理统计中的参数估计，即点估计(矩估计、极大似然估计)和区间估计(置信区间)部分的R语言实现
四种TEK数字示波器性能的比较
01-20
摘要： 对物理实验中常用的TEK公司的四种型号数字示波器延迟时间抖动进行测量， 给出的四种T E K 数字示波器延迟时间差异， 有助于设计实验测试系统时正确选型， 并作为测试系统中的一部分修正量对实验测量结果进行...
基础电子中的四种TEK数字示波器性能的比较
10-21
摘要： 对物理实验中常用的TEK公司的四种型号数字示波器延迟时间抖动进行测量， 给出的四种T E K 数字示波器延迟时间差异， 有助于设计实验测试系统时正确选型， 并作为测试系统中的一部分修正量对实验测量结果进行...
三种方法的比较
12-21
超市管理的目标基本以消费者对象而确定，提供的三种方法的比较能够给你学习和参考三种方法的比...该文档为三种方法的比较，是一份很不错的参考资料，具有较高参考价值，感兴趣的可以下载看看
Java多线程实现四种方式原理详解
12-17
前面两种可以归结为一类：无返回值，原因很简单，通过重写run方法，run方式的返回值是void，所以没有办法返回结果 后面两种可以归结成一类：有返回值，通过Callable接口，就要实现call方法，这个方法的返回值是...
多方法比较
lSky_Forever的博客
09-28 493
9*9乘法表 public class  test{ public static void main(String[] args) { for (int x=1;x 以x和y的初始值要设成1 { for (int y=1;y { System.out.print(y+"×"+x+"="+x*y+"\t");//从 行上面看，第一个数是变动的，而且是每一列加
R学习之统计算法与R优化包（极大似然估计） --（R语言编程）-----数模
热门推荐
默默小艺的博客
08-06 1万+
library(MASS);attach(geyser);hist(waiting,freq = F) # mnf&amp;amp;amp;lt;-function(pa,data){ x&amp;amp;amp;lt;-dnorm(data,pa[2],sqrt(pa[4])) y&amp;amp;amp;lt;-dnorm(data,pa[3],sqrt(pa[5])) pdf=pa[1]*x+(1-pa[1])*y l=sum(log(...
如何分析软件测试结果,软件测试结果对比
weixin_42493720的博客
07-24 818
软件测试结果对比测试结果对比为防止软件测试成绩不稳定，每款软件都经过三次测试，选择测试分数排名中间的成绩进行对比。对比测试之前，笔者也收到了蓝魔W17PRO的最新固件，处理器主频可提升至1.5GHz，并且针对性能方面有了进一步优化，笔者也在第一时间升级并进行软件测试。三款机型的各款软件测试结果对比如下所示。Quadrant测试成绩对比 Antutu测试成绩对比 图Nenamark1测试成绩对比 通...
Multiboot规范（2）
wuhui_gdnt的专栏
07-31 2865
4.      例子 警告：以下内容不是该规范文档的一部分，但为了未来的操作系统及启动载入器作者，而包括进来。 4.1        PC注意事项 考虑在Multiboot信息结构体中‘flags’参数的第0位，如果讨论中的启动载入器使用旧的BIOS接口，或最新的不可用（参
R语言中优化函数
辉的博客
04-01 1万+
利用Rdonlp2包求解有约束优化问题 安装代码： install.packages("Rdonlp2", repos="http://R-Forge.R-project.org") 使用可参考： https://blog.csdn.net/qq_27755195/article/details/51458659 ...
java多线程同步的四种方法_java中实现多线程的两种方法
weixin_36409041的博客
02-24 1883
java多线程有几种实现方法,都是什么?同步有几种实java中多线程的实现方法有两种：1.直接继承thread类；2.实现runnable接口；同步的实现方法有五种：1.同步方法；2.同步代码块；3.使用特殊域变量(volatile)实现线程同步；4.使用重入锁实现线程同步；5.使用局部变量实现线程同步 。JAVA多线程有哪几种实现方式JAVA多线程实现方式主要有三种：继承Thread类、实现Ru...
4种比较方法的示例
一心一信
06-17 431
见图：
Conditional Likelihood Maximisation: A Unifying Framework for Information Theoretic Feature Selecti
tangxianyu的博客
10-30 1225
目录1、文章信息2、主要思想2.1信息熵：2.2 基于互信息的滤波算法 1、文章信息 Title: Conditional Likelihood Maximisation: A Unifying Framework for Information Theoretic Feature Selection **Author:**Gavin Brown, Adam Pocock, Ming-Jie Z...
R语言统计方法
11-07
R语言初级统计方法，里面的方法都比较的全面详细具体，挺不错的
R语言中PORT方法
10-25
R统计语言的PORT方法，是nlminb，optimize等函数的具体算法来源
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值