Estimation du paramètre d'une distribution pour certains générer un nombre aléatoire

Question

Je viens juste de commencer par R pour résoudre mon problème statistique. Actuellement, je travaille pour estimer les paramètres d'une distribution en utilisant 200 nombres aléatoires (RN) que je génère en utilisant R. Je génère 200 RN dans 100 fois. Donc, cela signifie qu'il y aura 100 sortes de 200 RN et je vais estimer ces 100 types de RN. Cela signifie également qu'il y aura 100 types de résultats d'estimation. Voici donc le code que j'utilise pour générer la 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]) 
} 

Jusqu'à ici, j'ai réussi à générer 100 types de RN qui seront estimés. Cependant, le problème auquel je fais face maintenant est de savoir comment estimer ces 100 types de RN. Je peux seulement en estimer un. Voici le code que j'utilise pour estimer le paramètre avec le package maxLik et la méthode d'estimation est 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) 

Voici le résultat:

-------------------------------------------- 
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 
-------------------------------------------- 

Pour estimer l'autre 99 RN, je dois changer manuellement xi = X_sort[,k] pour k = 1,2, ..., 100, donc pour le deuxième RN, il devrait se transformer en X_sort[,2], et ainsi de suite jusqu'à la centième RN. Je pense que ce n'est pas efficace parce qu'il faut du temps pour les remplacer un par un. Donc, y a-t-il un moyen de modifier ce code pour qu'il ne soit pas long à estimer les autres RN?

Answer 1

Tout d'abord, je vous suggère de réécrire votre code de manière plus compacte.

1. Génération de nombres aléatoires. Il n'est pas nécessaire de générer 100 vecteurs de longueur 200 chacun, alors que nous pouvons générer un vecteur de longueur 100 * 200, puis écrire ce vecteur en colonne dans une matrice. Cela peut se faire de la manière suivante:

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

2. Calcul de la fonction de la matrice. En R, il est possible d'appliquer des fonctions vectorielles aux matrices. Donc, dans votre cas, il sera:

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

3. matrice en colonne de tri On peut facilement trier chaque colonne de la matrice en utilisant apply fonction. Le paramètre 2 signifie que nous le faisons par colonne (1 correspond aux rangées).

X_sort <- apply(X, 2, sort) 

4. Effectuer des estimations. Encore une fois, nous pouvons utiliser apply ici.

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

ensuite pour imprimer tous les résumés, vous pouvez effectuer les opérations suivantes:

lapply(estimations, summary)