估计一些生成随机数的分布参数

Question

我刚开始使用 R 来解决我的统计问题。目前，我正在使用我使用 R 生成的 200 个随机数 (RN) 来估计分布的参数。我在 100 次中生成了 200 个 RN。所以就是说会有100种200个RN我估计这100种RN。这也意味着将有100种估算结果。所以这是我用来生成 RN 的代码：

#Generate random numbers U~(0, 1)
rep <-100 #total replication
unif <-matrix(0, 200, rep)
for (k in 1: rep)
{
  unif[,k] <- runif(200, min = 0, max = 1)
}

# Based on the 100 kinds of generated random numbers that follow U ~ (0.1), I will generate again 100 kinds of random numbers which follow the estimated distribution:
# Define parameters
a <- 49.05 #1st parameter
b <- 3.148 #2nd parameter
c <- 0.145 #3rd parameter
d <- 0.00007181 #4th parameter

X <-matrix(0, 200, rep)
for (k in 1: rep)
{
  X[,k] <- a*(log(1-((log(1-((unif[,k])^(1/c))))/(a*d))))^(1/b)
}

# Sorting the generated RN from the smallest to the largest
X_sort <-matrix(0, 200, rep)
for (k in 1: rep)
{
  X_sort[,k] <- sort(X[,k])
}

到这里我已经成功地生成了 100 种将被估计的 RN。但是，我现在面临的问题是如何估算这100种RN。我只能估计一个。这是我使用 maxLik 包估计参数的代码，估计方法是 BHHH:

xi = X_sort[,1]
log_likelihood<-function(theta,xi){
  p1 <- theta[1] #1st parameter
  p2 <- theta[2] #2nd parameter
  p3 <- theta[3] #3rd parameter
  p4 <- theta[4] #4th parameter
  logL=log((p4*p2*p3*((xi/p1)^(p2-1))*(exp(((xi/p1)^(p2))+(p4*p1*(1-(exp((xi/p1)^(p2)))))))*((1-(exp((p4*p1*(1-(exp((xi/p1)^(p2))))))))^(p3-1))))
  return(logL)
}
library(maxLik);
# Initial parameters
a <- 49.05 #1st parameter
b <- 3.148 #2nd parameter
c <- 0.145 #3rd parameter
d <- 0.00007181 #4th parameter
m <- maxLik(log_likelihood, start=c(a,b,c,d), xi = xi, method="bhhh");
summary(m)

结果如下：

--------------------------------------------
Maximum Likelihood estimation
BHHH maximisation, 5 iterations
Return code 2: successive function values within tolerance limit
Log-Likelihood: -874.0024 
4  free parameters
Estimates:
      Estimate Std. error t value  Pr(> t)    
[1,] 4.790e+01  1.846e+00  25.953  < 2e-16 ***
[2,] 3.015e+00  1.252e-01  24.091  < 2e-16 ***
[3,] 1.717e-01  2.964e-02   5.793 6.91e-09 ***
[4,] 7.751e-05  6.909e-05   1.122    0.262    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------

要估计另外99个RN，我得手动改成xi = X_sort[,k] for k=1,2,...,100 所以第二个RN应该变成X_sort[,2] ，依此类推，直到第 100 个 RN。我认为这样效率不高，因为一个一个替换它们需要很长时间。那么有没有办法修改这段代码，让估计其他RN不用花很长时间呢？

Answer 1

首先，我建议您以更紧凑的方式重写代码。

1.生成随机数。 不需要生成100个长度为200的向量，我们可以生成一个长度为100*200的向量，然后将这个向量逐列写入矩阵。这可以通过以下方式完成：

rep <-100
n <-  200
unif <- matrix(runif(rep*n, min = 0, max = 1), n, rep)

2。矩阵的计算函数。 在 R 中，可以将矢量函数应用于矩阵。所以在你的情况下它将是：

X <- a*(log(1-((log(1-((unif)^(1/c))))/(a*d))))^(1/b)

3。按列矩阵排序 我们可以使用apply 函数轻松对矩阵的每一列进行排序。参数2表示我们按列进行（1代表行）。

X_sort <- apply(X, 2, sort)

4.执行估计。 同样，我们可以在此处使用 apply。

estimations <- apply(X_sort, 2, function(x) maxLik(log_likelihood, start=c(a,b,c,d),
                                          xi = x, method="bhhh"))

然后要打印所有摘要，您可以执行以下操作：

lapply(estimations, summary)

估计一些生成随机数的分布参数

Estimating parameter of a distribution for some generate random number

parameters

statistics

estimation

r

distribution