opgaver:Uge1
Pia Jensen (Talk | contribs) |
Pia Jensen (Talk | contribs) |
||
| Line 329: | Line 329: | ||
:<code>sd = zeros(1,10);</code> | :<code>sd = zeros(1,10);</code> | ||
| − | + | Et loop laver udregningerne: | |
| − | + | :<code>for n = 0:9</code> | |
| − | + | ||
| − | :<code>for n = | + | |
::<code> gn(n+1) = sum(t(n*bin+1:n*bin + bin))/bin;</code> | ::<code> gn(n+1) = sum(t(n*bin+1:n*bin + bin))/bin;</code> | ||
::<code> sd(n+1) = std(t(n*bin+1:n*bin + bin));</code> | ::<code> sd(n+1) = std(t(n*bin+1:n*bin + bin));</code> | ||
Revision as of 11:01, 10 April 2012
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.
Contents |
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:
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 Getting Started with MATLAB guiden, der ligger på Absalon. Kapitel 4 indeholder en masse plotning generelt. Prøv at slå funktionen errorbar 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å
t = 0:0.5:2;nu = [214 134 101 61 54];nu_usikk = sqrt(nu);
Herefter kan data plottes vha. MATLAB funktionen errorbar(x,y,e), der plotter $(x,y)$ data, hvor der er usikkerheden $e$ på $y$-dataen. En figur åbnes med
figure
der plottes med
errorbar(t,nu,nu_usikk,'ro')
og til sidst laves der en
hold on
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
axis([-0.5 2.5 40 240])title('Opgave 1 - Radioaktivt henfald')xlabel('Tid i timer')ylabel('Tælletal')
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. x = 0:0.01:1. Man kan også bruge linspace(x,y,n) 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. y = sin(x).
Man definerer linien for teorien rent numerisk ved at lave en tids-akse som man selv bestemmer hvor detaljeret er, f.eks.
tt = -0.25:0.01:2.25;
og så kan man udregne den teoretiske linie ved
nu_0 = 200;r = 0.693;nu_exp = nu_0 * exp(-r*tt);
Hvis man huskede sin hold on tidligere burde man nu bare kunne plotte oven i sit errorbar plot med
plot(tt,nu_exp,'b:')
Som en sidste krone på værket kan man tilføje en legend med
legend('Datapunkter','Teori-linie')
Prøv at lave plottet igen med forskellige layouts af punkterne og linien. F.eks. er der de forskellige farver r (rød), g (grøn), b (blå), y (gul), m (magenta) og k (sort); punkttyperne . (lille prik), x (krydser), o (store prikker) og * (stjerner); og linietyperne - (solid), : (prikker) og -- (stiplet), og selvfølgelig mange flere. Prøv at finde disse options i MATLAB hjælpen under LineSpec.
Du kan også lege med egenskaber som MarkerSize, LineWidth, MarkerEdgeColor og MarkerFaceColor for at gøre linien og punkterne mere tydelige.
Opgave 2 - Usikkerheden af et potensudtryk
Fra Taylor's regel (3.10) om usikkerheder af potensudtryk, ved vi at et $q = x^2$ har en relativ usikkerhed der er dobbelt så stor som den relative usikkerhed i $x$;
- $ \dfrac{\delta q}{q} = 2 \dfrac{\delta 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 Taylor's regel (3.18) gælde at
- $ \dfrac{\delta q}{q} = \sqrt{ \left( \dfrac{\delta x}{x} \right)^2 + \left( \dfrac{\delta x}{x} \right)^2 } = \sqrt2 \dfrac{\delta 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 randn gør.
Se på hvad antagelserne for at formel (3.18) 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. randn 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 (3.18) 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
x = 10 + randn(1,500);
Man kan nu helt manuelt se på den relative usikkerhed på x og på x i anden potens, ved at bruge de indbyggede MATLAB funktioner for standardafvigelse og gennemsnit, std og mean:
std(x)/mean(x)std(x.^2)/mean(x.^2)
Man finder at nummer to usikkerhed ovenfor ligner at man ganger den første med 2, ikke med $\sqrt{2}$, ved at sammenligne
2*std(x)/mean(x)sqrt(2)*std(x)/mean(x)
Og dermed er det altså også rent numerisk vist at den naive fremgangsmåde ikke er sand.
Opgave 3 - Flere tælletal
I Taylor kap. 3 lærte I 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 (bevises i Taylor kap. 11): 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 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$; dette resultat er bare kvadratrods-reglen fra Taylor kap. 3, med den ekstra information at "usikkerheden" faktisk er standardafvigelsen, der definerer grænserne hvorimellem man kan være cirka 68 % sikker på at det sande svar ligger.
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:
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 på disse tal vha. MATLAB.
Først defineres datasættet som en enkel liste med
nu = [10 13 8 15 8 13 14 13 19 8 13 13 7 8 6 8 11 12 8 7];
Derefter kan gennemsnittet og SD udregnes enten manuelt eller ved brug af mean funktionen. F.eks.
gn = sum(nu)/length(nu)std(nu)
Dermed har man gennemsnittet og datasættets faktiske standardafvigelse.
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)
sqrt(gn)
Denne standardafvigelse er faktisk et over-estimat, men den passer ok med den udregnede std(nu) 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:
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 l og et array for T. Nu kan du udregne g vha. den givne funktion, og lave statistik på denne.
Dataen læses ind med
l = [51.2 59.7 68.2 79.7 88.3];T = [1.448 1.566 1.669 1.804 1.896];
og tyngdeaccelerationskoefficienten udregnes ved
g = 4*pi^2*l./T.^2
Herefter kan man finde gennemsnittet, SD og SDOM af værdierne vha. de indbyggede MATLAB funktioner
g_mean = mean(g)sigma = std(g)sigma_mean = std(g)/sqrt(length(g))
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 er tæt på 10 gange større end hans usikkerhed. Bekræft denne sørgelige konklusion.
Diskrebansen findes bare som forskellen,
g_right = 979.6;discrepancy = g_mean - g_right
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å for-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,
x = 0.1:0.1:1.2;for i=1:length(x)g_new = 4*pi^2*(l+x(i))./T.^2;discr_new(i) = mean(g_new) + std(g_new)/sqrt(length(g_new)) - g_right;
end
Den nemmeste måde at se resultatet på er ved at plotte det
plot(x,discr_new,'.b')
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 discr_new udregningen ovenfor med
discr_new(i) = mean(g_new) - g_right;
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,
l = l + 1
og så udregner man sin $g$-værdi forfra, ligesom det blev gjort før:
g = 4*pi^2*l./T.^2;g_mean = mean(g)sigma_mean = std(g)/sqrt(length(g))
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
x = -2:0.05:6;
Derefter giver man sine parametre
X1 = 2;sigma1 = 1;X2 = 3;sigma2 = 0.3;
til sidst kan man så udregne sine funktionsværdier ved
y1 = 1/( sigma1*sqrt(2*pi) ) * exp(-(x-X1).^2/(2*sigma1^2));y2 = 1/( sigma2*sqrt(2*pi) ) * exp(-(x-X2).^2/(2*sigma2^2));
Man kunne også definere en funktion for sin Gauss-fordeling ved definitionen
MinGauss = @(x,a,b) 1/( b*sqrt(2*pi) ) * exp(-(x-a).^2/(2*b^2));
og med denne kan man bare skrive
y1 = MinGauss(x,X1,sigma1);y2 = MinGauss(x,X2,sigma2);
Til sidst plotter man de to funktioner sammen. Her kan man bruge at plot godt kan finde ud af at plotte flere funktioner af gangen:
plot(x,y1,'r-',x,y2,'b-')
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):
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
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,
t = [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];
Standardafvigelsen findes så simpelt med den indbyggede MATLAB funktion,
std(t)
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:
bin = 4;
Herefter startes der to tomme arrays op - et til gennemsnittene, og et til standardafvigelserne (zeros(n,m) laver et $n\times m$ array der indeholder nuller i alle elementer):
gn = zeros(1,10);sd = zeros(1,10);
Et loop laver udregningerne:
for n = 0:9gn(n+1) = sum(t(n*bin+1:n*bin + bin))/bin;sd(n+1) = std(t(n*bin+1:n*bin + bin));
end
Man kan så endelig se gn og sd ved simpelthen at printe dem ud via kommandolinien.
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 med
edges = 58:6:58+5*6;hist40 = histc(t,edges);hist10 = histc(gn,edges);
og de kan nu plottes med bar funktionen,
x = 1:6;figurebar(x,hist40,0.5,'y')hold onbar(x,hist10,0.3,'r')
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: opgaver_uge1_2_taylor.m
