Ajuster la distribution Gamma aux données en R en utilisant optim, ML

Question

Je suis un peu nouveau pour R. J'ai un jeu de données, qui inclut également des données de revenu familial et je dois ajuster une distribution gamma à ces données, en utilisant les estimations de vraisemblance maximum. Il est spécifiquement dit que nous avons besoin d'utiliser le paquet optim, et non fitdistr. Donc ceci est mon code:

t1 <- sum(log(newdata$faminc)) 
t2 <- sum(newdata$faminc) 
obs <- nrow(newdata) 
lh.gamma <- function(par) { 
    -((par[1]-1)*t1 - par[2]*t2 - obs*par[1]*log(par[2]) - obs*lgamma(par[1])) 
} 

#initial guess for a = mean^2(x)/var(x) and b = mean(x)/var(x) 
a1 <- (mean(newdata$faminc))^2/var(newdata$faminc) 
b1 <- mean(newdata$faminc)/var(newdata$faminc) 

init <- c(a1,b1) 
q <- optim(init, lh.gamma, method = "BFGS") 
q 

Vous avez également essayé de remplir juste des valeurs dans le vecteur init, et d'inclure cette partie de code;

dlh.gamma <- function(par){ 
    cbind(obs*digamma(par[1])+obs*log(par[2])-t2, 
    obs*par[1]/par[2]-1/par[2]^2*t1) 
} 

et le optim ressemblerait à ceci:

q <- optim(init, lh.gamma, dhl.gamma, method="BFGS") 

Rien de tout cela 'fonctionne'. D'abord, quand j'ai essayé le code sur les ordinateurs de l'école, cela m'a donné de très gros nombres pour les paramètres de forme et de vitesse, ce qui n'était pas possible. Maintenant, en essayant à la maison, je reçois ceci:

> q <- optim(init, lh.gamma, method = "BFGS") 
Error in optim(init, lh.gamma, method = "BFGS") : 
    non-finite finite-difference value [2] 
In addition: There were 50 or more warnings (use warnings() to see the first 50) 
> q 
function (save = "default", status = 0, runLast = TRUE) 
.Internal(quit(save, status, runLast)) 
<bytecode: 0x000000000eaac960> 
<environment: namespace:base> 

q n'est même pas 'créé'. Sauf quand j'inclus la partie dlh.gamma ci-dessus, mais ensuite je reçois juste des nombres énormes et pas de convergence.

Toute personne qui sait ce qui ne va pas ou que faire?

Edit:

> dput(sample(newdata$faminc, 500)) 
c(42.5, 87.5, 22.5, 17.5, 12.5, 30, 30, 17.5, 42.5, 62.5, 62.5, 
30, 30, 150, 22.5, 30, 42.5, 30, 17.5, 8.75, 42.5, 42.5, 42.5, 
62.5, 42.5, 30, 17.5, 87.5, 62.5, 150, 42.5, 150, 42.5, 42.5, 
42.5, 6.25, 62.5, 87.5, 6.25, 87.5, 30, 150, 22.5, 62.5, 42.5,  
150, 17.5, 42.5, 42.5, 42.5, 62.5, 22.5, 42.5, 42.5, 30, 62.5, 
30, 62.5, 87.5, 87.5, 42.5, 22.5, 62.5, 22.5, 8.75, 30, 30, 17.5, 
87.5, 8.75, 62.5, 30, 17.5, 22.5, 62.5, 42.5, 30, 17.5, 62.5, 
8.75, 62.5, 42.5, 150, 30, 62.5, 87.5, 17.5, 62.5, 30, 62.5, 
87.5, 42.5, 62.5, 30, 62.5, 42.5, 87.5, 150, 12.5, 42.5, 62.5, 
42.5, 62.5, 62.5, 150, 30, 87.5, 12.5, 17.5, 42.5, 62.5, 30, 
6.25, 62.5, 42.5, 12.5, 62.5, 8.75, 17.5, 42.5, 62.5, 87.5, 8.75, 
62.5, 30, 62.5, 87.5, 42.5, 62.5, 62.5, 12.5, 150, 42.5, 62.5, 
12.5, 62.5, 42.5, 62.5, 62.5, 87.5, 42.5, 62.5, 30, 42.5, 150, 
42.5, 30, 62.5, 62.5, 87.5, 42.5, 30, 62.5, 62.5, 42.5, 42.5, 
30, 62.5, 42.5, 42.5, 62.5, 62.5, 150, 42.5, 30, 42.5, 62.5, 
17.5, 62.5, 17.5, 150, 8.75, 62.5, 30, 62.5, 42.5, 42.5, 22.5, 
150, 62.5, 42.5, 62.5, 62.5, 22.5, 30, 62.5, 30, 150, 42.5, 42.5, 
42.5, 62.5, 30, 12.5, 30, 150, 12.5, 8.75, 22.5, 30, 22.5, 30, 
42.5, 42.5, 42.5, 30, 12.5, 62.5, 42.5, 30, 22.5, 42.5, 87.5, 
22.5, 12.5, 42.5, 62.5, 62.5, 62.5, 30, 42.5, 30, 62.5, 30, 62.5, 
12.5, 22.5, 42.5, 22.5, 87.5, 30, 22.5, 17.5, 42.5, 62.5, 17.5, 
250, 150, 42.5, 30, 42.5, 30, 62.5, 17.5, 87.5, 22.5, 150, 62.5, 
42.5, 6.25, 87.5, 62.5, 42.5, 30, 42.5, 62.5, 42.5, 87.5, 62.5, 
150, 42.5, 30, 6.25, 22.5, 30, 42.5, 42.5, 62.5, 250, 8.75, 150, 
42.5, 30, 42.5, 30, 42.5, 42.5, 30, 30, 150, 22.5, 62.5, 30, 
8.75, 150, 62.5, 87.5, 150, 42.5, 30, 42.5, 42.5, 42.5, 30, 8.75, 
42.5, 42.5, 30, 22.5, 62.5, 17.5, 62.5, 62.5, 42.5, 8.75, 42.5, 
12.5, 12.5, 150, 42.5, 42.5, 17.5, 42.5, 62.5, 62.5, 42.5, 42.5, 
30, 42.5, 62.5, 30, 62.5, 42.5, 42.5, 42.5, 22.5, 62.5, 62.5, 
62.5, 22.5, 150, 62.5, 42.5, 62.5, 42.5, 30, 30, 62.5, 22.5, 
62.5, 87.5, 62.5, 42.5, 42.5, 22.5, 62.5, 62.5, 30, 42.5, 42.5, 
8.75, 87.5, 42.5, 42.5, 87.5, 30, 62.5, 17.5, 62.5, 42.5, 17.5, 
22.5, 62.5, 8.75, 62.5, 22.5, 22.5, 22.5, 42.5, 17.5, 22.5, 62.5, 
42.5, 42.5, 42.5, 42.5, 42.5, 30, 30, 8.75, 30, 42.5, 62.5, 22.5, 
6.25, 30, 42.5, 62.5, 17.5, 62.5, 42.5, 8.75, 22.5, 30, 17.5, 
22.5, 62.5, 42.5, 150, 87.5, 22.5, 12.5, 62.5, 62.5, 62.5, 30, 
42.5, 22.5, 62.5, 87.5, 30, 42.5, 62.5, 22.5, 87.5, 30, 30, 22.5, 
87.5, 87.5, 250, 30, 62.5, 250, 62.5, 42.5, 42.5, 62.5, 62.5, 
42.5, 6.25, 62.5, 62.5, 62.5, 42.5, 42.5, 150, 62.5, 62.5, 30, 
150, 22.5, 87.5, 30, 150, 17.5, 8.75, 62.5, 42.5, 62.5, 150, 
42.5, 22.5, 42.5, 42.5, 17.5, 62.5, 17.5, 62.5, 42.5, 150, 250, 
22.5, 42.5, 30, 62.5, 62.5, 42.5, 42.5, 30, 150, 150, 42.5, 17.5, 
17.5, 42.5, 8.75, 62.5, 42.5, 42.5, 22.5, 150, 62.5, 30, 250, 
62.5, 87.5, 62.5, 8.75, 62.5, 30, 30, 8.75, 17.5, 17.5, 150, 
22.5, 62.5, 62.5, 42.5) 

La variable faminc est en 1000s

Edit2:

D'accord, le code est bon, mais maintenant j'essayer d'adapter la distribution sur l'histogramme en utilisant les éléments suivants :

x <- rgamma(500,shape=q$par[1],scale=q$par[2]) 
hist(newdata$faminc, prob = TRUE) 
curve(dgamma(x, shape=q$par[1], scale=q$par[2]), add=TRUE, col='blue') 

Il produit simplement une ligne bleue plate à l'axe des abscisses.

Answer 1

Vous avez des choses en cours que je n'ai pas réussi à faire, mais voici une démonstration de l'estimation. Commençons par générer des données (nous savons donc si l'optimisation fonctionne). J'ai seulement changé votre fonction d'optimisation ci-dessous, et j'ai utilisé Nelder-Mead au lieu du quasi-Newton.

set.seed(23) 
a <- 2 # shape 
b <- 3 # rate 

require(data.table) 
newdata <- data.table(faminc = rgamma(10000, a, b)) 

t1 <- sum(log(newdata$faminc)) 
t2 <- sum(newdata$faminc) 
obs <- nrow(newdata) 

llf <- function(x){ 
    a <- x[1] 
    b <- x[2] 
    # log-likelihood function 
    return(- ((a - 1) * t1 - b * t2 - obs * a * log(1/b) - obs * log(gamma(a)))) 
} 

# initial guess for a = mean^2(x)/var(x) and b = mean(x)/var(x) 
a1 <- (mean(newdata$faminc))^2/var(newdata$faminc) 
b1 <- mean(newdata$faminc)/var(newdata$faminc) 

q <- optim(c(a1, b1), llf) 
q$par 
[1] 2.024353 3.019376 

Je dirais que nous sommes assez proches.

Avec vos données:

(est <- q$par) 
[1] 2.21333613 0.04243384 

theoretical <- data.table(true = rgamma(10000, est[1], est[2])) 
library(ggplot2) 
ggplot(newdata, aes(x = faminc)) + geom_density() + geom_density(data = theoretical, aes(x = true, colour = "red")) + theme(legend.position = "none") 

Pas génial, mais raisonnable pour 500 obs.

Réponse à OP Edit 2:

Vous devriez regarder de plus près les fonctions que vous utilisez, curve accepte un argument de fonction, non valeurs vectorielles:

gamma_density = function(x, a, b) ((b^a)/gamma(a)) * (x^(a - 1)) * exp(-b * x) 
hist(newdata$faminc, prob = TRUE, ylim = c(0, 0.015)) 
curve(gamma_density(x, a = q$par[1], b = q$par[2]), add=TRUE, col='blue')