Jak interpretować macierz kowariancji z dopasowania krzywej?

15

Nie jestem zbyt dobry w statystyce, więc przepraszam, jeśli to proste pytanie. Dopasowuję krzywą do niektórych danych, a czasami moje dane najlepiej pasują do ujemnego wykładniczego w postaci zami(-bx)+do , a czasami dopasowanie jest bliższe . Czasami jednak oba z nich zawodzą i chciałbym wrócić do liniowego dopasowania. Moje pytanie brzmi: w jaki sposób mogę ustalić, który model najlepiej pasuje do określonego zestawu danych na podstawie wynikowej macierzy wariancji-kowariancji zwróconej z funkcji scipy.optimize.curve_fit () ? Uważam, że wariancja dotyczy jednej z przekątnych tej macierzy, ale nie jestem pewien, jak to interpretować.zami(-bx2))+do

AKTUALIZACJA: W oparciu o podobne pytanie mam nadzieję, że macierz wariancji-kowariancji może mi powiedzieć, który z trzech modeli najlepiej próbuję dopasować do danych (staram się dopasować wiele zestawów danych do jednego z tych trzech modeli).

Wynikowe macierze wyglądają następująco dla podanego przykładu:

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]]

Oto przykład tego, co robię:

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
źródło
Wspaniale jest, że łączysz się z tym pytaniem CV, a w konsekwencji z ważnym wątkiem komentarza (czarno-biały rolando2, Frank Harrell, ...), pytającym, czy właściwe jest wybranie modelu post facto na podstawie dopasowania. Być może lepiej jest skorzystać z wcześniejszej wiedzy o systemie, aby wybrać model.
Aman
To inne pytanie dotyczące CV może być pomocne: stats.stackexchange.com/questions/50830/…
Aman
Czy może to pomóc zrozumieć, jak interpretować macierz współzależności stats.stackexchange.com/questions/10795/… - Powiedziałbym, że wartość trzeciego modelu jest mniejsza, co wskazuje na mniejsze odchylenie.
user4581,

Odpowiedzi:

4

Dla wyjaśnienia, zmienna pcovfrom scipy.optimize.curve_fitjest szacunkową kowariancją oszacowania parametru, czyli luźno mówiąc, biorąc pod uwagę dane i model, ile informacji jest w danych, aby określić wartość parametru w danym modelu. Więc tak naprawdę nie mówi ci, czy wybrany model jest dobry, czy nie. Zobacz także to .

Problem, jakim jest dobry model, jest rzeczywiście trudnym problemem. Jak argumentowali statystycy

Wszystkie modele są błędne, ale niektóre są przydatne

Kryteria stosowane w porównaniu różnych modeli zależą od tego, co chcesz osiągnąć.

Na przykład, jeśli chcesz, aby krzywa była „jak najbliższa” danym, możesz wybrać model, który daje najmniejszą resztę . W twoim przypadku byłby to model funci oszacowane parametry, poptktóre mają najniższą wartość podczas obliczeń

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

Jeśli jednak wybierzesz model z większą liczbą parametrów, wartość resztkowa automatycznie spadnie , kosztem większej złożoności modelu. A więc wraca do tego, jaki jest cel modelu.

hakanc
źródło