Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Como resolver em python o problema de uma regressão linear simples segmentada?

+1 voto
440 visitas
perguntada Mai 20, 2016 em Ciência da Computação por CaioP (21 pontos)  
Compartilhe

1 Resposta

0 votos
respondida Mai 20, 2016 por CaioP (21 pontos)  
editado Mai 24, 2016 por CaioP

A solução do problema passa pela estrutura usual dos problemas de programação dinâmica, uma formulação recursiva da função valor. Tal formulação permite que se resolva o problema combinando soluções de subproblemas.

No caso particular desse problema, a ideia é testar as diferentes partições/quebras da amostra, tendo como valor a ser minimizado o erro quadrático, e levando-se em conta um custo arbitrário C (input) acrescido a cada quebra da amostra, conforme a equação abaixo:

\(
V(j)=min_{1\leq i\leq j}[SSR(i,j) + C+ V(i-1)]
\)

Em que SSR(i,j) representa a soma dos erros quadráticos da regressão simples aplicada na partição da amostra que vai do elemento i ao elemento j da amostra. Dessa forma, i* é o primeiro "ponto de quebra" ótimo, a partir de uma amostra de tamanho j. Naturalmente, o caso base dessa formulação é V(0)=0.

A implementação abaixo, em python, segue 1 e a resposta que se obtém, armazenada no dicionário arg_min deve ser interpretada sequencialmente. Em outras palavras, para cada tamanho da amostra (key), mostra-se o tamanho seguinte, isto é, após a retirada do primeiro "ponto de quebra" ótimo para aquele tamanho, e assim sucessivamente...

from __future__ import division

class OLS:
    def __init__(self,a,b,SE): #classe OLS, armazena coeficientes e erro quadratico
        self.a=a 
        self.b=b 
        self.SE=SE

def results (x,y): #roda um OLS simples e retorna um objeto da classe OLS
    n = len(x)
    if len(x)>1:
        soma_x=0
        soma_y=0
        soma_xy=0
        soma_xx=0
        for i in range(0,n):
            soma_x=soma_x + x[i]
            soma_y=soma_y + y[i]
            soma_xy=soma_xy + x[i]*y[i]
            soma_xx=soma_xx + x[i]*x[i]
        a = (n*soma_xy-soma_x*soma_y)/(n*soma_xx-soma_x*soma_x)
        b = (soma_y-a*soma_x)/n
        SE = 0
        for i in range(0,n):
            SE = SE + (y[i]-a*x[i]-b)*(y[i]-a*x[i]-b)
        results = OLS(a,b,SE)
        return results
    else:
        a = 0
        b = 0
        SE = 1000000
        results = OLS(a,b,SE)
        return results

def get_results(x,y): #armazena o erro quadratico de todas as particoes da amostra
    n = len(x)
    matrix_results = [[1000000 for l in range(0,n+1)] for k in range(0,n+1)]
    for j in range(0,n+1):
        for i in range(0,j+1):
            aux_x = x[i:j]
            aux_y = y[i:j]
            aux_r = results(aux_x,aux_y)
            matrix_results[i][j]=aux_r.SE
    return matrix_results

arg_min = {}

def opt(x,y,C,length): #funcao valor recursiva
    if length <= 0:
        return 0
    else:
        e = get_results(x,y)
        mini = 1000000
        for i in range(1,length+1):
            aux = e[i-1][length] + C + opt(x,y,C,i-1)
            if aux < mini:
                mini = aux
                global arg_min
                arg_min[length]=i-1
        return mini

x = [1,2,3,4,5,6]
y = [2,8,13,32,41,67]

print get_results(x,y)
a = []
for i in range(0,len(x)+1):
    print i
    a.append(opt(x,y,0,i))
print arg_min #arg_min indica o argumento otimo para a "primeira etapa" daquele valor
comentou Mai 20, 2016 por danielcajueiro (5,326 pontos)  
Caio, acho que seria legal se você pudesse documentar melhor sua solução. No momento, a documentação está disponível no blog que você citou. Entretanto, não temos como saber se essa documentação estará disponível lá no próximo mês ou próximo ano.
...