我对R语言还比较新手。我有一组数据，其中包括家庭收入数据，我需要使用最大似然估计对这些数据拟合Gamma分布。特别指出，我们需要使用optim包，而不是fitdistr包。所以这是我的代码：

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]))}#初始猜测a = mean^2(x)/var(x) 和 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

我也尝试过直接在init向量中填入数值，并加入以下代码：

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

然后optim的代码看起来像这样：

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

这些尝试都没有成功。首先，当我在学校的电脑上运行代码时，形状和速率参数得到了非常大的数值，这是不可能的。现在，在家里尝试，我得到了以下结果：

> 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)> qfunction (save = "default", status = 0, runLast = TRUE) .Internal(quit(save, status, runLast))<bytecode: 0x000000000eaac960><environment: namespace:base>

q甚至没有被创建。除了当我加入上面的dlh.gamma部分时，但我再次得到了非常大的数值，并且没有收敛。

有谁知道哪里出了问题/该怎么做？

编辑：

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

faminc变量的单位是千元

编辑2：

好的，代码是正确的，但现在我尝试使用以下代码在直方图上拟合分布：

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

它只是在x轴上产生了一条平坦的蓝色线…

回答：