2016-10-03 20 views
0

所以我想最大限度地为伽玛泊松似然函数和我编程式的R为以下几点:二元函数的最大化 - R的代码

lik<- function(x,t,a,b){  
    for(i in 1:n){ 
     like[i] = 
     log(gamma(a + x[i]))-log(gamma(a)) 
      -log(gamma(1+x[i] + x[i]*log(t[i]/b)-(a+x[i])*log(1+t[i]/b) 
     } 
     return(sum(like)) 
} 

其中X和T是数据,我有n个数据行。

我需要a和b来同时解决。 R中是否存在内置函数?还是我需要硬编码算法来解决方程组? [我宁可不]我知道optimize()解决了1个变量,fminbnd()也解决了这个问题。我试图复制FindMaximum()在mathematica中的行为。在一个完美的世界,我想工作代码东西这样的:

optimize(f=lik, a>0, b>0, x=x, t=t, maximum=TRUE, iteration=5000) 
$maximum 
    a 150 
    b 6 

感谢。

+0

是的,这可能是我想要的,我刚刚过了一个小时google搜索,并在此方案中R. –

+0

没有一(1)看'?optim',要记住,*减少*默认情况下(在帮助文件中看到'fnscale'的描述); (2)'dnbinom'可能有帮助; (3)所以可能'MASS :: fitdistr' –

回答

1

optim' s第一个参数可以是参数的向量。所以,你可以尝试这样的事:

lik <- function(p=c(1,1), x, t){ 
    # In the body of the function replace a by p[1] and b by p[2] 
} 

optim(c(1,1), lik, method = c("L-BFGS-B"), x=x, t=t, control=list(fnscale=-1)) 
0

这样结束了工作出来的解决方案是:

attempt2d <- optim(
    par = c(sumx/sumt, 1), fn = lik, data = data11, 
    method = "L-BFGS-B", control = list(fnscale = -1, trace=TRUE), 
    lower=0.1, upper = 170 
    ) 

但是我的参数跑出去170,实质上意味着我的伽马参数是天道酬勤。因为伽马()比较快地命中无限。在mathematica中,解为a = 169和b = 16505,而R在170处最大值远不及170.在某些情况下,已知的解决方案超出了170个解决方案。