私の統計的な問題を解決するための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を推定するのに時間がかからないようにする方法はありますか?
はあなたの助けのためにありがとうございました。これは本当に私の問題を解決する – crhburn