Mercurial > hg > ltpda

comparison m-toolbox/test/straight_line_fit/straight_line_fit.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % Make some data that is a straight-line then fit with curvefit and
2 % lscov.
3 %
4 %
5 mc
6
7 m = 1.3;
8 c = 4;
9 x = 1:30;
10 ca = ao(plist('dtype', 'cdata', 'yvals', c.*ones(size(x))));
11 xa = ao(plist('dtype', 'cdata', 'yvals', x));
12
13 n = 1*randn(size(x));
14 sn = 2;
15 na = sn.*ao(plist('dtype', 'cdata', 'yvals', n));
16 ya = m.*xa + ca + na;
17
18 % Fit with straightLineFit
19
20 fy = straightLineFit(ya);
21 sigma = find(fy.procinfo, 'Sigma');
22 params = find(sigma.procinfo, 'P')
23 Pstd = find(sigma.procinfo, 'STD')
24 Pcov = find(sigma.procinfo, 'cov')
25
26
27 % theoretical errors
28 vn = var(na.y);
29 N = length(x);
30 delta = N.*sum(x.^2) - sum(x).^2;
31
32 em = sqrt(N.*vn / delta)
33 ec = sqrt(vn .* sum(x.^2./delta))
34
35 params
36 Pstd = 100.*find(sigma.procinfo, 'STD')./params
37
38
39 %% curvefit needs an xy data
40 xy = ya.convert('to xydata');
41 pl = plist('Function', 'P(1).*Xdata + P(2)', 'P0', [1 0]);
42 b = curvefit(xy, pl);
43
44
45 %% With weights of 1/sigma^2
46
47 m = 1.3;
48 c = 4;
49 x = 1:30;
50 ca = ao(plist('dtype', 'cdata', 'yvals', c.*ones(size(x))));
51 xa = ao(plist('dtype', 'cdata', 'yvals', x));
52
53 n = 1*randn(size(x));
54 sn = 1;
55 na = sn.*ao(plist('dtype', 'cdata', 'yvals', n));
56 ya = m.*xa + ca + na;
57
58 % Fit with straightLineFit
59 vn = var(na.y);
60 W = 1./vn*ones(size(x));
61 fy = straightLineFit(ya, plist('Weights', W));
62 sigma = find(fy.procinfo, 'Sigma');
63 params = find(sigma.procinfo, 'P')
64 Pstd = 100*find(sigma.procinfo, 'STD')./params
65 Pcov = find(sigma.procinfo, 'cov')
66
67 % theoretcical errors
68 y = ya.y.';
69 delta = (sum(W).*sum(W.*(x.^2)) - sum(W.*x).^2 );
70 tc = (sum(W.*(x.^2)) .* sum(W.*y) - sum(W.*x).*sum(W.*x.*y))./delta
71
72 find(sigma.procinfo, 'STD')
73 em = sqrt(sum(W)./delta)
74 ec = sqrt(sum(W.*(x.^2)) ./ delta)
75
76 params
77 Pstd = 100.*find(sigma.procinfo, 'STD')./params
78
79
80 %% curvefit needs an xy data
81 xy = ya.convert('to xydata');
82 pl = plist('Function', 'P(1).*Xdata + P(2)', 'P0', [1 0], 'Yerr', vn);
83 b = curvefit(xy, pl);
84
85 %%
86
87 M = [];
88 C = [];
89 dM = [];
90 dC = [];
91
92 noise = logspace(-2,1,100);
93
94 m = 1.3;
95 c = 4;
96 x = 1:30;
97 ca = ao(plist('dtype', 'cdata', 'yvals', c.*ones(size(x))));
98 xa = ao(plist('dtype', 'cdata', 'yvals', x));
99
100 for k=1:length(noise)
101
102 k
103
104 n = noise(k)*randn(size(x));
105 na = ao(plist('dtype', 'cdata', 'yvals', n));
106 ya = m.*xa + ca + na;
107
108 % Fit with straightLineFit
109
110 fy = straightLineFit(ya);
111 sigma = find(fy.procinfo, 'Sigma');
112 params = find(sigma.procinfo, 'P');
113 std = find(sigma.procinfo, 'STD');
114 cov = find(sigma.procinfo, 'cov');
115
116 M = [M params(1)];
117 C = [C params(2)];
118 dM = [dM std(1)];
119 dC = [dC std(2)];
120
121
122 end
123
124 iplot(fy,ya, plist('YerrL', {sigma,[]}, 'YerrU', {sigma,[]}))
125
126 %% Plot
127
128 figure
129 subplot(2,1,1)
130 shadedplot(noise, M-dM, M+dM, [0.6 0.6 0.6], [1 1 1]);
131 set(gca, 'XScale', 'log')
132 hold on
133 semilogx(noise, M);
134 xscale('log')
135 xlabel('Noise level (sigma)');
136 ylabel('gradient');
137
138 subplot(2,1,2)
139 shadedplot(noise, C-dC, C+dC, [0.6 0.6 0.6], [1 1 1]);
140 set(gca, 'XScale', 'log')
141 hold on
142 semilogx(noise, C);
143 xscale('log')
144 xlabel('Noise level (sigma)');
145 ylabel('Intercept');
146
147
148 % END