Przykład fitowania nieliniowego oraz sporzadzania wykresów ze słupkami błędów

Directory[]

SetDirectory["L:/home/andrzej/"]

dane   =  
        Import["tmp/dane/Si_errorbars.txt", "Table"];

Przypominam, że do elementów listy odwołujemy się za pomocą podwójnych nawiasów kwadratowych

dane[[1]]

dane[[5]][[2]]

Pakiet Graphics zawiera liczne instrukcje graficzne m. in. do rysowania słupków błędów

<< Graphics`

Zmienna wykres zawiera tablicę która będzie użyta jako punkty wykresu (kolumny 1 i 2 z danych). Pierwszy wiersz (o numerze 1) został pominięty, gdyż zawiera punkt (0.0, 0.0) na którym załamuje się algorytm fitowania nieliniowego .

pomiary  =  
        Table[{dane[[i]][[1]], dane[[i]][[2]]}, {i, 
            1, 100}];

Tu przepisujemy kolumnę 3 zawierającą błędy:

bledy  =  
        Table[dane[[i]][[3]], {i, 1, 100}];

Tak wyglądają dane:

wykres1 = ListPlot[pomiary]

Kazdy wykres mozna zapisac w fowolnym formacie np. PNG. (Trzeba usunąć komentarz (* ... *) !)

Niżej podano sposób sporządzenia standardowego wykresu ze słupkami błędu:

dane1  =  
        Table[{{dane[[i]][[1]], dane[[i]][[2]]}, 
            ErrorBar[dane[[i]][[3]]]}, {i, 1, 100}];

wykres2  = 
      MultipleListPlot[dane1, PlotRange -> {{0, 6}, {0, 0.7}}, 
        SymbolShape -> {PlotSymbol[Triangle, 1.5]}]

Export["tmp/w1.png",  %,  "PNG", 
      ImageResolution -> 300]

W celu dofitowania nieliniowej funkcji ładujemy pakiet Statistics

<< Statistics`

Teraz mozemy użyć funkcji NonlinearFit aby dopasować naszą nieliniową w parametrach funkcję. Instrukcja Drop[pomiary,1]
wyrzuca pierwszy punkt  {0.0,0.0,0.0} powodujący błąd algorytmu fitującego (sprawdz co się dzieje bez tej funkcji
aby ,,oswoić'' się z komunikatami o błędach)

pomiary1  =  Drop[pomiary, 1];

wykres3  =  
      NonlinearFit[pomiary1, a x^b Exp[-c x], x, {a, b, c}]

Analogicznie jak w przypadku funkcji Fit / Regress mamy 2 funkcje NoninearFit / NonlinearRegress różniące się szczegółowością wyświetlanej analizy.

NonlinearRegress[pomiary1, a x^b Exp[-c x], 
      x, {a, b, c}]

Tu fitowanie uwzględniające błędy:

NonlinearRegress[pomiary1, a x^q Exp[-b x], 
      x, {{a, 50.1, 200}, {b, 1, 3}, {q, 2, 4}}, 
      Weights -> Drop[bledy, 1]]

wykres4  =  
      NonlinearFit[pomiary1, a x^q Exp[-b x], 
        x, {{a, 50.1, 200}, {b, 1, 3}, {q, 2, 4}}, 
        Weights -> Drop[bledy, 1]]

Fit biorący pod uwagę błędy (wykres4) wygląda prawie identycznie jak ten z równymi wagami (wykres3)

wykres5  =  
      Plot[{wykres3,  wykres4}, {x, 0, 6}, 
        PlotStyle -> {Red, Blue}]

Jak widać funkcje te różnią się:

Plot[(wykres3 - wykres4), {x, 0, 10}, 
      PlotRange -> All]

I wreszcie fit razem z punktami "pomiarowymi"

Show[wykres2, wykres5]

Odczytujemy dofitowane parametry a,b,c ("Estimate") oraz ich błąd standardowy ("Asymptotic SE")
zapisane w tablicy "ParameterCITable"  generowanej przez NonlinearRegress

parametry = ParameterCITable/. NonlinearRegress[
          pomiary1, a x^b Exp[-c x], 
          x, {a, b, c}]

aFit  =  parametry[[1]][[1]][[1]]
aBlad  =  parametry[[1]][[1]][[2]]
bFit  =  parametry[[1]][[2]][[1]]
bBlad  =  parametry[[1]][[2]][[2]]
cFit  =  parametry[[1]][[3]][[1]]
cBlad  =  parametry[[1]][[3]][[2]]

Na koniec dodajemy kilka szczegółów aby wyniki prezentowaly sie elegancko. Symbol "<>" łączy łańcuchy tekstowe, natomiast ToString[] zamienia wyrażenie matematyczne na tekst. NumberForm[] wypisuje zadana liczbe cyfr znaczacych, co pozwoli uniknac poprawiania sprawozdania ;-) Prowadzący irytuja sie widzac ze a=16.7275 podczas gdy blad wynosi ponad 0.2; pozostałe cyfry "a" nie niosą żadnej informacji a  uzyskane wyniki mogą stać się nieczytelne.

\(\(\(Opis\  = \ \*""Y= a (X^b) \!\(e\^\(\(-c\)\\ X\)\) \n"" <> "a=" <> ToString[NumberForm[aFit, 3]] <> 
        ToString[\(±NumberForm[aBlad, 
              1]\)] <> "
"\  <> \ "b=" <> 
        ToString[NumberForm[bFit, 3]] <> 
        ToString[\(±NumberForm[bBlad, 
              1]\)] <> "
" <> 
		"c=" <> 
        ToString[NumberForm[cFit, 3]] <> 
        ToString[\(±NumberForm[cBlad, 
              1]\)]\)\(
\)
    \)\)

final  = 
      Show[wykres2, wykres5, Frame -> True, 
        FrameLabel -> {"X", "Y", "Pomiary + fit", ""}, 
        FrameTicks -> {Automatic, {0.2, 0.4, 0.6}, None, None}, 
        Epilog -> 
          Text[Opis, {5, 0.5}, Background -> GrayLevel[0.7]], 
        TextStyle -> {FontFamily -> "Times", 
            FontSize -> 10}]

Wreszcie zapisujemy otrzymany wykres w kilku formatach.

Export["tmp/final.eps", 
      final, "EPS"]
Export["tmp/final.jpeg", final, "JPEG", 
      ImageResolution -> 300]
Export["tmp/final.png", final, "PNG", 
      ImageResolution -> 300]