Оценка параметра распределения для некоторого порождающего случайного числа

Question

Я просто новый в R для решения моей статистической задачи. В настоящее время я работаю над оценкой параметров распределения с использованием 200 случайных чисел (RN), которые я генерирую с использованием R. Я генерирую 200 RN в 100 раз. Таким образом, это означает, что будет 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]) 
} 

Up здесь я сумел создать 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] для k = 1,2, ..., 100, поэтому для второго RN он должен превратиться в X_sort[,2] и так далее до сотого 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)