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
