乱数を生成するための分布のパラメータを推定する

Question

私の統計的な問題を解決するためのRの新機能です。現在、私はRを使って生成する200の乱数（RN）を使って分布のパラメータを推定するように働いています。私は200回のRNを100回生成します。 200種類のRNが100種類あり、この100種類のRNが推定されます。 100種類の推定結果があることも意味します。

#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]) 
} 

ここ

アップ私が推定されるRNの100種類を生成することができた： だからここで私はRNを生成するために使用するコードです。しかし、今私が直面している問題は、この100種類のRNをどのように見積もるかです。私は1つしか推定することができません。ここで私は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とすると、2番目のRNはX_sort[,2]になり、100番目のRNまで続きます。私はこれを効率的ではないと思います。なぜなら、それを一つずつ置き換えるのに時間がかかるからです。それで、このコードを修正して他のRNを推定するのに時間がかからないようにする方法はありますか？

Answer 1

まず、コードをよりコンパクトに書き直すことをお勧めします。

1.乱数を生成する。長さ200のベクトルをそれぞれ100個生成する必要はありませんが、長さ100 * 200のベクトルを生成し、このベクトルを列方向に行列に書き込むことができます。これは、次の方法で行うことができます

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

2.行列の機能を計算します。 Rではベクトル関数を行列に適用することができます。だからあなたの場合には、それは次のようになります。私たちは、簡単にapply機能を使用して行列の各列を並べ替えることができをソート

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

3列単位行列。パラメータ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)