In jtrecenti/tensorglm:

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

$f$ é a log-verossimilhança

Newton raphson:

Queremos fazer a derivada da log-verossimilhanca ($f$) igual a zero

$$ f'(\beta) = 0 $$

$$ \beta_{n+1} = \beta_n - \frac{f'(x_n)}{f''(x_n)} $$

utiliza teorema do valor médio (cálculo 1)
usa segunda derivada
converge mais rápido

Gradient descent:

$$ \beta_{n+1} = \beta_{n} - \alpha f'(\beta_n) $$

library(tensorflow)
library(tensorglm)

Normal

m_lm <- lm(carat ~ x + z, data = ggplot2::diamonds)
m_lm_tf <- glmtf(carat ~ x + z, data = ggplot2::diamonds, n = 1000, lr = 0.7)
m_lm_tf <- glmk(carat ~ x + z, data = ggplot2::diamonds, epochs = 100,
                family = gaussian, lr = 0.6)



coef(m_lm)
unlist(m_lm_tf)

Gamma

form <- mpg ~ disp
d <- mtcars

m_gamma <- glm(form, data = d, family = Gamma(link = 'log'))
m_gamma_tf <- glmtf(form, data = d, family = Gamma(link = 'log'))

coef(m_gamma)
unlist(m_gamma_tf)

library(magrittr)
  m_gamma_k <- glmk(form, data = d, family = Gamma(link = 'log'),
                  batch_size = 10, lr = 0.01)

m_gamma
m_gamma_k

Poisson

d <- tibble::tibble(
  x = runif(1e5, min = 0, max = 1), 
  y = rpois(1e5, exp(1.5 + 0.1 * x))
)

m_poisson <- glm(y ~ x, data = d, family = poisson)
m_poisson_tf <- glmtf(y ~ x, data = d, family = poisson,
                      lr = 0.05, n = 200)

m_poisson_k <- glmk(y ~ x, data = d, family = poisson,
                    batch_size = 100, epochs = 10,
                    lr = 0.01)


m_poisson
m_poisson_k

round(coef(m_poisson), 5)
unlist(m_poisson_tf)

Binomial

N <- 10000
x1 <- rnorm(N)
x2 <- rnorm(N)
eta <- 1 + 2 * x1 + 3 * x2
logit <- 1 / (1 + exp(-eta))
x <- cbind(1, x1, x2)
y <- rbinom(N, 1, logit)
d <- data.frame(y, x1, x2)

y <- rnorm(n, eta)

fit2 <- lm(y ~ x)
fit <- .lm.fit(x, y)


# qr
#

# coeficientes
# solve(t(x) %*% x) %*% t(x) %*% matrix(y, ncol = 1)

set.seed(1)
N <- 10000
X <- matrix(rnorm(3*N), ncol = 3)
y <- rnorm(N, apply(X, 1, sum))

coef(lm(y ~ X+0))
lm.fit(X, y)$coefficients
.lm.fit(X, y)$coefficients
qr.solve(X, y)
solve(qr(X), y)
solve(t(X) %*% X) %*% t(X) %*% y



library(tensorflow)

set.seed(1)
N <- 10000
X <- matrix(rnorm(3*N), ncol = 3)
y <- rnorm(N, apply(X, 1, sum))

lm_qr_tf <- function(X, y) {
  X_ <- tf$to_float(X)
  y_ <- tf$to_float(matrix(y, ncol = 1))
  QR <- tf$qr(X_)
  qy <- tf$matmul(QR$q, y_, transpose_a = TRUE)
  beta <- tf$matrix_solve(QR$r, qy)
  sess <- tf$Session()
  sess$run(beta)
}

lm_qr_tf(X, y)


m <- microbenchmark::microbenchmark(
  qr.solve(X, y),
  lm_qr_tf(X, y)
)




solve(qr.R(qr(X))) %*% t(qr.Q(qr(X))) %*% y


fit$coefficients

head(fit$qr)

all.equal(fit$residuals, as.numeric(resid(fit2)))




summary(resp2) <- glm(y ~ x1 + x2, data = d, family = binomial())
summary(resp2_2) <- glm2(y ~ x1 + x2, data = d, family = binomial())

resp3 <- glmtf(y ~ x1 + x2, data = d, family = binomial(),
               lr = 0.01, n = 10000, verbose = TRUE,
               scale = FALSE)

resp <- glmk(y ~ x1 + x2, data = d, family = binomial(), 
             lr = 2,
             epochs = 100)
print(resp)