Точность алгоритма для регрессии Пуассона (glm, R)

Question

Я хочу сделать регрессию Пуассона, но мне нужна моя функция регрессии, чтобы работать быстрее, чем glm и, по крайней мере, с такой точностью. Рассмотрим следующий эксперимент:

## Here is some "data": 
da = data.frame(matrix(c(0,1,212,1,0,200,1,1,27), nrow = 3, byrow = TRUE)) 
names(da) = c("c1", "c2", "c") 

## I want to do a Poisson regression of c on c1 and c2 and an intercept. 

## Here is my function that uses optim for Poisson regression with the data da to find the intercept term: 
zglm2 = function(precision = 1){ #predictors = best.terms, data = ddat, normalized = normalized 
    # The design matrix 
    M = as.matrix(cbind(rep(1, nrow(da)), da[,1:2])) 
    # Initialize beta, the coefficients 
    beta = rep(0, 3) 
    # State the log-likelihood (up to a constant) for the data da and parameter beta: 
    neg.pois.log.like.prop = function(beta){ 
    log.lambda = M%*%beta # log-expected cell counts under poisson model 
    return(-sum(-exp(log.lambda) + da$c*log.lambda))} 
    # State the gradient of the log-likelihood: 
    grad.fun = function(beta){a = exp(M%*%beta)-da$c; return(t(a)%*%M)} 
    # Estimate the MLE 
    beta = optim(beta, neg.pois.log.like.prop, method = "BFGS", gr = grad.fun, control = list(reltol = precision*sqrt(.Machine$double.eps)))$par 
    return(beta[1])} 

## Here are two ways of estimating the intercept term: 
# Method 1 
zglm2(precision = 1) 
# Method 2 
as.numeric(glm(c ~ 1+c1+c2, data = da, family = poisson)$coefficients[1]) 

Моя функция, zglm2 использует подпрограмму R в optim найти максимальное решение правдоподобия для регрессионной задачи Пуассона (для этого особого случая). zglm2 принимает аргумент precision; значения этого аргумента, которые меньше 1, делают optim превышать критерии завершения по умолчанию для достижения большей точности.

К сожалению, результаты метода 1 и метода 2 слишком разные (для моих целей); 7.358 против 7.359. Предоставление меньшего значения, например 0,01, для аргумента precision приводит к разумному соглашению двух методов, что приводит меня к подозрению, что функция Rочень точная.

Итак, вот мой вопрос: что определяет уровень точности в результате glm? Возможно, в качестве подзапроса, какой алгоритм использует glm, чтобы найти максимальную вероятность (я врыл в исходный код, но это нелегко для меня).

Answer 1

Я с трудом веря, что «вы закопали в код» очень глубоко, так как есть аа параметр «контроль» и функция, которая доставляет argumetn к glm:

?glm 
# control 
#  a list of parameters for controlling the fitting process. 
#  For glm.fit this is passed to glm.control.