opgaver:Uge3

Start med at udregne selve den teoretiske linie y_model, ved at definere A, B og din x-akse. Definér derefter en sigma du vil bruge som din usikkerhed. Så kan din $y$ udregnes som y_model plus en vektor der består af Gaussisk fordelte tal med spredning sigma og centrum nul.

Start med at definere de to konstanter, din x-akse og den usikkerhed du gerne vil have på punkterne, her f.eks.

A_model = 0;

B_model = 1;

sigma = 0.5;

x = -2:0.1:2;

Nu kan du udregne modellens y-værdi direkte ved

y_model = A_model + B_model*x;

Du kan nu lave en model med støj på punkterne med spredning givet ved sigma du definerede tidligere, ved at bruge randn() funktionen

y = y_model + sigma*randn(size(y_model));

Nu kan du plotte din model sammen med dine "målepunkter" (bemærk at for at få errorbars på er man nødt til at lave en vektor med samme længde som $x$ og $y$, som indeholder $\sigma$ for hver $y$-værdi - men da $\sigma$ er den samme for alle punkterne kan man bare bruge ones() funktionen),

Modellen vist som en linie, sammen med genererede data med støj.

figure

errorbar(x,y,sigma*ones(size(y_model)),'bo','MarkerFaceColor', [0.4 0.4 0.8],...

'MarkerEdgeColor',[0.4 0.4 0.8])

hold on

plot(x,y_model,'r-','LineWidth',2)

xlabel('Datapunkter x')

ylabel('Datapunkter y')

En sådan figur er vist her til højre.

Modellen og datapunkterne med støj, nu med en fit-linie (grøn).

Regressionsparametrene kan udregnes med

N = length(x);

Sy = sum(y);

Sx = sum(x);

Sxx = sum(x.^2);

Sxy = sum(x.*y);

Delta = N*Sxx-Sx^2;

A_fit = ( Sxx*Sy - Sx*Sxy )/Delta;

B_fit = ( N*Sxy - Sx*Sy )/Delta;

og fit-linien er så fundet til at være

y_fit = A_fit + B_fit*x;

Denne linie plottes oven på model-linien og punkterne med

plot(x,y_fit,'g-','LineWidth',2)

og så får man figuren til højre. Den reducerede $\chi^2$ findes til sidst med koden

chi2 = sum((y-y_fit).^2 ./ (sigma*ones(size(y)).^2);

chi2red = chi2/(N-2);

Først laves en ny serie y-modelværdier:

C_model = 0.3;

y_model_ny = A_model + B_model*x + C_model*x.^2;

Dernæst laves nye værdier med støj:

y_ny = y_model_ny + sigma_y.*randn(size(y_model));

Vi kan derefter plotte de nye værdier sammen med de gamle samt både den gamle og nye model samt fittet til de lineære data. Du kan evt. supplere med et lineært fit til de nye data.

figure

errorbar(x,y,sigma_y,'bo','MarkerFaceColor', [0.4 0.4 0.8],...

'MarkerEdgeColor',[0.4 0.4 0.8])

hold on

errorbar(x,y_ny,sigma_y,'bo','MarkerFaceColor', [0.4 0.4 0.4],...

'MarkerEdgeColor',[0.4 0.4 0.4])

plot(x,y_model,'r-','LineWidth',2)

plot(x,y_model_ny,'k-','LineWidth',2)

y_fit = A_fit + B_fit*x;

plot(x,y_fit,'g-','LineWidth',2)

xlabel('Datapunkter x')

ylabel('Datapunkter y')

hold off

Herefter udregnes det nye reducerede $\chi^2$-værdi.

chi2_ny = sum((y_ny-y_fit).^2 ./ sigma_y.^2);

chi2red_ny = chi2_ny/(N-2);

Til sidst plottes residuals.

figure

bar(x,y_ny-y_fit)

hold on

bar(x,y-y_fit,'r','BarWidth',0.5)

Start med at bruge den kode du lavede til spørgsmål a, og lav så et loop uden om det hele. I slutningen af loopet skal du gemme A, B og din reducerede $\chi^2$ i et array - så du kan plotte alle værdierne bagefter.

Man kan genbruge al sin kode fra spørgsmål a til at lave sit loop. Definitionen af A_model, B_model, sigma, x og y_model får lov til at blive, og så skal man forberede sit loop med antal gange man vil køre det, og tre lister til at gemme de resultater man får,

N_iter = 100;

A_liste = zeros(1,N_iter);

B_liste = zeros(1,N_iter);

chi2_liste = zeros(1,N_iter);

Herefter starter man sit loop lige før udregningen af y, således at hvert loop laver en ny række punkter med anderledes støj.

for i = 1:N_iter

y = y_model + sigma*randn(size(y_model));

...

Udregningen af regressionsparametrene er helt som før. Til sidst skal man gemme sine resultater i sine arrays,

...

A_liste(i) = A_fit;

B_liste(i) = B_fit;

chi2_liste(i) = chi2red;

end

Disse arrays kan man nu se lidt på. F.eks. kan man lave histogrammer eller korrelationsplots med hist eller scatter.

Modellen og datapunkterne med støj, plottet med fit-linien (grøn).

Den automatiske fitning i MATLAB kan bruges med koden

[res good] = fit(x',y','poly1')

her er poly1 den indbyggede lineære funktion, men man kan også selv skrive sin funktion ind (se på MATLAB hjælp siden for fit). Fittet kan nu plottes sammen med datapunkterne og model-linien med

figure

errorbar(x,y,sigma*ones(size(y),'ro')

hold on

plot(x,y_model,'b-')

plot(x,res(x),'g-')

xlabel('Datapunkter x')

ylabel('Datapunkter y')

Denne figur kan ses til højre.