Generatcja danych

Główną wadą prowadzenia obliczeń numerycznych jest to, że nie zawsze wiadomo jaki powinien być wynik. W związku z tym stosuje się bardzo prostą metodę: Jako danych wejściowych do algorytmu używa się danych „syntetycznych". Wygenerowanych w sposób algorytmiczny i takich, dla których można przewidzieć wynik. Bardzo często dane są zaburzane. Najczęściej zaburzenie jest generowane (pseudo)losowo i znamy jego rozkład, średnią (ta, bardzo często, równa jest zero) i wariancję.

Bardzo często będziemy potrzebowali danych w formie:

$$y_i=f(x_i)+\zeta_i,\quad i=1,2,\ldots,N; \quad \zeta_i, \mbox{— losowe}$$

Zazwyczaj (albo najprościej) $x_i=x_0+i\Delta$, gdzie $\Delta=(x_N - x_0)/N$ (czasami dopuszczamy, żeby $i$ zaczynało się od zera).

Aby wygenerować wartość losową — należy skorzystać z generatora liczb pseudolosowych dostarczanych przez używany język programowania. W przypadku C może to być GSL — GNU Scientific Library¹. Zarówno matlab jak i Mathematica mają swoje własne biblioteki.

Mathematica

1. Generacja $x_i=x_0+i\Delta$. Można w tym celu użyć polecenia Table:

   data1 = Table[x,{x0,xN,Delta}]
   

   na przykład

   data1 = Table[i, {i,0,10,0.1}]
   

   Zmienna (złożona) data1 zawierać będzie 101 wartości z zakresu od 0 do 10.

2. Generacja $f(x)$.

   • najpierw definiujemy funkcję; W języku Mathematici nie jest to specjalnie skomplikowane:

     f[x_] := x + Sin[x]
     

     (zwracam uwagę na zapis x_; oznacza on, że $x$ jest zmienną niezależną funkcji oraz symbol := oznaczajacy definicję, a nie podstawienie). Cechą charakterystyczną funkcji jest to, że mogą działać ona na zmiennych (stałych) prostych:

     f[4]
     4 + Sin[4]
     

     (mathematica, tak długo jak się da, liczy symbolicznie)

     N[f[4]]
     3.2432
     

     również, gdy argumentem jest lista

     N[f[{3, 4}]]
     {3.14112, 3.2432}
     

   • następnie nakładamy funkcję na nasze dane:

     data2 = f[data1]
     

   • możemy (przy okazji) zobaczyć wykres:

     ListLinePlot[Transpose[{data1, data2}]]
     

     (funkcja Transpose[] oznacza transpozycję; transpozycja jest potrzebna ze względu formę danych dla polecenie ListLinePlot; macierz utworzona poleceniem ${\text{data1},\text{data2}}$ to dwie kolumny $x$-ów i $y$-ków. ListLinePlot wymaga wektora, którego elementami są pary $(x_i,y_i)$; transpozycja to załatwia.)

     

3. Wreszcie możemy nałożyć zaburzenia. Zazwyczaj o zaburzeniach będziemy zakładać że mają rozkład normalny² o średniej zero i niezbyt wielkim odchyleniu standardowym. Funkcja RandomVariate pozwala wygenerować zadaną liczbę danych o wskazanym rozkładzie. Pierwszym argumentem funkcji jest rodzaj rozkładu, drugim — ilość liczb.

   data3 = RandomVariate[NormalDistribution[], 101]
   

   (brak parametrów oznacza średną zero i wariancję 1, rozkładu normalnego). Łatwo zgadnąć, że data2 + data3 to będzie wektor zaburzonych wartości naszej funkcji, czyli $f(x)+\zeta$. Poniżej wykres, przy czym dane zostały zaburzone zmienną o rozkłądzie $N(0,0.1)$:

   

Matlab

W matlabie wszystko praktycznie tak samo, tylko (chyba) nieco prościej:

1. Generacja $x_i=x_0+i\Delta$. Można w tym celu użyć operatora „:" w formie Start:przyrost:Koniec:

   data1=0:0.1:10;
   

   (gdy nie dodamy średnika, matlab wypisze wszystkie wartości.

2. Teraz funkcja. Ze zdefiniowaniem funkcji jest pewien kłopot. Funkcje powinny znajdować się w plikach zewnętrznych, ale można spróbować czegoś takiego³:

   f=@(x) sin(x)+x;
   

   Aby sprawdzić czy działa piszemy f(data1) albo data2=f(data1).

3. Generacja liczb losowych odbywa się za pomocą funkcji random. Wywołanie jej jest proste:

   data3 = random(name,A,B,...)
   

   gdzie (name) to nazwa rozkładu (dla rozkładu normalnego będzie to norm albo Normal), A, B,… to parametry rozkładu. Parametr może być jeden, może nie być go wcale, w przypadku rozkładu normalnego parametry powinny być dwa: średnia rozkładu i wariancja. Uwaga: nazwę rozkładu podajemy jako tekst.

   random('Normal', 0.0, 1.0)
   

   wygeneruje jedną wartość. Na końcu wszystkich parametrów może pojawić się dodatkowy parametr mówiący ile wartości generować (czy raczej jaką tablicę wartościami wypełnić).

   data3 = random('Normal', 0.0, 0.1, [1,101])
   

   Żeby zobaczyć dane oryginalne i zaburzone (data2+data3) używamy polecenia plot:

   plot(data1, data2)
   plot(data1, data2+data3)
   

   i uzyskujemy: