Как получить BIC / AIC-график для выбора количества основных компонентов в Python или R

Question

Я хочу получить график, подобный этому, для выбора количества компонентов в PCA:

Однако я застрял, пытаясь вручную кодировать BIC / AIC. Есть ли пакеты в R или Python, которые могут помочь мне получить это? Любой пример кода очень помог бы.

Спасибо

Answer 1

Здесь ссылка на некоторый пример кода R, который вычисляет AIC и BIC, а также выбор переменных вперед / назад / пошагово.Вся заслуга принадлежит Джо Хардин.Для удобства я воспроизведу часть приведенного ниже кода, слегка отредактированного для форматирования:

> sat.data <- read.table("sat.csv", header=T, sep=",")
> attach(sat.data)
> sat.n <- nrow(sat.data) # be careful with missing values!!
> ltakers <- log(takers) # variable is quite right skewed

AIC и BIC в R

Метод 1:

> sat.lm0 <- lm(sat ~ 1)
> summary(sat.lm0)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 948.45 10.21 92.86 <2e-16 ***
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

Residual standard error: 71.5 on 48 degrees of freedom

> sat.sse0 <- sum(resid(sat.lm0) ^2)
> sat.n + sat.n*log(2*pi) + sat.n * log(sat.sse0 / sat.n) + 2 * (1+1)
[1] 560.4736
> AIC(sat.lm0, k=2)
[1] 560.4736
> sat.n + sat.n * log(2*pi) + sat.n*log(sat.sse0/sat.n) + log(sat.n)*(1+1)
[1] 564.2573
> AIC(sat.lm0, k=log(sat.n))
[1] 564.2573

Метод 2:

> sat.lm1 <- lm(sat ~ ltakers)
> summary(sat.lm1)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1112.408 12.386 89.81 <2e-16 ***
ltakers -59.175 4.167 -14.20 <2e-16 ***
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

Residual standard error: 31.41 on 47 degrees of freedom
Multiple R-squared: 0.811, Adjusted R-squared: 0.807
F-statistic: 201.7 on 1 and 47 DF, p-value: < 2.2e-16

> sat.sse1 <- sum(resid(sat.lm1) ^2)
> sat.n + sat.n*log(2*pi) + sat.n * log(sat.sse1 / sat.n) + 2 * (2+1)
[1] 480.832
> AIC(sat.lm1, k=2)
[1] 480.832
> sat.n + sat.n * log(2*pi) + sat.n*log(sat.sse1/sat.n) + log(sat.n) * (2+1)
[1] 486.5075
> AIC(sat.lm1, k=log(sat.n))
[1] 486.5075