Ajustando duas curvas com regressão linear / não linear

Preciso ajustar duas curvas (que devem pertencer a funções cúbicas) em um conjunto de pontos com o JuMP.

Eu já fiz o ajuste de uma curva, mas estou lutando para ajustar 2 curvas no mesmo conjunto de dados.

Eu pensei que se eu puder distribuir pontos em curvas - então, se cada ponto puder ser usado apenas uma vez - eu posso fazê-lo como abaixo, mas não funcionou. (Eu sei que posso usar coisas muito mais complicadas, quero simplificar.)

Esta é uma parte do meu código atual:

# cubicFunc is a two dimensional array which accepts cubicFunc[x,degree]

@variable(m, mult1[1:4]) // 0:3 because it's cubic
@variable(m, mult2[1:4]) // 0:3 because it's cubic

@variable(m, 0 <= includeIn1[1:numOfPoints] <= 1, Int)
@variable(m, 0 <= includeIn2[1:numOfPoints] <= 1, Int)

# some kind of hack to force one of them to 0 and other one to 1
@constraint(m, loop[i in 1:numOfPoints], includeIn1[i] + includeIn2[i] == 1)

@objective(m, Min, sum( (yPoints - cubicFunc*mult1).*includeIn1 .^2 ) + sum( (yPoints - cubicFunc*mult2).*includeIn2 .^2 ))

Mas dá vários erros, dependendo do que estou tentando; *includeIn1e, .*includeIn1não funciona, tentei fazê-lo via, @NLobjectivemas isso me deu gritos ~ 50 linhas de erros etc.

A minha ideia é realista? Posso entrar no código?

Qualquer ajuda será muito apreciada. Muito obrigado.

Você pode anotar o problema, por exemplo:

using JuMP, Ipopt

m = Model(with_optimizer(Ipopt.Optimizer))

@variable(m, mult1[1:4])
@variable(m, mult2[1:4])
@variable(m, 0 <= includeIn1[1:numOfPoints] <= 1)
@variable(m, 0 <= includeIn2[1:numOfPoints] <= 1)

@NLconstraint(m, loop[i in 1:numOfPoints], includeIn1[i] + includeIn2[i] == 1)

@NLobjective(m, Min, sum(includeIn1[i] * (yPoints[i] - sum(cubicFunc[i,j]*mult1[j] for j in 1:4)) ^2 for i in 1:numOfPoints) +
                     sum(includeIn2[i] * (yPoints[i] - sum(cubicFunc[i,j]*mult2[j] for j in 1:4)) ^2 for i in 1:numOfPoints))

optimize!(m)

Dadas as restrições includeIn1e includeIn2serão 1ou estão 0ótimas (se não forem, isso significa que não importa a qual grupo você atribui o ponto), portanto, não precisamos construí-los para serem binários. Também uso o solucionador não linear, pois o problema não parece ser possível reformular como tarefa de otimização linear ou quadrática.

No entanto, dou o código acima apenas como exemplo, como você pode anotá-lo. A tarefa que você formulou não possui um mínimo local exclusivo (que é global então), mas vários mínimos locais. Portanto, o uso de solucionadores convexos não lineares padrão suportados pelo JuMP encontrará apenas um ótimo local (não necessariamente global). Para procurar ótimas globais, você precisa mudar para solucionadores globais como, por exemplo, https://github.com/robertfeldt/BlackBoxOptim.jl .