opgaver:Uge1

Disse statistik-øvelser er redigerede udgaver af opgaverne 3.4, 3.25, 4.6, 4.28, 5.11 og 5.31 fra R. J. Taylor's An Introduction to Error Analysis.

Det er meningen af så meget som muligt af opgaverne skal laves i MATLAB.

Opgave 1 - Tælletal
Når en prøve med radioaktive atomer henfalder, vil antallet af radioaktive atomer falde, og prøvens radioaktivitet vil falde proportionalt med dette. For at undersøge denne effekt måler en fysiker på partiklerne der udsendes fra en radioaktiv prøve i løbet af to timer. Hun tæller antallet af partikler der udsendes i løbet af en 1 minut lang periode med halve timers intervaller, med følgende resultater:


 * {| cellpadding="4"


 * Tid gået, $t$ (timer): || 0.0 || 0.5 || 1.0 || 1.5 || 2.0
 * Antal tællinger, $\nu$, på 1 minut: || 214 || 134 || 101 || 61 || 54
 * }
 * }
 * }

Spørgsmål a
Brug MATLAB til at plotte antallet af tællinger imod den tid der er gået. Inkludér errorbars for at vise usikkerheden på tallene. (Negligér usikkerheder i den tid der er gået).

For en hjælp til hvordan man plotter data, se "MATLAB® Primer" (tidligere kaldet "Getting Started with MATLAB") guiden, der ligger på Absalon. Kapitel 4 indeholder en masse plotning generelt. Prøv at slå funktionen  op i MATLAB hjælpen bagefter.

For at plotte datapunkterne med usikkerheder skal data først indlæses. Usikkerhederne er kvadratroden af tælletallet, så



Herefter kan data plottes vha. MATLAB funktionen, der plotter $(x,y)$ data, hvor der er usikkerheden $e$ på $y$-dataen. En figur åbnes med



der plottes med



og til sidst laves der en



for at blive ved med at plotte oven i det allerede lavede (ellers vil de næste kommandoer overskrive errorbar plottet). Jeg sætter nu også akser, titel og akselabels med



Spørgsmål b
En teori forudsiger at antallet af udsendte partikler burde falde eksponentielt ved $\nu = \nu_0 \exp{(-rt)}$, hvor der (i dette tilfælde) gælder at $\nu_0 = 200$ og $r = 0.693 \,\text{hr}^{-1}$. Plot denne forventede kurve på din graf oven på din tidligere graf, og kommentér på hvor godt data ser ud til at passe med den teoretiske forudsigelse.

For at plotte en funktion skal man først definere sin førsteakse som en vektor, f.eks. . Man kan også bruge  funktionen til at lave en liste fra $x$ til $y$ med $n$ punkter. Antallet af punkter definerer opløsningen på det plot man vil lave.

Herefter kan man udregne sine værdier til andenaksen som sin funktion af førsteakseværdierne, f.eks. .

Man definerer linjen for teorien rent numerisk ved at lave en tids-akse som man selv bestemmer hvor detaljeret er, f.eks.



og så kan man udregne den teoretiske linje ved



Hvis man huskede sin  tidligere burde man nu bare kunne plotte oven i sit errorbar plot med



Som en sidste krone på værket kan man tilføje en legend med



Prøv at lave plottet igen med forskellige layouts af punkterne og linjen. F.eks. er der de forskellige farver  (rød),   (grøn),   (blå),   (gul),   (magenta) og   (sort); punkttyperne   (lille prik),   (krydser),   (store prikker) og   (stjerner); og linjetyperne   (solid),   (prikker) og   (stiplet), og selvfølgelig mange flere. Prøv at finde disse options i MATLAB hjælpen under.

Du kan også lege med egenskaber som,  ,   og   for at gøre linjen og punkterne mere tydelige.

Bonus
I uge 3 skal I lære hvordan man fitter en funktion til datapunkter i MATLAB - så i stedet for at plotte en given kurve vil I skulle plotte den bedst fittede kurve til data.

Lad os sige at data allerede er indlæst i variablene  og   (med usikkerhed  ), som ovenfor. Disse er også allerede plottet som det blev gjort i spørgsmål a.

For at kunne fitte en linje skal man selvfølgelig først definere sin fitte-kurve, der i dette tilfælde vil være $\nu = \nu_0 \exp{(-rt)}$, hvor $\nu_0$ og $r$ nu er fitte-parametre. En funktion kan defineres på mange måder i MATLAB, men når man fitter med simple funktionsudtryk i MATLAB er det nemmest simpelthen bare at fodre fitte-rutinen med den direkte.



Bemærk at  betyder den transponerede af. Funktionen  kræver nemlig søjlevektorer som data, ikke rækkevektorer, som både   og   er defineret som. Objektet  indeholder nu information om fittet, og man kan plotte kurven med de bedste værdier af   og   med



man kan også få parametrene ud ved at kalde variablene  og. Det er også muligt at gå dybere ned i  objektet, men dette skal ikke uddybes her.

Bemærk også at ovenstående ikke fodrer MATLAB med et start-gæt på de to fitte-parametre. Dette bør man normalt gøre, da programmet ellers vælger tilfældige startpunkter, og her ikke altid vil konvergere. Man kan gøre det ved at bruge  option'en, og dette vil vi komme ind på senere.

Opgave 2 - Usikkerheden af et potensudtryk
Fra Barlow's regel (4.10) om usikkerheder af funktionsudtryk, ved vi at et $q = x^2$ har en relativ usikkerhed der er dobbelt så stor som den relative usikkerhed i $x$;


 * $ \dfrac{\sigma_q}{q} = 2 \dfrac{\sigma_x}{x} . $

Overvej nu det følgende (forkerte) argument: Vi kan tænke på $x^2$ som $x$ gange $x$, så


 * $ q = x \times x; $

derfor vil der ifølge Barlow's regel (4.14) gælde at


 * $ \dfrac{\sigma_q}{q} = \sqrt{ \left( \dfrac{\sigma_x}{x} \right)^2 + \left( \dfrac{\sigma_x}{x} \right)^2 }

= \sqrt2 \dfrac{\sigma_x}{x}. $

Denne konklusion er forkert. Forklar hvorfor.

Bonus: Prøv at undersøge dette i MATLAB ved at sammenligne den relative usikkerhed (SD divideret med gennemsnitsværdien) for en række af tilfældige tal med den relative usikkerhed på samme række af tal kvadreret. Du kan evt. læse i MATLAB-hjælpen hvad funktionen  gør.

Se på hvad antagelserne for at formel (4.14) gælder faktisk er - passer dette med hvad der arbejdes med i denne opgave?

For at se på det i MATLAB kan man generere en række tilfældige tal vha. funktionen. Slå den op i MATLAB hjælpen og se hvad den kan.

Konklusionen i den fejlagtige udledning er forkert fordi man har brugt udtrykket i formel (4.14) uden at overholde dennes antagelser - nemlig at de variable man ganger sammen er uafhængige. Og nu er $x$ jo fuldkommen afhængig af sig selv - hvorfor konklusionen altså kun kan være forkert.

For at se på dette i MATLAB kan man starte med at lave en liste med normalfordelte tal



Bemærk at der lægges 10 til, da vi gerne vil dividere med gennemsnittet senere (og at dividere med nul er jo ikke så godt!).

Man kan nu helt manuelt se på den relative usikkerhed på  og på   i anden potens, ved at bruge de indbyggede MATLAB funktioner for standardafvigelse og gennemsnit,   og  :



Man finder at nummer to usikkerhed ovenfor ligner at man ganger den første med 2, ikke med $\sqrt{2}$, ved at sammenligne



Og dermed er det altså også rent numerisk vist at den naive fremgangsmåde ikke er sand.

Opgave 3 - Flere tælletal
I Barlow kap. 3.3.1 kan I læse at i et tælleeksperiment er usikkerheden på en tælling givet af "kvadratrods-reglen" til at være kvadratroden af tallet. Denne regel kan gøres mere præcis med følgende udsagn: Hvis man tager adskillige tællinger


 * $ \nu_1, \nu_2, \ldots, \nu_N $

af antallet $\nu$ af tilfældige hændelser der sker inden for tidsrummet $T$, så gælder der: (1) det bedste estimat for det sande gennemsnit af hændelser inden for tidsrummet $T$ er gennemsnittet $\overline\nu = \sum \nu_i/N$ af målingerne, og (2) standardafvigelsen (også kaldt SD, for Standard Deviation) af de observerede tal er approksimativt lig med kvadratroden af dette samme bedste estimat; altså, usikkerheden i hver måling er $\sqrt{\overline\nu}$. I specialtilfældet hvor der kun måles én gang, hvor man får et antal tællinger $\nu$, er det bedste estimat $\nu$, mens usikkerheden er kvadratroden $\sqrt\nu$.

En fysiker bruger en Geiger-tæller til at måle antallet af kosmiske partikler der kommer til hans laboratorie i givne to-sekunders intervaller. Han tæller dette tal 20 gange med følgende resultater:


 * {| cellpadding="2"


 * 10, || 13, || 8, || 15, || 8, || 13, || 14, || 13, || 19, || 8,
 * 13, || 13, || 7, || 8, || 6, || 8, || 11, || 12, || 8,  || 7.
 * }
 * }
 * }

Spørgsmål a
Find gennemsnittet og standardafvigelsen (SD) på disse tal vha. MATLAB.

Først defineres datasættet som en enkel liste med



Derefter kan gennemsnittet og SD udregnes enten manuelt eller ved brug af  funktionen. F.eks.



Dermed har man gennemsnittet og datasættets faktiske standardafvigelse.

Husk ikke at forvirre standardafvigelsen SD med standardafvigelsen på gennemsnittet (også kaldt SDOM, for Standard Deviation Of the Mean). Hvis man her vil finde SDOM, skal man dividere med kvadratroden af antallet af punkter,



Spørgsmål b
Standardafvigelsen på tallene skulle gerne være ca. lig med kvadratroden på deres gennemsnit. Hvor godt passer dette?

Kvadratroden af gennemsnittet udregnes simpelthen med (ved brug af løsningen fra spørgsmål a)



Denne standardafvigelse er faktisk et over-estimat, men den passer ok med den udregnede  fra ovenfor.

Opgave 4 - Systematiske fejl
Systematiske fejl kommer nogle gange fra at fysikeren helt uvidende måler det forkerte. Her er et eksempel: En studerende forsøger af måle $g$ ved at bruge et pendul lavet af en stålkugle der hænger i en let snor (se figuren). Han måler fem forskellige længder af pendulet og de tilhørende perioder $T$ som følger:


 * {| cellpadding="4"


 * Længde, $l$ (cm): || 51.2 || 59.7 ||  68.2 ||  79.7 ||  88.3
 * Periode, $T$ (s): || 1.448 || 1.566 || 1.669 || 1.804 || 1.896.
 * }
 * }
 * }

Spørgsmål a
For hver datapar udregner han $g$ som $g = 4 \pi^2 l/T^2$. Han udregner derefter gennemsnittet af de fem værdier, deres SD, og deres SDOM. Ved at antage at alle hans fejl er tilfældige kan han bruge SDOM som sin endelige usikkerhed, og skriver sit resultat på standardformen gennemsnit $\pm$ SDOM. Find hans resultat for $g$ vha. MATLAB.

Start med at opskrive data i to 1D arrays - et array for  og et array for. Nu kan du udregne  vha. den givne funktion, og lave statistik på denne.

Dataen læses ind med



og tyngdeaccelerationskoefficienten udregnes ved



Herefter kan man finde gennemsnittet, SD og SDOM af værdierne vha. de indbyggede MATLAB funktioner



Spørgsmål b
Han sammenligner nu sit resultat med den accepterede værdi $g = 979.6 \,\text{cm}/\text{s}^2$ og forfærdes over at se at hans diskrepans (afvigelse fra den accepterede værdi) er tæt på 10 gange større end hans usikkerhed. Bekræft denne sørgelige konklusion.

Diskrebansen findes bare som forskellen,



Spørgsmål c
Efter at have gennemtjekket alle sine udregninger konkluderer han at han må have overset en systematisk fejl. Han er helt sikker på at der ikke var problemer med målingen af perioden $T$, så han spørger sig selv: Hvor stor skulle en systematisk fejl i længden $l$ være for at grænserne for den totale usikkerhed lige præcis inkluderer den accepterede værdi $979.6 \,\text{cm}/\text{s}^2$? Vis med MATLAB at svaret er ca. 1 cm.

Prøv at lave en udregning hvor du ændrer din værdi for $l$ en lille smule, og se hvor langt du så skal ændre den for at få en $g$ der passer med den korrekte $g$. Det kan evt. være en god idé at se lidt på -løkker i MATLAB, så programmet kan iterere for dig.

For at finde hvor stor fejlen skulle være for at den rigtige værdi er inden for SDOM, kan man lave et loop der lægger lidt ekstra til $l$, og så udregner forskellen mellem den fundne $g$ og den rigtige,



Den nemmeste måde at se resultatet på er ved at plotte det



Man kan se at den skærer 0 ca. omkring 1 cm (se figuren), hvilket netop var hvad der skulle vises.

Man kunne selvfølgelig også sammenligne den teoretiske værdi for $g$ med gennemsnits-værdien for sin udregnede $g$, hvilket selvfølgelig bare skal gøres ved at erstatte  udregningen ovenfor med



Spørgsmål d
Dette resultat ville betyde at hans længdemålinger har en systematisk fejl på omkring en centimeter - en konklusion som han først afviser som absurd. Mens han stirrer på pendulet kommer han dog i tanke om at 1 cm er ca. radius af kuglen, og at de længder han har målt var længderne af snoren. Da den korrekte længde af et pendul er afstanden fra fastgørelsespunktet til centrum af kuglen (se figuren), er hans målinger altså systematisk blevet målt forkert med radius af kuglen. Han bruger derfor en skydelære til at måle kuglens diameter, der viser sig at være 2.00 cm. Lav de nødvendige korrektioner til hans data og udregn hans endelige resultat for $g$ med dennes usikkerhed.

For at rette den systematiske fejl kan man bare lægge den ene cm til sine $l$-værdier,



og så udregner man sin $g$-værdi forfra, ligesom det blev gjort før:



Opgave 5 - Gauss-fordelingen
Plot Gauss-fordelingen (også kaldt normalfordelingen),


 * $ G_{X,\sigma}(x) = \dfrac{1}{\sigma\sqrt{2\pi}} \exp{\left( -\dfrac{(x-X)^2}{2\sigma^2} \right)}, $

i den samme MATLAB figur for de følgende to sæt af parametre: $X=2$, $\sigma=1$, og $X=3$, $\sigma=0.3$. Tilføj en titel, labels på akserne, og en legend til dit plot, og sammenlign de to grafer.

Først definerer man sin $x$-akse, f.eks. ved



Derefter giver man sine parametre



til sidst kan man så udregne sine funktionsværdier ved



Man kunne også definere en funktion for sin Gauss-fordeling ved definitionen



og med denne kan man bare skrive



Til sidst plotter man de to funktioner sammen. Her kan man bruge at  godt kan finde ud af at plotte flere funktioner af gangen:



Opgave 6 - Binning
I den følgende liste er der 40 målinger $t_1,\ldots,t_{40}$ af tiden det tager en sten at falde fra et vindue til jorden (alle i hundrededele sekunder):


 * {| cellpadding="2"


 * 63|| 58|| 74|| 78|| 70|| 74|| 75|| 82|| 68|| 69
 * 76|| 62|| 72|| 88|| 65|| 81|| 79|| 77|| 66|| 76
 * 86|| 72|| 79|| 77|| 60|| 70|| 65|| 69|| 73|| 77
 * 72|| 79|| 65|| 66|| 70|| 74|| 84|| 76|| 80|| 69
 * }
 * 72|| 79|| 65|| 66|| 70|| 74|| 84|| 76|| 80|| 69
 * }
 * }
 * }

Spørgsmål a
Brug MATLAB til at udregne standardafvigelsen $\sigma_t$ for alle 40 målinger.

Først og fremmest skal datasættet lige indlæses,



Standardafvigelsen findes så simpelt med den indbyggede MATLAB funktion,



Spørgsmål b
Udregn gennemsnittene $\overline t_1,\ldots,\overline t_{10}$ af fire målinger af gangen (f.eks. af de fire målinger i hver af de 10 kolonner, eller af fire nabo-tal ad gangen). Du kan tænke på dataen som om den kom fra 10 eksperimenter, hvor man i hvert eksperiment fandt gennemsnittet af fire tidsmålinger. Givet resultatet i del (a), hvad forventer du så at standardafvigelsen på de 10 gennemsnit $\overline t_1,\ldots,\overline t_{10}$ er? Udregn den.

For at binne med fire målinger af gangen er det nemmest at arbejde med et loop der gør det automatisk for dig. For at bruge et generelt antal punkter at binne med defineres tallet først:



Herefter startes der to tomme arrays op - et til gennemsnittene, og et til standardafvigelserne ( laver et $n\times m$ array der indeholder nuller i alle elementer):



Et loop laver udregningerne:



Man kan så endelig se  og   ved simpelthen at printe dem ud via kommandolinjen.

Ud fra del (a) ville man forvente at SD er fordelt omkring den SD man fandt tidligere, og det er også det der sker.

Spørgsmål c
Plot histogrammer for både de 40 individuelle målinger $t_1,\ldots,t_{10}$ og de 10 gennemsnit $\overline t_1,\ldots,\overline t_{10}$. Brug de samme skalaer og bin størrelser for begge plot, så de kan sammenlignes. Hvad forventer du af de to histogrammer, og ser de faktisk sådan ud?

For at plotte histogrammer der kan sammenlignes skal man bruge histc, hvor man kan vælge bin-grænser manuelt. Histogrammerne laves ved først at definere grænserne for bins'ene ved



For at kunne plotte bins'ene med en korrekt førsteakse skal man også finde centrum af hvert bin. Dette gøres med



Herefter udregnes histogrammerne med



og de kan nu plottes med  funktionen,



Man kan tydeligt se at hist10 er mere centreret, som forventet.

Samlede løsninger

 * Et samlet dokument med MATLAB kode til at løse alle opgaverne ovenfor kan hentes her: [[Media:opgaver_uge1_2_taylor.m|opgaver_uge1_2_taylor.m]]