Como interpreto a matriz de covariância a partir de um ajuste de curva?

15

Eu não sou muito bom em estatística, então peço desculpas se esta for uma pergunta simplista. Estou ajustando uma curva de alguns dados, e às vezes os meus dados melhor se encaixa uma exponencial negativa na forma , e às vezes o ajuste é mais perto de um * e ( - b * x 2 ) + c . No entanto, às vezes os dois falham e eu gostaria de voltar a um ajuste linear. Minha pergunta é: como posso determinar qual modelo se ajusta melhor a um determinado conjunto de dados da matriz de variância-covariância resultante retornada doumae(-bx)+cumae(-bx2)+cfunção scipy.optimize.curve_fit () ? Acredito que a variação esteja em uma das diagonais dessa matriz, mas não sei como interpretar isso.

ATUALIZAÇÃO: Com base em uma pergunta semelhante , espero que a matriz de variância-covariância possa me dizer qual dos três modelos que estou tentando melhor se ajusta aos dados (estou tentando ajustar muitos conjuntos de dados a um desses três modelos).

As matrizes resultantes são assim para o exemplo dado:

pcov_lin 
[[  2.02186921e-05  -2.02186920e-04]
 [ -2.02186920e-04   2.76322124e-03]]
pcov_exp
[[  9.05390292e+00  -7.76201283e-02  -9.20475334e+00]
 [ -7.76201283e-02   6.69727245e-04   7.90218415e-02]
 [ -9.20475334e+00   7.90218415e-02   9.36160310e+00]]
pcov_exp_2 
[[  1.38338049e-03  -7.39204594e-07  -7.81208814e-04]
 [ -7.39204594e-07   8.99295434e-09   1.92970700e-06]
 [ -7.81208814e-04   1.92970700e-06   9.14746758e-04]]

Aqui está um exemplo do que estou fazendo:

import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
import scipy.optimize

def exp_func(x, a, b, c):
    return a * np.exp(-b * x) + c

def exp_squared_func(x, a, b, c):
    return a * np.exp(-b * x*x*x) + c

def linear_func(x, a, b):
    return a*x + b

def main():
    x = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20], np.float)
    y = np.array([1, 1, 1, 1, 0.805621, 0.798992, 0.84231, 0.728796, 0.819471, 0.570414, 0.355124, 0.276447, 0.159058, 0.0762189, 0.0167807, 0.0118647, 0.000319948, 0.00118267, 0, 0, 0], np.float)

    p0 = [0.7746042467213462, 0.10347274384077858, -0.016253458007293588]
    popt_lin, pcov_lin      = scipy.optimize.curve_fit(linear_func, x, y)
    popt_exp, pcov_exp      = scipy.optimize.curve_fit(exp_func, x, y)
    popt_exp_2, pcov_exp_2  = scipy.optimize.curve_fit(exp_squared_func, x, y)

    plt.figure()
    plt.plot(x, y, 'ko', label="Original data")
    plt.plot(x, linear_func(x, *popt_lin), 'r-', label='linear')
    plt.plot(x, exp_func(x, *popt_exp), 'b-', label='exponential')
    plt.plot(x, exp_squared_func(x, *popt_exp_2), 'g-', label='exponential squared')
    plt.legend()
    plt.show()

if __name__ == '__main__':
    main()
Jason Martens
fonte
É ótimo que você faça um link para essa pergunta do CV e, consequentemente, para o importante tópico de comentários (b / w rolando2, Frank Harrell, ...) questionando se é apropriado escolher o modelo pós-fato com base no ajuste. Talvez seja melhor usar o conhecimento prévio do sistema para escolher o modelo.
Aman
Esta outra pergunta no CV pode ser útil: stats.stackexchange.com/questions/50830/…
Aman
Isso pode ser útil para entender como interpretar a matriz de co-variância stats.stackexchange.com/questions/10795/… - Eu diria que o valor do terceiro modelo é menor, indicando menor desvio.
user4581

Respostas:

4

Como um esclarecimento, a variável pcovde scipy.optimize.curve_fité a covariância estimada da estimativa de parâmetro, ou seja, falando livremente, dados os dados e um modelo, quanta informação há nos dados para determinar o valor de um parâmetro no modelo fornecido. Portanto, ele realmente não diz se o modelo escolhido é bom ou não. Veja também isso .

O problema que é um bom modelo é realmente um problema difícil. Como argumentado por estatísticos

Todos os modelos estão errados, mas alguns são úteis

Portanto, os critérios a serem usados ​​na comparação de modelos diferentes dependem do que você deseja alcançar.

Por exemplo, se você quiser uma curva que seja o "mais próximo possível" dos dados, poderá selecionar um modelo que dê o menor resíduo . No seu caso, seria o modelo funce os parâmetros estimados poptque têm o menor valor ao calcular

numpy.linalg.norm(y-func(x, *popt))

No entanto, se você selecionar um modelo com mais parâmetros, o residual diminuirá automaticamente , com o custo de uma complexidade maior do modelo. Então, volta ao que é o objetivo do modelo.

hakanc
fonte