Analiza danych pomiarowych

Regresja liniowa

Przykład 1: dopasowanie linii prostej z użyciem FindFit

dane=Table[{x, 2 x +3 + RandomReal[{-1,1}]},{x,Sort@RandomReal[{-12,12},7]}]

w1=ListPlot[dane]

Manipulate[Show[w1,Plot[a x + b,{x,-12,12}]],{a,-5,5},{b,-5,5}]

fit1=FindFit[dane, a x +b, {a,b},x]

w2=Plot[a x+b/.fit1,{x,-12,12}]

Show[w1,w2]

Przykład 2: dopasowanie linii prostej z użyciem LinearModelFit

Bardziej zaawansowane informacje (błędy współczynników, korelacji itp.) można uzyskać za pomocą polecenia LinearModelFit:

model=LinearModelFit[dane,x,x]

Polecenie powyższe, ale także wiele innych, w tym zintegrowane bazy danych, produkują strukturę danych, zawierającą sporo informacji. Listę rekordów ukrytych w strukturze można wyświetlić korzystając ze standardowego pola “Properties”, zawierającego listę pozostałych:

model["Properties"]

W zależności od potrzeby, wyszukujemy na liście interesujące nas dane, powiedzmy “ParameterTable”, które to pole zawiera tabelkę dopasowanych parametrów:

model["ParameterTable"]

Podane błędy są na poziomie ,,1σ’’ ( czyt: jeden sigma ), co oznacza, iż prawdopodobieństwo, że współczynniki prostej różnią się od podanych bardziej niż podany błąd, wynosi około 68%. Przypominam, że dane zostały ,,wygenerowane’’ dla prostej 2 x+3.

Jeżeli interesuje nas wiele pól, lub nie wiemy pod jaką nazwą ukrywają się interesujące nas dane, warto postąpić następująco. Zapisujemy listę pół do pewnej zmiennej, np: p:

p=model["Properties"]

Obliczmy ile pól jest dostępnych:

num=Length[p]

Nastepnie, użyjemy Manipulate:

Manipulate[{p[[ii]],model[p[[ii]]]},{ii,1,num,1}]

Natrafiamy np. na ciekawe pole ,,MeanPredicitionBands’’, pokazujące jak daleko możemy ekstrapolować dofitowaną prostą, i ile rzeczywiście poza zbadanym zakresem jest to linia prosta:

model["MeanPredictionBands"]

w3=Plot[model["MeanPredictionBands"],{x,-12,12}]

Show[w1,w3]

Na koniec trochę kosmetyki:

wsp=model["BestFitParameters"]

err=model["ParameterErrors"]

txt0="y = a x + b"

txt1="b = "<>ToString[NumberForm[wsp[[1]],3]]<>"±"<>ToString[NumberForm[err[[1]],2]]

txt2="a = "<>ToString[NumberForm[wsp[[2]],3]]<>"±"<>ToString[NumberForm[err[[2]],2]]

Show[w2,w1,Epilog->Text[txt0<>"\n"<>txt2<>"\n"<>txt1,{7,7}]]

?NumberForm

Przykład 3: fit liniowy do danych z podanymi błędami

Demonstrację fitu z uwzględnieniem błędów ,,pomiarowych’’ rozpoczniemy od przygotowania ,,sztucznych’’ błędów. Liczba punktów pomiarowych:

num=7

Wartości na osi poziomej (losowe):

xTBL=Sort@RandomReal[{-5,5},num]

Tutaj spreparujemy ,,błędy pomiarowe’’. Aby zilutrować wagę pomiarów z małym błedem, pierwszy i ostatni punkt będzie miał ręcznie przypisaną dokładność, natomiats pozostałe losowe błedy z przedziału od -1 do 1:

rnd=Join[{10^-4},RandomReal[{-1,1},num-2],{10^-4}]

Generujemy wartości na osi Y dodając sztuczne błędy:

yTBL=(2#+3&)[xTBL]+err

Za błędy pomiarowe przyjmijmy wartość bezwzględną dodanych losowo lub ręcznie liczb:

err=Abs[rnd]

dane={xTBL,yTBL}//Transpose

w1=ListPlot[dane]

Każdemu z punktów pomiarowych przypisujemy wagę, odwrotnie proporcjonalną do kwadratu błędu. Innymi słowy, punkt tym ważniejszy, im dokładniej został zmierzony:

wagi=1/err^2

Porównajmy wyniki fitu z uwzględnieniem i bez uwzględnienia informacji o błędach:

LinearModelFit[dane,x,x,Weights->wagi]

LinearModelFit[dane,x,x]

LinearModelFit[dane,x,x,Weights->wagi,VarianceEstimatorFunction->(1&)]

Jak widać powyżej, dodanie dwóch punktów ,,zmierzonych’’ bardzo precyzyjnie (błąd (10^(-8))) wystarczyło aby wyznaczyć prostą. Jest to zgodne z intuicją. Bez tej informacji prosta (jej współczynniki) są wyznaczone z dużymi błędami.

model1=LinearModelFit[dane,x,x]

model2=LinearModelFit[dane,x,x,Weights->wagi]

model1["ParameterTable"]

model2["ParameterTable"]

Na koniec trochę kosmetyki. Wykresy wyglądają bardziej elegancko, gdy naniesiemy na nie tzw. słupki błędów. Słuzy do tego pakiet ErrorBarPlots:

Needs["ErrorBarPlots`"]

Table[{{xTBL[[ii]],yTBL[[ii]]},ErrorBar[err[[ii]]]},{ii,1,num}]

w4=ErrorListPlot[%]

w5=Plot[model2["Function"][x],{x,-5,5}]

Show[w4,w5]

Regresja nieliniowa

Przykład 1: dopasowanie wzoru nieliniowego  - regresja symboliczna numerycznego rozwiązania równania Lane-Emdena

Na początek przygotujemy ,,dane pomiarowe’’. Motywację daje brak rozwiązania pewnego równania różniczkowego:

DSolve[{y''[x]+2/x y'[x]+y[x]^5==0,y[0]==1,y'[0]==0},y,x]//AbsoluteTiming

Możemy rozwiązać powyższy problem mumerycznie. Niestety, z powodu dzielenia przez zero najprostszy sposób nie działa:

NDSolve[{y''[x]+2/x y'[x]+y[x]^5==0,y[0]==1,y'[0]==0},y,x]

Dlatego omijamy ten problem zastępując zera ,,bardzo małą liczbą’’. Jak pamiętamy z wcześniejszych wykładów, najmniejszą liczbą z punktu widzenia procesora może być np: $MachineEpsilon, lub bezpieczniej jej pierwiastek Sqrt[$MachineEpsilon]:

rozw=NDSolve[{y''[x]+2/x y'[x]+y[x]^5==0,y[Sqrt[$MachineEpsilon]]==1,y'[Sqrt[$MachineEpsilon]]==0},y,{x,0,8}][[1]]

f=y/.rozw

Szczegóły powyższej procedury nie są w tym momencie istotne. Ważne jest, że posiadamy pewne dane numeryczne, które chcemy analizować:

dane=Table[{x,f[x]},{x,Sort@RandomReal[{0,8},16]}]

Analizę danych zawsze zaczynamy od ich obejrzenia. W tym wypadku punkty układają sie na pewnej krzywej.

w0=ListPlot[dane]

W praktyce studenckiej i badawczej, zwykle znamy wzór (np: z teorii ) jaki mógłby opisywać dane. W powyższym przypadku jest to funkcja zadana wzorem:

Manipulate[Plot[a(1+x^p/b)^q,{x,0,8}],{{a,1},0,2},{b,1,12},{p,-2,2},{q,-2,2}]

W obecnej chwili metody znajdywania współczynników takiej krzywej są bardzo zaawansowane, i automatcznie radzą sobie z typowymi przypadkami. Mimo to dobrą praktyką, jest ręczne znalezienie wstępnych wartości. Często zdaża się, że algorytm dopasuje totalnie bzdurną krzywą, co widać na wykresie natychmiast. Manipulate jest tu idealnym narzędziem:

Manipulate[
Show[
w0,Plot[a(1+x^p/b)^q,{x,0,8}]
]
,{{a,1},0,2},{b,1,12},{p,-2,2},{q,-2,2}]

Poniżej mamy niemiłą niespodziankę: algorytm z jakiegoś powodu usiłuje zbadać przypadek b<0, kiedy to pojawia się liczba ujemna pod pierwiastkiem:

FindFit[dane,a(1+x^p/b)^q,{{a,1.1},{b,3.3},{p,2},{q,-1.7}},x]

W tym przypadku FindFit radzi sobie lepiej, gdy sam dobierze wartości startowe:

fit=FindFit[dane,a(1+x^p/b)^q,{a,b,p,q},x]

w1=Plot[a(1+x^p/b)^q/.fit,{x,0,8}]

Show[w1,w0]

Ponownie, gdy potrzebujemy bardziej szczegółowych danych, np: błedów współczynników dopasowania (a jest to raczej reguła niż wyjatek!), mamy do dyspozycji powiązane polecenie NonlinearModelFit:

model=NonlinearModelFit[dane,a(1+x^p/b)^q,{a,b,p,q},x]

model["ParameterTable"]

Powyższe wyniki silnie sugerują, że faktyczne wartości współczynników to a=1, b=3, p=2 i q=-1/2. Sprawdźmy to:

a(1+x^p/b)^q/.model["BestFitParameters"]

Rationalize[%,0.001]

y''[x]+2/x y'[x]+y[x]^5==0/.y->Function[x,Evaluate[%]]

Simplify[%]

Fitując wzór do ,,danych’’ numerycznych rozwiązaliśmy ściśle r. różniczkowe. Jest to przykład udanej regresji symbolicznej.

Przykład 2: dofitowanie modelu w postaci całki (odległości kosmologiczne)

Należy podkreślić, że program typu Mathematica potrafi operować właściwie dowolnym modelem, który mozna wyrazić w języku matematyki. Dlatego nie jesteśmy w żaden sposób ograniczeni jawnymi wzorami typu Y=a Log[1+Sin[ b X]] itp. Równie dobrze może to być nierozwiązywalne analitycznie równanie różniczkowe, całka niewyrażająca się przez funkcje specjalne i elementarne, obiekt typu Root a nawet DifferentialRoot itp. itd.

Możemy wyobrazić sobie np: wyznaczanie współczynnika oporu powietrza na podstawie pozycji fotografowanej w locie piłki, czy badanie wahadła fitując rozwiązanie nieliniowego równania.

Poniższy przykład pokazuje jak znaleźć parametry modelu komologicznego, dysponując danymi w postaci przesunięcia ku czerwieni z (miara odległości) oraz jasności supernowych typu Ia. Wzór na odległość ma postac całki:

Hold@Integrate[1/Sqrt[(1+zz)^2(1+ Ωm zz)-zz(zz+2)Ωλ],{zz,0,z}]//TraditionalForm

Całka ta jest co prawda możliwa do obliczenia (metoda identyczna jak w artykule P. Macha https://arxiv.org/abs/1205.3480 ),
ale w praktyce zmuszeni jesteśmy użyć metod numerycznych:

Integrate[1/Sqrt[(1+zz)^2(1+ Ωm zz)-zz(zz+2)Ωλ],{zz,0,z}]//AbsoluteTiming (* Nie doczekałem się niczego...*)

Potrzebne dane są dostępne w internecie:

(* Import["https://braeburn.pha.jhu.edu/~ariess/R06/Davis07_R07_WV07.dat"] *)
(*Import["https://supernova.lbl.gov/Union/figures/SCPUnion_mu_vs_z.txt"];*)
Import["https://supernova.lbl.gov/Union/figures/SCPUnion2_mu_vs_z.txt"];
Short[%,20] (* Powyższze nie działa bez dostępu do internetu !!! *)!

Wyrzucany 5 pierwszych linijek, stanowiących nagłówek:

sn= Drop[Import["https://supernova.lbl.gov/Union/figures/SCPUnion2_mu_vs_z.txt","Table"] ,5] ;

sn= Drop[Import["C:\\Dysk Google\\SCPUnion2_mu_vs_z.txt","Table"] ,5] ;

Dane mają postać 4 kolumn, pierwsza to identyfikator, druga to z (wartośc na osi poziomej), trzecia to pewna miara odległości, a czwarta to błąd pomiaru.

sn[[1]]

Kluczowe dane to kolumna 2 i 3:

zm={Transpose[sn][[2]],Transpose[sn][[3]]}//Transpose;

ListPlot[zm]

Sprawdzamy i zapisujemy ilość danych (liczbą pomiarów):

num=Dimensions[zm][[1]]

Wzór konwertujący odległość do standardowych jednostek (Mpc, megaparseki):

zdL=Table[{zm[[ii,1]],(10^(-5+(m/5)))/.m->zm[[ii,2]]},{ii,1,num}]//Sort;

Sprawdzamy jaki zakres z obejmują dane (wartość maks. i min.):

z0=zdL[[1]] //First(* najniższe z *)

z1=zdL //Last //First(* najwyższe z *)

Ładujemy pakiet który umożliwia rysowanie słupków błędów:

Needs["ErrorBarPlots`"]

Błędy pomiaru odległości:

err=Table[((2^(-5+((sn[[ii,3]])/5))) (5^(-6+((sn[[ii,3]])/5))) Log[10])*sn[[ii,4]],{ii,1,num}];

zdLerr=Table[{{zdL[[ii,1]],zdL[[ii,2]]},ErrorBar[err[[ii]]]},{ii,1,num}];

Oto dane pomiarowe wraz z błędami:

ErrorListPlot[zdLerr]

Naszym zadaniem jest ustalenie, jakie wartości parametrów Ωm, Ωλ, H0 w poniższym wzorze najlepiej odtwarzają powyższe dane:

c/H0(1+z)/Sqrt[Ωm+Ωλ-1]Sin[Sqrt[Ωm+Ωλ-1] NIntegrate[1/Sqrt[(1+zz)^2(1+ Ωm zz)-zz(zz+2)Ωλ],{zz,0,z}] ]

Stałe te opisują parametry Wszechświata. Sprawę kompilkuje fakt, iż w rzeczywistości mamy 3 różne wzory (model zamknięty, otwarty i płaski). Proste wstawienie takiego wzoru do FindFit nie zadziała. Trudności mają jednak tylko i wyłącznie charakter techniczny.

c=299792.458 (* prędkość światła w km/s *)

Użycie skomplikowanego modelu wymaga użycia kilku technik (trików):

a) wzór należy zapisać w postaci funkcji parametrów, które chcemy dofitować

b) wynikiem działania funkcji powinna być inna funkcja, tym razem zależna tylko od wartości na osi poziomej

c) parametry powinny mieć dodany wzorzec NumberQ, co pozwala na uniknięcie błedów wywołania funkcji (jest ona obliczana tylko gdy argumentem jest liczba; w trakcie analizy wstępnej/symbolicznej funkcja nie jest obliczana !)

d) w wielu przypadkach użycie funkcji pamiętacjącej ( f[x_]:=f[x]=... ) może drastycznie przyspieszyć obliczenia !

Wygląda to następująco:

dL[Ωm_?NumberQ,Ωλ_?NumberQ,H0_?NumberQ]:=(dL[Ωm,Ωλ,H0]=Interpolation[Table[{z,Piecewise[{{c/H0(1+z)/Sqrt[1-Ωm-Ωλ]Sinh[Sqrt[1-Ωm-Ωλ] NIntegrate[1/Sqrt[(1+zz)^2(1+ Ωm zz)-zz(zz+2)Ωλ],{zz,0,z}] ],Ωm+Ωλ<1},{c/H0(1+z)NIntegrate[1/Sqrt[(1+zz)^2(1+ Ωm zz)-zz(zz+2)Ωλ],{zz,0,z}] ,Ωm+Ωλ==1},{c/H0(1+z)/Sqrt[Ωm+Ωλ-1]Sin[Sqrt[Ωm+Ωλ-1] NIntegrate[1/Sqrt[(1+zz)^2(1+ Ωm zz)-zz(zz+2)Ωλ],{zz,0,z}] ],Ωm+Ωλ>1}}]},{z,0,3,0.25}]])

Rozbieżmy definicję na części składowe:

(* deklaracja funkcji tylko od parametrów dofitowanych/szukanych *)
dL[Ωm_?NumberQ,Ωλ_?NumberQ,H0_?NumberQ]:= 
(
(* ta część odpowiada za zapamiętywanie raz obliczonych całek *)
dL[Ωm,Ωλ,H0]=
(* zwracaną wartością musi być pewna czysta funkcja; wynikiem działania polecenia Interpolation jest właśnie taki obiekt *)
Interpolation[
(* generowanie tabeli wartości funkcji *)
Table[{z,
(* Piecewise pozwala zapisać trzy wzory, działające w różnych zakresach zmiennych *)
Piecewise[
(*trzy wzory *)
{{c/H0(1+z)/Sqrt[1-Ωm-Ωλ]Sinh[Sqrt[1-Ωm-Ωλ]
(*całkowanie numeryczne*)
 NIntegrate[1/Sqrt[(1+zz)^2(1+ Ωm zz)-zz(zz+2)Ωλ],{zz,0,z}] ],Ωm+Ωλ<1},{c/H0(1+z)NIntegrate[1/Sqrt[(1+zz)^2(1+ Ωm zz)-zz(zz+2)Ωλ],{zz,0,z}] ,Ωm+Ωλ==1},{c/H0(1+z)/Sqrt[Ωm+Ωλ-1]Sin[Sqrt[Ωm+Ωλ-1] NIntegrate[1/Sqrt[(1+zz)^2(1+ Ωm zz)-zz(zz+2)Ωλ],{zz,0,z}] ],Ωm+Ωλ>1}}]},{z,0,3,0.25}]])

Powyższy przykład jest nieco skomplikowany ( ale za to realistyczny! ). Generalnie, zalecany schemat tworzenia funkcji dla FindFit jest nastepujący:

funkcja[ parametr1_?NumberQ, parametr2_?NumberQ]
:=
funkcja[parametr1, parametr2]
=
czysta funkcja jednej zmiennej

Ostatnia linijka może być wygenerowana w różny sposób:

a) Interpolation - gdy wewnątrz obliczamy jakieś punkty, np: numeryczne wartości całki, rozwiązanie równania, itp.)
b) y/.NDSolve[ rownanie i war pocz. , y, x] - automatycznie tworzy f. interpolującą
itp. itd.

Wywołanie tak napisanej funkcji składa się z 2 faz:

1) dL[OmegaM,OmegaΛ,Hubble] - wynikiem jest ,,czysta’’ funkcja
2) dL[OmegaM,OmegaΛ,Hubble][z] - obliczamy wartość czystej funkcji dL[OmegaM,OmegaΛ,Hubble] dla argumentu z

Pozostała część analizy to już sama przyjemność. Niestety, w przypadku gdy mierzymy się  z problemem nowym, funkcję podobną do powyższej musimy napisać sami. Jej końcowa postać zaprezentowana powyżej  nie jest oczywista. W rzeczywistości miała ona numer 3  i była jedną z pięciu przetestowanych. Pozostałe także działały, ale wolniej lub gorzej.

model1=NonlinearModelFit[zdL,dL[OmegaM,OmegaΛ,Hubble][z],{{OmegaM,0.3},{OmegaΛ,0.8},{Hubble,50}},z]  ;

model1["ParameterTable"]

model2=NonlinearModelFit[zdL,dL3[OmegaM,OmegaΛ,68][z],{{OmegaM,0.3},{OmegaΛ,0.8}},z]  ;

model2["ParameterTable"]

model3=NonlinearModelFit[zdL,dL3[1-OmegaΛ,OmegaΛ,72][z],{{OmegaΛ,0.8}},z] ;

model3["ParameterTable"]

model3["BestFitParameters"]

model4=NonlinearModelFit[zdL,dL3[0,0,H][z],{{H,50}},z] ;

model4["ParameterTable"]

m1=Plot[Tooltip[dL3[OmegaM,OmegaΛ,Hubble][z]/.model1["BestFitParameters"],model1["ParameterTable"]],{z,0,2},PlotStyle->Red];

m2=Plot[Tooltip[dL3[OmegaM,OmegaΛ,68][z]/.model2["BestFitParameters"],model2["ParameterTable"]],{z,0,2},PlotStyle->Blue];

m3=Plot[Tooltip[dL3[1-OmegaΛ,OmegaΛ,68][z]/.model3["BestFitParameters"],model3["ParameterTable"]],{z,0,2},PlotStyle->Black];

Show[ErrorListPlot[zdLerr],m1,m2,m3]

Manipulate[Show[ErrorListPlot[zdLerr],Plot[dL3[OmegaM,OmegaΛ,Hubble][z],{z,0,2}]],{Hubble,40,90,2},{{OmegaM,1.0},0,2,0.1},{{OmegaΛ,0},-2,2,0.2},ControlPlacement->{Right}]

Znaczenie wykresów w skali logarytmicznej

W pewnych sytuacjach spotykamy się z danymi, o których z punktu widzenia teoretycznego nic nie wiemy. Tradycyjną metodą szybkiej i skutecznej analizy takich danych są wykresy w których jedna lub obie osie zostały zlogarytmowane. Jeżeli na wykresie logarytmicznym dane układają się wzdłuż linii prostej, podpowiada nam to jaki wzór być może będzie pasował do danych lub ich części.

Wykres logartymiczny

funs={Tooltip[Sqrt[x],"(x^n)"] ,Tooltip[x^3,"(x^n)"],Tooltip[Exp[x],"(a^x)"],Tooltip[2 x+3,"Liniowa"],Tooltip[Log10[x],"Log"]}

Na wykresie w którym na osi pionowej mamy logarytm wartości funkcji, funkcje wykładnicze (typu (a^x) ) pojawiają się jako linie proste:

LogPlot[funs,{x,0,12},PlotStyle->{{Black},{Black,Dotted},{Black,DotDashed},{Black,Dashed}}]

Poniższy manipulator ewidentnie demonstruje, że każda funkcja wykładnicza na wykresie LogPlot pojawia się jako linia prosta:

Manipulate[ LogPlot[a b^(c x),{x,0,12}],{a,0.1,12},{b,0.1,4},{c,-4,4}]

Wykres log-log

Na wykresie podwójnie logarytmicznym, funkcje potęgowe (typu (x^n) ) są liniami prostymi:

LogLogPlot[funs,{x,0.01,12},PlotStyle->{{Black},{Black,Dotted},{Black,DotDashed},{Black,Dashed}}]

Manipulate[LogLogPlot[a x^n,{x,0.01,12}],{a,0.1,12},{n,-2,2,1/4}]

Wykres ,,półlogarytmiczny’’

LogLinearPlot[funs,{x,0.01,12},PlotStyle->{{Black},{Black,Dotted},{Black,DotDashed},{Black,Dashed}}]

Import/eksport danych

Pliki tekstowe: Import

Najważniejsza w praktyce jest umiejętność operowania plikami tekstowymi. Typowo, dane mają postać kilku kolumn liczb. Poniżej wczytujemy zawartość pewnego pliku:

dane=Import["G:\\DyskGoogle\\Dysk Google\\pair.dat"] ;

Jeżeli dane mają duży rozmiar, to należy zablokować wyświetlanie ich zawartości na ekranie średnikiem. Skrócone wyświetlanie (podgląd) umożliwia polecenie Short:

dane//Short

Inny sposób sprawdzenia czy wczytano liczby w spodziewanej formie (3 liczby na wiersz), to wyświetlenie pierwszego elementu:

dane[[1]]

lub ostatniego:

dane[[-1]]

Ilość wierszy zapisujemy do zmiennej num:

num=Length[dane]

Aby zobaczyć pełną strukturę danych używamy Dimensions:

Dimensions[dane]

Powyżej widać, że dane składają się z 301 list, a każda z nich ma 3 elementy. Powiedzmy, że chcemy narysować zależność drugiej kolumny od pierwszej. Można przepisać pojedyncze elementy:

dane1=Table[{dane[[ii,1]],dane[[ii,2]]},{ii,1,num}];

ListPlot[dane1,PlotRange->All,Joined->True]

Lub posłużyć się funkcjami operującymi na kolumnach i wierszach:

ListPlot[dane[[;;,1;;2]],PlotRange->All]

Jeszcze inny sposób:

dane1a={Transpose[dane][[1]],Transpose[dane][[2]]}//Transpose;

ListLogPlot[dane1a,PlotRange->All,Joined->True]

Mathematica usiłuje odgadnąć jak wczytać plik na podstawie jego rozszerzenia. Często prowadzi to problemów, gdy rozszerzenie jest (z punktu widzenia programu) błędne. Polecenie Import posiada drugi argument, określający jaka jest zawartość pliku. Kwalifikator zawsze ma postaś łańcucha tekstowego, dlatego podajemy go w cudzysłowach. Np: “TSV” oznacza tab separated values, czyli wartości rozdzielone znakiem tabulatora. Pełna lista jest dostępna w Help-ie.

Import["G:\\DyskGoogle\\Dysk Google\\pair.dat","TSV"] ;

Short[%]

Często dział użycie “Table”, czyli tabeli wartości rozdzielonych spacjami lub tabulatorami (tzw. Whitespace, https://pl.wikipedia.org/wiki/Znaki_niedrukowalne ), daje spodziewany efekt:

Import["G:\\DyskGoogle\\Dysk Google\\pair.dat","Table"] ;

Short[%]

Niestety, zdarza się, że Mathematica nie jest w stanie odczytać poprawnie pliku, gdyż jego struktura jest błędna lub bardzo skomplikowana. Wtedy zmuszeni jesteśmy do ,,ręcznego’’ rozbioru łancuchów tekstowych. Pojedyncze linie wczytujemy kwalifikatorem “Lines”:

linie=Import["G:\\DyskGoogle\\Dysk Google\\pair.dat","Lines"] ;

Short[%]

Poniżej widzimy, że dane mają postać liczb rozdzielonych znakiem “\t”, czyli tabulatorem (klawisz nad Caps Lock).

linie[[1]]

FullForm[%]

Tekst możemy analizować na wiele sposobów. Np: polecenie StringSplit pozwala rozbić tekst w zadanych miejscach:

StringSplit[%,"\t"]

Dane powyżej mają postać 3 osobnych pól, ale nadal jest to tekst. Aby zamienić go na liczbę możemy użyć np: ToExpression:

%[[1]]//FullForm

ToExpression[%]

Operowanie plikami i łańcuchami tekstowymi to rozległy temat. Niektóre zastosowania można znaleźć na stronie ćwiczeń z bioinformatyki https://ribes.if.uj.edu.pl/bioinfo/.

Eksportowanie wyrażeń

w=ListLogLogPlot[dane1a,Joined->True,InterpolationOrder->0]

Export["G:\\test.png",w,ImageResolution->300]

Export["G:\\wykres.m",w]

Import["G:\\wykres.m"]

Compress[w]

Export["G:\\wykres.txt",%]

Import["G:\\wykres.txt"]

Uncompress[%%]

Zintegrowane bazy danych (wymaga połączenia internetowego!)

Przykłady baz danych

Bardzo ciekawą funkcją programu Mathematica jest możliwość pobierania danych dotyczących wielu dziedzin nauki oraz naszego życia. Ilość tych baz jest spora. Dla przykładu, regularnie korzystam z:

AstronomicalData[]//Short (* dane obiektów astronomicznych *)

IsotopeData[]//Short (* dane jąder atomowych *)

GenomeData["HOXA*"] (* geny m .in. człowieka*)

ChemicalData["AminoAcids"]

WeatherData["Krakow","Temperature"]

PolyhedronData["Octahedron"]

KnotData["Trefoil"]

Przykład 1: sposób użycia bazy danych na przykładzie energii wiązania stabilnych izotopów

Podstawowy schemat pracy z bazami danych jest następujący. Po pierwsze potrzebujemy listy zawartych w bazie obiektów. Robimy to wypisując nazwę bazy danych, bez żadnych argumentów:

IsotopeData[]//Short

iso=IsotopeData[];

Kolejny krok, to wydobycie informacji o tym, co jest zapisane w bazie. Służy do tego specjalne pole “Properties” (podobnie jak np: w NonlinearModelFit ):

p=IsotopeData["Properties"]

Teraz możemy przeglądać bazę:

Manipulate[
{iso[[k]],p[[i]],
IsotopeData[iso[[k]],p[[i]]]
},{k,1,Length[iso],1},{i,1,Length[p],1}]

To co teraz zrobimy z tymi danymi zależy wyłącznie od naszej pomysłowości i umiejętności. Na przykład, odtworzymy słynny w astrofizyce i fizyce jądrowej wykres pokazujący energię wiazania stabilnych jąder atomowych:

num=Length[iso]

stabilne={}

For[k=1,k≤num,k++,
If[IsotopeData[iso[[k]],"Stable"],AppendTo[stabilne,{IsotopeData[iso[[k]],"AtomicNumber"],IsotopeData[iso[[k]],"NeutronNumber"],IsotopeData[iso[[k]],"BindingEnergy"]}]
]
]

BE=ListPointPlot3D[stabilne,Filling->Bottom,PlotRange->{7.4,8.8},AxesLabel->{"Z","N","E [MeV]"},
Ticks->{Automatic,Automatic,{7.5,8,8.5,8.8}},BoxStyle->Thick]

Przykład 2: rozkład liczby miast w zależności od liczby mieszkańców w Polsce i Niemczech

Dla przykładu, spróbujemy odpowiedzieć na pytanie: jak wygląda rozkład liczby miast w zależności od liczby mieszkańców w Polsce. Czy rozkład ten jest taki sam jak np: w Niemczech?

Quit[]

miasta=CityData[]

miastaPL=Select[miasta,#[[3]]=="Poland"&];

numPL=Length[miastaPL]

CityData["Properties"]

Sort[miastaPL,CityData[#1,"Population"]>CityData[#2,"Population"]&]//Short

popPL=Table[CityData[miastaPL[[i]],"Population"],{i,1,numPL-1}];Short[popPL]

Total[popPL]

?Count

f=Function[{population,s},Count[population,q_/;q>s]]

wykresPL=LogLogPlot[f[popPL,s],{s,10^4,10^6}]

miastaDE=Select[miasta,#[[3]]=="Germany"&]

popDE=Table[CityData[miastaDE[[i]],"Population"],{i,1,Length[miastaDE]}];

wykresDE=LogLogPlot[f[popDE,s],{s,10^4,10^6}]

w1=Show[wykresDE,wykresPL]

datPL=Table[{10^lgs,f[popPL,10^lgs]},{lgs,4,6,1/3}]//N

fitPL=FindFit[datPL, a x^p,{a,p},x]

datDE=Table[{10^lgs,f[popDE,10^lgs]},{lgs,4,6,1/3}]//N

fitDE=FindFit[datDE, a x^p,{a,p},x]

w2=LogLogPlot[{a x^p/.fitDE,a x^p/.fitPL},{x,10^4,10^6}]

Show[w1,w2]

w3=ListLogLogPlot[datPL]

Show[w1,w2,w3]

NonlinearModelFit[datPL,a x^p,{a,p},x]["ParameterTable"]

datPL2=Table[{lgs,Log10[f[popPL,10^lgs]]},{lgs,4,6,1/3}]//N

LinearModelFit[datPL2,x,x][1]

Show[ListPlot[datPL2],Plot[LinearModelFit[datPL2,x,x][X],{X,4,6}]]

LinearModelFit[datPL2,x,x]["ParameterTable"]