Como escolher valores iniciais para o ajuste de mínimos quadrados não lineares

12

A pergunta acima diz tudo. Basicamente, minha pergunta é para uma função de ajuste genérico (poderia ser arbitrariamente complicada) que não será linear nos parâmetros que estou tentando estimar, como escolher os valores iniciais para inicializar o ajuste? Estou tentando fazer mínimos quadrados não lineares. Existe alguma estratégia ou método? Isso foi estudado? Alguma referência? Algo além de adivinhação ad hoc? Especificamente, no momento, uma das formas de ajuste com as quais estou trabalhando é uma forma gaussiana e linear com cinco parâmetros que estou tentando estimar, como

y=Ae(xBC)2+Dx+E

onde (dados de abscissa) e y = log 10x=log10y=log10 (dados ordenados), o que significa que no espaço de log-log meus dados se parecem com uma linha reta mais um solavanco que eu estou aproximando por um gaussiano. Não tenho teoria, nada para me guiar sobre como inicializar o ajuste não-linear, exceto, talvez, gráficos e globos oculares como a inclinação da linha e qual é o centro / largura do solavanco. Mas eu tenho mais de cem desses ajustes para fazê-lo, em vez de fazer gráficos e adivinhar, eu preferiria alguma abordagem que possa ser automatizada.

Não consigo encontrar nenhuma referência, na biblioteca ou online. A única coisa em que consigo pensar é escolher aleatoriamente os valores iniciais. O MATLAB oferece para escolher valores aleatoriamente entre [0,1] distribuídos uniformemente. Portanto, com cada conjunto de dados, eu executo o ajuste inicializado aleatoriamente mil vezes e, em seguida, escolho aquele com o r 2 mais altor2 ? Alguma outra (melhor) ideia?


Adendo # 1

Primeiro, aqui estão algumas representações visuais dos conjuntos de dados apenas para mostrar a vocês que tipo de dados eu estou falando. Estou postando os dados em sua forma original sem qualquer tipo de transformação e, em seguida, sua representação visual no espaço de log-log, uma vez que esclarece alguns dos recursos dos dados enquanto distorce outros. Estou postando uma amostra de dados bons e ruins.

bons dados log-log de bons dados dados incorretos dados incorretos

Cada um dos seis painéis de cada figura mostra quatro conjuntos de dados plotados em vermelho, verde, azul e ciano e cada conjunto de dados possui exatamente 20 pontos de dados. Estou tentando encaixar cada um deles com uma linha reta mais uma gaussiana por causa dos solavancos vistos nos dados.

A primeira figura é alguns dos bons dados. A segunda figura é o gráfico log-log dos mesmos dados válidos da figura um. A terceira figura é alguns dos dados incorretos. A quarta figura é o gráfico log-log da figura três. Há muito mais dados, esses são apenas dois subconjuntos. A maioria dos dados (cerca de 3/4) é boa, semelhante aos bons dados que mostrei aqui.

Agora, alguns comentários, por favor, tenham paciência comigo, pois isso pode demorar, mas acho que todos esses detalhes são necessários. Vou tentar ser o mais conciso possível.

Originalmente, eu esperava uma lei de energia simples (ou seja, linha reta no espaço de log-log). Quando plotei tudo no espaço de log-log, vi um aumento inesperado em torno de 4,8 mHz. A colisão foi minuciosamente investigada e foi descoberta em outros trabalhos também, por isso não é que erramos. Está fisicamente lá e outros trabalhos publicados mencionam isso também. Então, acabei de adicionar um termo gaussiano à minha forma linear. Observe que esse ajuste deveria ser feito no espaço de log-log (daí minhas duas perguntas, incluindo esta).

Agora, depois de ler a resposta de Stumpy Joe Pete para outra pergunta minha (que não está relacionada a esses dados) e de ler isso e isso e fazer referência a elas (material de Clauset), percebo que não devo caber no log-log espaço. Então agora eu quero fazer tudo no espaço pré-transformado.

Pergunta 1: Observando os bons dados, ainda acho que um linear mais um gaussiano no espaço pré-transformado ainda é uma boa forma. Gostaria muito de ouvir de outras pessoas que têm mais experiência em dados o que pensam. Gaussian + linear é razoável? Devo fazer apenas um gaussiano? Ou uma forma completamente diferente?

Perguntas 2: Qualquer que seja a resposta à pergunta 1, eu ainda precisaria (provavelmente) de mínimos quadrados não lineares, por isso ainda precisaria de ajuda com a inicialização.

Nos dados em que vemos dois conjuntos, preferimos muito capturar a primeira colisão em torno de 4-5 mHz. Portanto, não quero adicionar mais termos gaussianos, e nosso termo gaussiano deve estar centrado no primeiro solavanco, que quase sempre é o maior solavanco. Queremos "mais precisão" entre 0.8mHz e cerca de 5mHz. Não nos importamos muito com as frequências mais altas, mas também não queremos ignorá-las completamente. Então, talvez algum tipo de pesagem? Ou B pode ser inicializado sempre em torno de 4.8mHz?

Os dados de abscissa são a frequência em unidades de milihertz, denotados por . A ordenada dados é um coeficiente que são computação, denotar que por L . Portanto, nenhuma transformação de log e o formulário éfL

L=Ae(fBC)2+Df+E.
  • é frequência, é sempre positivo.f
  • é um coeficiente positivo. Então, estamos trabalhando no primeiro quadrante.L
  • , a amplitude sempre deve ser positiva também, eu acho, porque estamos apenas lidando com solavancos. Quando olho para os dados, sempre vejo picos e sem vales. Parece que em todos os dados existem vários inchaços em frequências mais altas. O primeiro solavanco é sempre muito maior que os outros. Em dados bons, os solavancos secundários são muito fracos, mas em dados ruins (painéis 2 e 5, por exemplo), os solavancos secundários são fortes. Então, na verdade, nãotemos um vale, mas sim dois solavancos. Significando que a amplitude A > 0 . E como nos preocupamos principalmente com o primeiro pico, mais uma razão para A ser positivo.AA>0A
  • B
  • é a largura do solavanco. Eu imagino que é simétrico em torno de zero significadoCCC
  • D
  • ELELf=0

Ae(B/C)2+E.

EEf=0

L

Perguntas 3: O que vocês acham extrapolar dessa maneira neste caso? Quaisquer prós / contras? Alguma outra idéia para extrapolação? Novamente, nos preocupamos apenas com as frequências mais baixas, extrapolando entre 0 e 1mHz ... às vezes, frequências muito pequenas, próximas de zero. Eu sei que este post já está empacotado. Fiz essa pergunta aqui porque as respostas podem estar relacionadas, mas se vocês preferirem, posso separar essa pergunta e fazer outra mais tarde.

Por fim, aqui estão dois conjuntos de dados de amostra, mediante solicitação.

0.813010000000000   0.091178000000000   0.012728000000000
1.626000000000000   0.103120000000000   0.019204000000000
2.439000000000000   0.114060000000000   0.063494000000000
3.252000000000000   0.123130000000000   0.071107000000000
4.065000000000000   0.128540000000000   0.073293000000000
4.878000000000000   0.137040000000000   0.074329000000000
5.691100000000000   0.124660000000000   0.071992000000000
6.504099999999999   0.104480000000000   0.071463000000000
7.317100000000000   0.088040000000000   0.070336000000000
8.130099999999999   0.080532000000000   0.036453000000000
8.943100000000001   0.070902000000000   0.024649000000000
9.756100000000000   0.061444000000000   0.024397000000000
10.569000000000001   0.056583000000000   0.025222000000000
11.382000000000000   0.052836000000000   0.024576000000000
12.194999999999999   0.048727000000000   0.026598000000000
13.008000000000001   0.045870000000000   0.029321000000000
13.821000000000000   0.041454000000000   0.067300000000000
14.633999999999999   0.039596000000000   0.081800000000000
15.447000000000001   0.038365000000000   0.076443000000000
16.260000000000002   0.036425000000000   0.075912000000000

A primeira coluna são as frequências em mHz, idênticas em cada conjunto de dados. A segunda coluna é um bom conjunto de dados (dados bons, figura um e dois, painel 5, marcador vermelho) e a terceira coluna é um conjunto de dados incorretos (dados ruins, figura três e quatro, painel 5, marcador vermelho).

Espero que isso seja suficiente para estimular uma discussão mais esclarecida. Obrigado a todos.

Ponto fixo
fonte
+1 para obter informações adicionais, mas agora isso parece uma nova pergunta. Aliás, se você deseja excluir a anterior agora, acho que tudo bem, parece que você já cobriu as informações adicionais que ela possuía.
Glen_b -Reinstar Monica
@Glen_b Por que é isso? Por que parece uma nova pergunta? Quanto à velha pergunta, todos nós merecemos pontos ;-D e a antiga tem dois votos positivos, alguma maneira de fundi-la com isso para que eu possa manter esses dois votos também?
Fixed Point
Bem, para iniciantes, agora você está perguntando sobre o que deve caber, em vez de especificar o que caber, como antes. Existem várias outras diferenças, algumas das quais eu consideraria bastante significativas. Vou tentar mudar minha resposta, mas acho que essa poderia ser a pergunta e resposta originais e as novas partes em que você está perguntando outras coisas podem ser novas. Vou deixar isso para o seu julgamento no momento.
Glen_b -Reinstala Monica
@Glen_b Justo, fiz perguntas extras. Portanto, as perguntas ainda são: tenho alguns dados que quero ajustar usando a forma linear + gaussiana. Posso fazer melhor do que a inicialização aleatória?
Fixed Point
Acho que minha resposta atual mostra que - em pelo menos algumas circunstâncias - você pode fazer melhor, e o @whuber sugere algo ainda mais simples do que o meu processo. Eu poderia voltar e ver como o que eu tenho desempenha nos seus dados, mas, mesmo como estão agora, dá uma idéia de como configurar esses pontos de partida.
Glen_b -Reinstala Monica

Respostas:

10

Se houvesse uma estratégia que fosse boa e geral - uma que sempre funcionasse - ela já seria implementada em todos os programas de mínimos quadrados não lineares e os valores iniciais seriam um problema.

Para muitos problemas específicos ou famílias de problemas, existem algumas abordagens muito boas para iniciar valores; alguns pacotes vêm com bons cálculos de valor inicial para modelos não-lineares específicos ou com abordagens mais gerais que geralmente funcionam, mas podem precisar de ajuda com funções mais específicas ou entrada direta de valores iniciais.

Explorar o espaço é necessário em algumas situações, mas acho que sua situação provavelmente será tal que estratégias mais específicas provavelmente valerão a pena - mas projetar uma boa exige bastante conhecimento do domínio que provavelmente não possuiremos.

Para o seu problema específico, é provável que possam ser projetadas boas estratégias, mas não é um processo trivial; quanto mais informações você tiver sobre o tamanho e a extensão provável do pico (os valores típicos dos parâmetros e típicosx dariam alguma idéia), mais poderá ser feito para projetar uma boa opção de valor inicial.

yx

UMA sempre positivo? etc.

Alguns dados de amostra ajudariam - casos típicos e casos difíceis, se você puder.


Edit: Aqui está um exemplo de como você pode se sair bem se o problema não for muito barulhento:

Aqui estão alguns dados gerados a partir do seu modelo (os valores da população são A = 1,9947, B = 10, C = 2,828, D = 0,09, E = 5):

dados nls

Os valores iniciais que eu pude estimar são
(Como = 1,658, Bs = 10,001, Cs = 3,053, Ds = 0,0881, Es = 5,026)

O ajuste desse modelo inicial é assim:

nlstart

Os passos foram:

  1. Ajuste uma regressão de Theil para obter uma estimativa aproximada de D e E
  2. Subtraia o ajuste da regressão de Theil
  3. Use LOESS para ajustar uma curva suave
  4. Encontre o pico para obter uma estimativa aproximada de A e o valor x correspondente ao pico para obter uma estimativa aproximada de B
  5. Pegue os ajustes LOESS cujos valores y são> 60% da estimativa de A como observações e ajuste um quadrático
  6. Use o quadrático para atualizar a estimativa de B e estimar C
  7. A partir dos dados originais, subtraia a estimativa do valor gaussiano
  8. Ajuste outra regressão de Theil aos dados ajustados para atualizar a estimativa de D e E

Nesse caso, os valores são muito adequados para iniciar um ajuste não linear.

Eu escrevi isso como Rcódigo, mas o mesmo poderia ser feito no MATLAB.

Eu acho que coisas melhores do que isso são possíveis.

Se os dados forem muito barulhentos, isso não funcionará bem.


Edit2: Este é o código que usei no R, se alguém estiver interessado:

gausslin.start <- function(x,y) {

  theilreg <- function(x,y){
    yy <- outer(y, y, "-")
    xx <- outer(x, x, "-")
    z  <- yy / xx
    slope     <- median(z[lower.tri(z)])
    intercept <- median(y - slope * x)
    cbind(intercept=intercept,slope=slope)
  }

  tr <- theilreg(x,y1)
  abline(tr,col=4)
  Ds = tr[2]
  Es = tr[1]
  yf  <- y1-Ds*x-Es
  yfl <- loess(yf~x,span=.5)

  # assumes there are enough points that the maximum there is 'close enough' to 
  #  the true maximum

  yflf   <- yfl$fitted    
  locmax <- yflf==max(yflf)
  Bs     <- x[locmax]
  As     <- yflf[locmax]

  qs     <- yflf>.6*As
  ys     <- yfl$fitted[qs]
  xs     <- x[qs]-Bs
  lf     <- lm(ys~xs+I(xs^2))
  bets   <- lf$coefficients
  Bso    <- Bs
  Bs     <-  Bso-bets[2]/bets[3]/2
  Cs     <- sqrt(-1/bets[3])
  ystart <- As*exp(-((x-Bs)/Cs)^2)+Ds*x+Es

  y1a <- y1-As*exp(-((x-Bs)/Cs)^2)
  tr  <- theilreg(x,y1a)
  Ds  <- tr[2]
  Es  <- tr[1]
  res <- data.frame(As=As, Bs=Bs, Cs=Cs, Ds=Ds, Es=Es)
  res
}

.

# population parameters: A = 1.9947 , B = 10, C = 2.828, D = 0.09, E = 5
# generate some data
set.seed(seed=3424921)
x  <- runif(50,1,30)
y  <- dnorm(x,10,2)*10+rnorm(50,0,.2)
y1 <- y+5+x*.09 # This is the data
xo <- order(x)

starts <- gausslin.start(x,y1)
ystart <- with(starts, As*exp(-((x-Bs)/Cs)^2)+Ds*x+Es)
plot(x,y1)
lines(x[xo],ystart[xo],col=2)
Glen_b -Reinstate Monica
fonte
3
+1. Repetir o ajuste mil vezes e escolher o melhor (se bem entendi) parece uma idéia estranha: mínimos quadrados não lineares devem convergir se o modelo for razoável para os dados e houver bons valores iniciais. Naturalmente, o segundo é o que você está perguntando. Mas parece pessimista sugerir que você pode precisar escolher diferentes valores iniciais para cada ajuste.
Nick Cox
1
@NickCox Tudo se resume à variedade de problemas encontrados - se bem me lembro de postagens anteriores, o OP recebe um grande número desses problemas, mas não me lembro de ter visto detalhes suficientes para fazer boas sugestões antes, embora tenha investido um pouco de tempo brincando com possíveis abordagens (que não renderam nada definitivo o suficiente para postar). O OP provavelmente possui o tipo de conhecimento de domínio que pode gerar bons valores iniciais que resolvem seus problemas quase sempre.
Glen_b -Reinstar Monica
1
Bastante. Eu perdi o post anterior no stats.stackexchange.com/questions/61724/...
Nick Cox
3
|A|BA>0CA1/4A>0A<0
2
BB
6

Existe uma abordagem geral para ajustar esse tipo de modelo não linear. Isso envolve reparametrizar os parâmetros lineares com valores da variável dependente, digamos o primeiro, o último valor da frequência e um bom ponto no meio, digamos o sexto ponto. você pode manter esses parâmetros fixos e resolver o parâmetro não linear na primeira fase da minimização e, em seguida, minimizar os 5 parâmetros gerais.

Schnute e eu descobrimos isso por volta de 1982 ao ajustar modelos de crescimento para peixes.

http://www.nrcresearchpress.com/doi/abs/10.1139/f80-172

No entanto, não é necessário ler este artigo. Devido ao fato de os parâmetros serem lineares, é simplesmente necessário configurar e resolver um sistema linear de equações 3x3 para usar a parametrização estável do modelo.

M

M=(exp(-((x(1)-B)/C)2)x(1)1exp(-((x(6)-B)/C)2)x(6)1exp(-((x(n)-B)/C)2)x(n)1)
Onde n=20nesse caso. 1 O código está escrito no AD Model Builder, que agora é de código aberto. Faz diferenciação automática e suporta muitas coisas que facilitam a otimização não linear, como fases em que diferentes parâmetros são mantidos fixos.
DATA_SECTION
  init_int n
  int mid
 !! mid=6;
  init_matrix data(1,n,1,3)
  vector x(1,n)
  vector y(1,n)
 !! x=column(data,1);
 !! y=column(data,3);   //use column 3
PARAMETER_SECTION
  init_number L1(3)     //(3) means estimate in phase 3
  init_number Lmid(3)
  init_number Ln(3)

  vector L(1,3)
  init_number log_B       // estimate in phase 1
  init_number log_C(2)    // estimate in phase 2 
  matrix M(1,3,1,3);
  objective_function_value f
  sdreport_vector P(1,3)
  sdreport_number B
  sdreport_number C
  vector pred(1,n);
PROCEDURE_SECTION
  L(1)=L1;
  L(2)=Lmid;
  L(3)=Ln;
  B=exp(log_B);
  C=exp(log_C);
  M(1,1)=exp(-square((x(1)-B)/C));
  M(1,2)=x(1);
  M(1,3)=1;
  M(2,1)=exp(-square((x(mid)-B)/C));
  M(2,2)=x(mid);
  M(2,3)=1;
  M(3,1)=exp(-square((x(n)-B)/C));
  M(3,2)=x(n);
  M(3,3)=1;

  P=solve(M,L);  // solve for standard parameters 
                 // P is vector corresponding to A,D,E

  pred=P(1)*exp(-square((x-B)/C))+P(2)*x+P(3);
  if (current_phase()<4)
    f+=norm2(y-pred);
  else
    f+=0.5*n*log(norm2(y-pred))  //concentrated likelihood

O modelo é adequado em três fases. Apenas na fase 1BÉ estimado. Aqui, o mais importante é queC é grande o suficiente para que o modelo possa "ver" para onde mover Bpara. Na fase 2, ambosB e C são estimados. Finalmente, na fase 3, todos os parâmetros são estimados. No gráfico, a linha verde é o ajuste após a fase 1 e a linha azul é o ajuste final.

insira a descrição da imagem aqui

Para o seu caso com dados incorretos, ele se encaixa facilmente e as estimativas (usuais) de parâmetros são:

         estimate    std dev
A      2.0053e-01 5.8723e-02
D      1.6537e-02 4.7684e-03
E     -1.8197e-01 7.3355e-02
B      3.0609e+00 5.0197e-01
C      5.6154e+00 9.4564e-01]
Dave Fournier
fonte
Dave, isso é interessante, mas levanta algumas questões. Exatamente o que você quer dizer com "esse tipo de modelo não linear"? A pergunta começa com a referência a uma "função de ajuste genérico", mas sua descrição se refere apenas a "cinco parâmetros gerais".
whuber
Quero dizer modelos como o vonbertalanffy, logístico ou exponencial duplo, por exemplo. Em todos os casos, o modelo é linear em alguns parâmetros e não linear em outros. As pessoas geralmente tentam transformá-las para obter parametrizações mais estáveis, concentrando-se nos parâmetros não lineares. No entanto, esta é a abordagem errada. É a parametrização linear que deve ser modificada. Por exemplo, para uma logística de 4 parâmetros, o modelo é linear na assíntota superior e inferior, mas, em vez de usar esses parâmetros, deve-se usar os valores previstos para o menor e o maior ind. var.
Dave fournier
@davefournier Obrigado por responder e apontar para o seu artigo. Seu trabalho parece um pouco difícil de ser lido, mas a técnica parece interessante, portanto, mal pode esperar para lê-lo.
Ponto fixo
2

Se você precisar fazer isso várias vezes, sugiro que você use um algoritmo evolutivo na função SSE como um front-end para fornecer os valores iniciais.

Por outro lado, você pode usar o GEOGEBRA para criar a função usando controles deslizantes para os parâmetros e brincar com eles para obter valores iniciais.

OU os valores iniciais dos dados podem ser estimados por observação.

  1. D e E vêm da inclinação e interceptação dos dados (ignorando o gaussiano)
  2. A é a distância vertical do máximo do Gaussiano a partir da estimativa da linha Dx + E.
  3. B é o valor x do máximo do valor gaussiano
  4. C é metade da largura aparente do Gaussiano
Adrian O'Connor
fonte
1

Para valores iniciais, você pode fazer um ajuste mínimo de quadrados. Sua inclinação e interceptação seriam os valores iniciais para D e E. O maior residual seria o valor inicial para A. A posição do maior residual seria o valor inicial para B. Talvez alguém possa sugerir um valor inicial para sigma.

No entanto, mínimos quadrados não lineares sem derivar qualquer tipo de equação mecanicista do conhecimento do assunto são negócios arriscados, e realizar muitos ajustes separados torna as coisas ainda mais questionáveis. Existe algum conhecimento sobre o assunto por trás da equação proposta? Existem outras variáveis ​​independentes relacionadas às diferenças entre os 100 ajustes separados? Pode ser útil incorporar essas diferenças em uma única equação que ajustará todos os dados de uma só vez.

Emil Friedman
fonte