Podstawy analizy danych

Wstęp

Podczas badań naukowych (i nie tylko) wielokrotnie napotykamy podobny schemat działania z danymi, które są uzyskane z różnych źródeł (z eksperymentu, z baz danych, publikacji itp.). Zwykle sprowadza się on do trzech podstawowych kroków: wczytania danych, ich obróbki (analizy) oraz stworzenia informacji wyjściowych w postaci nowych danych, tabel, wykresów itp. W poniższych przykładach spróbujemy przejść przez te etapy w bardzo uproszczonym schemacie.

Wczytywanie danych

Dane, które uzyskujemy są w bardzo rozmaitych formatach. Mogą to być dane tesktowe (np. w formacie ASCII, CSV), binarne (w jakimś szerszym standardzie lub specyficznym, np. dla używanych urządzeń) ewentualnie pochodzące z komunikacji z bazą danych (np. MySQL). W przykładach skupimy się na dwóch pierwszych typach, ponieważ są one najczęściej spotykane w czasie pracy w laboratoriach.

Dane tekstowe

Przykładem danych tekstowych jest plik uzyskany z analizatora Tukan, m.in. używanego na kilku stanowiskach na zaawansowanej pracowni fizyki w ćwiczeniach z fizyki jądrowej. Plik został zebrany podczas pomiaru kalibracyjnego ze źródłem $^{152}$Eu i ma postać następującą

$SPEC_ID:
Ge
$SPECTRUM_DESC:

$DATE_MEA:
3/13/2023 09:21:06 AM
$MEAS_TIM:
2069 s   2111 s
$DATA:
    0   8192
         0              0
         1              0
         2              0
         3              0
(...)
      8189              0
      8190              0
      8191              0

Na początku znajduje się nagłówek zawierający informacje o dacie, czasie trwania pomiaru i podobnym rzeczom. Dane zaczynają się w tym przypadku od linii 11 zawierającej 0 0. Mamy tu dwie kolumny, pierwsza to numer kanału, druga liczba zliczeń zarejestrowana w danym kanale. Jest to przykład widma uzyskanego ze spektrometru promieniowania gamma, ale wiele tego i innego typu urządzeń jako wynik pomiaru zwraca nam jakieś dane w postaci tabel.

Zaczniemy od najprostszej metody polegającej na ręcznym wczytaniu pliku i zinterpretowaniu jego zawartości. Nasza funkcja load_tukan będzie brała informacje o pliku do wczytania oraz o długości nagłówka, który w ogólności może się zmieniać. Następnie otworzymy plik w trybie tekstowym

datafile = open(filename, 'r')

i wczytamy go linia po linii

for line in datafile:
    ...

Opuszczamy linie tak długo, aż miniemy długość nagłówka. Prawdziwe dane zawierają dwie (lub jedną w innej wersji plików) kolumny z liczbami całkowitymi. Wczytana linia to napis typu string. Każdą linię oczyścimy z białych znaków na początku i na końcu i podzielimy na kolejne wyrazy w kolumnach

line.strip().split()

Teraz wystarczy przerobić kolejne wyrazy na typ Integer i dodać do tabeli z danymi. Na koniec tabelę przerobimy na macierz numpy.

Cały kod wygląda następująco.

import numpy

def load_tukan(filename, header):
    i = 0
    datafile = open(filename, 'r')
    data = []
    for line in datafile:
        i += 1
        if i <= header:
            continue
        row = line.strip().split()
        for i in range(len(row)):
            row[i] = int(row[i])
        data.append(row)
    return numpy.array(data)


data = load_tukan('eu152.lst', 11)
print(data)

Ten sam efekt można osiągnąć używając funkcji numpy.loadtxt. Ale nie zawsze pliki są skonstruowane w przyjazny sposób i automatyczne wczytywanie może zawieść. Na przykład nie poradzi sobie z brakującymi danymi, które są zastąpione jakimś znakiem tekstowym (np. -, ?), ponieważ oczekuje kolumn w dobrze ustalonym typie. Dlatego tego typu funkcja stworzona krok po kroku, czasem jest niezbędna. W naszym przykładzie numpy bez problemu jednak wczyta dane.

data = numpy.loadtxt("eu152.lst", skiprows=10)
print(data)

W obu wersjach powinniśmy móc już narysować widmo

import matplotlib.pyplot as plt

plt.plot(data[:, 0], data[:, 1], ds='steps-mid')
plt.show()

Dane binarne

Format binarny oznacza, że dane są zapisywane do pliku w takiej postaci w jakiej występują w pamięci komputera (reprezentacji bitowej), a nie w postaci tekstowej. Takie pliki trudno czytać i interpretować, jeżeli nie znamy dokładnie struktury danych, ale za to zajmują znacznie mniej miejsca. Zapisanie jakiejś liczby w pliku tekstowym zajmuje tyle bajtów ile ma ona znaków, ale jeżeli użyjemy np. reprezentacji Int32, to zajmie ona zawsze 4 bajty niezależnie od długości i nie będzie potrzeby dodawania spacji, przecinków czy innych znaków oddzielających wartości, które także zajmują miejsce. Czasami pliki używają kompresji danych, specjalnych rozmiarów pól lub innych sztuczek i wtedy byłoby dobrze, żeby istniała biblioteka obsługująca dany format. Popularnym formatem binarnym z wieloma możliwościami obsługiwanym przez wiele języków jest na przykład HDF5 (ale nie będziemy się nim tu zajmować).

W naszym przykładzie plik zawiera dane z analizatora posiadającego 3 karty po 16 kanałów. Dane są umieszczone w pliku dla kolejnych modułów i dla kolejnych kanałów, czyli: moduł 1: kanał 1, kanał 2, ... kanał 16, moduł 2: kanał 1, ... Każdy kanał zawiera widmo o 32768 kanałach, a każdy kanał to liczba typu UInt32 (4 bajty).

Zaczynamy od otwarcia pliku w trybie binarnym

datafile = open("cs137.bin", "rb")

Dane można czytać poleceniem read podając długość wczytanych informacji w bajtach, np.

datafile.read(4)

wczyta 4 kolejne bajty. Możemy też wczytać jednocześnie całe widmo z danego kanału

datafile.read(4 * 32768)

Wcztane w ten sposób dane na razie są w postaci tablicy bajtów. Aby przerobić je na właściwą postać musimy znać ich postać i użyć polecenia struct.unpack. Postać danych wpisujemy kodami w postaci liter: np. h oznacza liczbę typu Int16, d liczbę Float64 itd. (pełna lista jest w dokumentacji). Nasz typ to I (liczba Int32 bez znaku) powtórzony 32768 razy. Dane przerobione na listę możemy od razu narysować (osobny wykres dla każdego modułu).

import struct
import matplotlib.pyplot as plt

datafile = open("cs137.bin", "rb")
for m in range(3):
    plt.figure(m+1, (8, 6))
    for c in range(16):
        plt.subplot(4, 4, c+1)
        v = struct.unpack('I' * 32768, datafile.read(4 * 32768))
        plt.plot(v, label="{}:{}".format(m, c))
        plt.legend()
datafile.close()

plt.show()

Okazuje się, że wszystkie kanały są puste poza modułem 0, kanał 1 (tylko jeden detektor był podłączony do urządzenia).

Dopasowywanie krzywych

Jedną z często wykonywanych czynności jest dopasowanie krzywej (funkcji) do danych w celu ich opisania, porówania z modelem, czy określenia pewnych parametrów. Na początek spróbujemy dopasowania prostej do sztucznie stworzonych danych, żeby zobaczyć jak to się robi. Potem spróbujemy dopasować funkcję Gaussa do linii we wczytanych wcześniej danych.

Zaczniemy od "wyprodukowania" danych. Zrobimy zestaw punktów układających się wzdłuż prostej o znanych parametrach, ale aby było trochę ciekawiej, dodamy niewielki losowy rozrzut (z rozkładu normalnego).

xdata = numpy.linspace(0, 10, 11)
ydata = 0.2 * xdata + 1.5 + 0.1 * numpy.random.normal(size=xdata.size)

Najważniejszą częścią jest oczywiście dopasowanie. Funkcja curve_fit z pakietu scipy.optimize używa metody nieliniowej najmniejszych kwadratów. Musimy podać funkcję, którą chcemy opisać dane, listę położeń x, y, oraz początkowe wartości (których oczywiście musimy podać tyle ile bierze funkcja). Dla bardziej skomplikowanych funkcji musimy wybrać w miarę sensowne wartości, w pobliżu oczekiwanych wartości. W innym przypadku algorytm szukający minimum (domyślnie: Levenberga–Marquardta) potrafi zdryfować w inne lokalne minimum. Funkcja curve_fit zwraca dwie macierze, pierwsza zawiera znalezione parametry, druga jest macierzą kowariancji, która zawiera kowariancję dla wszystkich par parametrów. Z definicji kowariancji wiemy, że na diagonalii, gdzie mamy kowariancję paremetrów samych z sobą, dostaniemy kwadraty odchyleń standardowych. W uproszczonej analizie często jest to wystarczająca informacja utożsamiana z niepewnością (co jest dobrym przybliżeniem, o ile pozadiagonalne elementy są małe).

import numpy
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

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

xdata = numpy.linspace(0, 10, 11)
ydata = 0.2 * xdata + 1.5 + 0.1 * numpy.random.normal(size=xdata.size)
plt.plot(xdata, ydata, 'o')

popt, pcon = curve_fit(f, xdata, ydata, p0=[1, 0])

print('a =', popt[0], '+/-', numpy.sqrt(pcon[0][0]))
print('b =', popt[1], '+/-', numpy.sqrt(pcon[1][1]))

plt.plot(xdata, f(xdata, *popt), '--')

plt.tight_layout()
plt.show()

W dalszej części zabawy użyjemy dane wczytywane w pierwszym przykładzie. W okolicach kanału 270 znajduje się tam największa linia. Linie w tego typu widmach mają kształt krzywej Gaussa, ale dodatkowo wokół nich znajduje się tło. Dla w miarę niewielkiego obszaru możemy je przybliżyć funkcją liniową.

data = numpy.loadtxt("eu152.lst", skiprows=10)
plt.subplot(1, 2, 1)
plt.plot(data[:, 0], data[:, 1],  ds="steps-mid")

plt.subplot(1, 2, 2)
plt.plot(data[:, 0], data[:, 1],  ds="steps-mid")
plt.xlim(250, 300)
plt.ylim(0, None)
plt.tight_layout()
plt.show()

Zdefinujemy funkcję N zawierającą funkcję Gaussa i liniowe tło. Do danych w wybranym zakresie trzeba wybrać rozsądne początkowe parametry. Położenie środka łatwo odczytać z wykresu (270), wartość odchylenia standardowego oszacujemy na około 2 (szerokość linii w połowie wysokości jest równa mniej więcej 2 razy $\sigma$) i wreszcie parametr A oznaczający całkę pod krzywą ustalimy na równy dwa razy liczbie zliczeń w największym kanale (co powinno być wystarczające bliskim). Tło spróbujemy ustalić na stałe a0 = 2000, a1 = 0.

import numpy
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def N(x, x0, s, A, a0, a1):
    return (A / numpy.sqrt(2 * numpy.pi * s**2)
            * numpy.exp(-(x - x0)**2 / (2 * s**2)) + a0 + a1 * x)

data = numpy.loadtxt("eu152.lst", skiprows=10)
plt.plot(data[:, 0], data[:, 1],  ds="steps-mid")

x0 = 260
x1 = 280

popt, pcon = curve_fit(N, data[x0:x1, 0], data[x0:x1, 1], 
        p0=[270, 2.0, 240000, 2000, 0])

names = ['x0', 's', 'A', 'a0', 'a1']

for i, name in enumerate(names):
    print('{} = {:.3f} +/- {:.3f}'.format(name, popt[i], 
                                        numpy.sqrt(pcon[i][i])))

xf = numpy.linspace(x0, x1, 1000)
plt.plot(xf, N(xf, *popt), '-')
plt.xlim(x0 - 10, x1 + 10)

plt.tight_layout()

plt.show()

Można jeszcze z ciekawości sprawdzić, co się stanie, gdy parametry początkowe wybierzemy bardzo źle.

popt, pcon = curve_fit(N, data[x0:x1, 0], data[x0:x1, 1], 
        p0=[1.0, 2.0, 3.0, 4.0, 5.0])

for i, name in enumerate(names):
    print('{} = {:.3f} +/- {:.3f}'.format(name, popt[i], 
                                        numpy.sqrt(pcon[i][i])))

xf = numpy.linspace(x0, x1, 1000)
plt.plot(data[:, 0], data[:, 1],  ds="steps-mid")
plt.plot(xf, N(xf, *popt), '-')
plt.xlim(x0 - 10, x1 + 10)

plt.tight_layout()

plt.show()

Ale wystarczy podać wartości nieco lepsze wartości niektórych z nich (x0, A), żeby rozwiązanie zbiegło do właściwego minimum.

popt, pcon = curve_fit(N, data[x0:x1, 0], data[x0:x1, 1], 
        p0=[270.0, 1.0, 1000.0, 0.0, 0.0])

for i, name in enumerate(names):
    print('{} = {:.3f} +/- {:.3f}'.format(name, popt[i], 
                                        numpy.sqrt(pcon[i][i])))

xf = numpy.linspace(x0, x1, 1000)
plt.plot(data[:, 0], data[:, 1],  ds="steps-mid")
plt.plot(xf, N(xf, *popt), '-')
plt.xlim(x0 - 10, x1 + 10)

plt.tight_layout()

plt.show()

Przedstawiona tu analiza używa otwartych przedziałów parametrów. Można także dopasowywać parametry przez curve_fit, innym algorytmem (pod "spodem" ukrywa się algorytm trust-region-reflective), który pozwala na wpisanie pewnych przedziałów, wynikających np. z fizycznych lub logicznych ograniczeń. W przykładzie dopasowywania linii podamy pewną listę rozsądnych warunków:

• środek linii musi się znaleźć w zadanym przedziale (pomiędzy x0 i x1)

• szerokość linii (sigma) jest pomiędzy 0, a szerokością przedziału

• pole powierzchni linii (A) jest pomiędzy 0 i całkowitą liczbą zliczeń w zadanym przedziale (data[x0:x1, 1].sum())

• parametry liniowego tła są dowolne (pomiędzy $\pm\infty$)

Warunki są podawane w postaci tabeli zawierającej listę dolnych ograniczeń dla poszczególnych parametrów i drugą listę górnych ograniczeń. W tym przypadku nie musimy podać początkowej estymacji parametrów.

popt, pcon = curve_fit(N, data[x0:x1, 0], data[x0:x1, 1], 
        bounds=[(x0, 0, 0, -numpy.inf, -numpy.inf), (x1, x1-x0, data[x0:x1, 1].sum(), numpy.inf, numpy.inf)])

for i, name in enumerate(names):
    print('{} = {:.3f} +/- {:.3f}'.format(name, popt[i], 
                                        numpy.sqrt(pcon[i][i])))

xf = numpy.linspace(x0, x1, 1000)
plt.plot(data[:, 0], data[:, 1],  ds="steps-mid")
plt.plot(xf, N(xf, *popt), '-')
plt.xlim(x0 - 10, x1 + 10)

plt.tight_layout()

plt.show()

Podane granice przedziałów nie były bardzo restrykcyjne, ale jak widać w zupełości wystarczyły, aby algorytm zbiegł do właściwego rozwiązania.

Jak można się spodziewać, przedstawione tu przykłady to wierzchołek góry lodowej dotyczącej analizy danych, ale mam nadzieję, że ułatwi zaczęcie jej poznawania.