python – 对大型数据集运行逻辑回归的有效方法是什么(2亿个变量为2亿个)?

我目前正在尝试运行逻辑回归模型.我的数据有两个变量,一个响应变量和一个预测变量.问题是我有2亿观察.我正在尝试运行逻辑回归模型,但即使在亚马逊上的EC2实例的帮助下,在R / Stata / MATLAB中也是如此.我认为问题在于如何在语言本身中定义逻辑回归函数.还有另一种快速运行逻辑回归的方法吗?目前我遇到的问题是我的数据很快就填满了它正在使用的空间.我甚至尝试使用高达30 GB的RAM无济于事.任何解决方案都将受到欢迎.

解决方法:

如果您的主要问题是能够在给定计算机内存限制的情况下估计logit模型,而不是估计的快速性,则可以利用最大似然估计的可加性并为ml编写自定义程序.Logit模型只是一个使用逻辑分布估计最大似然.您只有一个自变量的事实简化了这个问题.我已经模拟了下面的问题.您应该从以下代码块中创建两个do文件.

如果您在整个数据集中没有加载任何问题 – 您不应该使用,我的模拟仅使用〜2演出的RAM使用2亿个障碍和2个变量,尽管里程可能会有所不同 – 第一步是将数据集分解为易于管理的作品.例如:

depvar =你的因变量(0或1s)
indepvar =您的自变量(某些数值数据类型)

cd "/path/to/largelogit"

clear all
set more off

set obs 200000000

// We have two variables, and independent variable and a dependent variable.
gen indepvar = 10*runiform()
gen depvar = .

// As indpevar increases, the probability of depvar being 1 also increases.
replace depvar = 1 if indepvar > ( 5 + rnormal(0,2) )
replace depvar = 0 if depvar == .

save full, replace
clear all

// Need to split the dataset into managable pieces

local max_opp = 20000000    // maximum observations per piece

local obs_num = `max_opp'

local i = 1
while `obs_num' == `max_opp' {

    clear

    local h = `i' - 1

    local obs_beg = (`h' * `max_opp') + 1
    local obs_end = (`i' * `max_opp')

    capture noisily use in `obs_beg'/`obs_end' using full

    if _rc == 198 {
        capture noisily use in `obs_beg'/l using full
    }
    if _rc == 198 { 
        continue,break
    }

    save piece_`i', replace

    sum
    local obs_num = `r(N)'

    local i = `i' + 1

}

从这里开始,尽量减少内存使用量,关闭Stata并重新打开它.当您创建如此大的数据集时,即使清除数据集,Stata也会为开销等分配一些内存.你可以在保存完全后和清除之后键入内存以查看我的意思.

接下来,您必须定义自己的自定义ml程序,该程序将在程序中一次一个地输入每个部分,计算每个部分的每个观察的对数似然并将它们加起来并将它们全部加在一起.您需要使用d0 ml方法而不是lf方法,因为使用lf的优化例程需要将所有数据加载到Stata中.

clear all
set more off

cd "/path/to/largelogit"

// This local stores the names of all the pieces 
local p : dir "/path/to/largelogit" files "piece*.dta"

local i = 1
foreach j of local p {    // Loop through all the names to count the pieces

    global pieces = `i'    // This is important for the program
    local i = `i' + 1

}

// Generate our custom MLE logit progam. This is using the d0 ml method 

program define llogit_d0

    args todo b lnf 

    tempvar y xb llike tot_llike it_llike

quietly {

    forvalues i=1/$pieces {

        capture drop _merge
        capture drop depvar indepvar
        capture drop `y'
        capture drop `xb'
        capture drop `llike' 
        capture scalar drop `it_llike'

        merge 1:1 _n using piece_`i'

        generate int `y' = depvar

        generate double `xb' = (indepvar * `b'[1,1]) + `b'[1,2]    // The linear combination of the coefficients and independent variable and the constant

        generate double `llike' = .

        replace `llike' = ln(invlogit( `xb')) if `y'==1    // the log of the probability should the dependent variable be 1
        replace `llike' = ln(1-invlogit(`xb')) if `y'==0   // the log of the probability should the dependent variable be 0

        sum `llike' 
        scalar `it_llike' = `r(sum)'    // The sum of the logged probabilities for this iteration

        if `i' == 1     scalar `tot_llike' = `it_llike'    // Total log likelihood for first iteration
        else            scalar `tot_llike' = `tot_llike' + `it_llike' // Total log likelihood is the sum of all the iterated log likelihoods `it_llike'

    }

    scalar `lnf' = `tot_llike'   // The total log likelihood which must be returned to ml

}

end

//This should work

use piece_1, clear

ml model d0 llogit_d0 (beta : depvar = indepvar )
ml search
ml maximize

我刚刚运行了上面两段代码并收到了以下输出:

这种方法的优点和缺点:
优点:

– “max_opp”大小越小,内存使用率越低.如上所述,我从未使用过模拟器超过~1 gig.

– 您最终得到无偏估计,整个数据集的估计的完整对数似然,正确的标准误差 – 基本上所有对推理都很重要.

缺点:

 
– 你在内存中保存的内存,你必须牺牲CPU时间.我使用带有i5处理器的Stata SE(一个核心)在我的个人笔记本电脑上运行它,它让我一夜之间.

– Wald Chi2统计数据是错误的,但我相信你可以根据上面提到的正确数据来计算它

– 你没有得到像logit那样的Psudo R2.

要测试系数是否真的与标准logit相同,请将obs设置为相对较小的值,100000,并将max_opp设置为1000.运行我的代码,查看输出,运行logit depvar indepvar,查看输出,它们与我在上面的“缺点”中提到的相同.将obs设置为与max_opp相同将更正Wald Chi2统计数据.

上一篇:Linux上的Stata 13:如何使用export delimited处理文件路径中的空格?


下一篇:Stata 中Mata的st_view函数 - 小心使用!