Conint

Conint Commit Details

Date:2012-05-15 22:13:12 (6 years 7 months ago)
Author:Michael Baudin
Commit:1
Message:Initial import.
Changes:
A/macros/cleanmacros.sce
A/help/en_US/conint_bernoullip.xml
A/etc/conint.quit
A/macros/names
A/help/en_US/conint_normmu.xml
A/tests/unit_tests/bernoullip.tst
A/tests/unit_tests/normmu.tst
A/macros/lib
A/help/en_US/conint_normvar.xml
A/help/en_US/conint_bernoullipnumber.xml
A/help/builder_help.sce
A/macros/conint_normvar.sci
A/help/en_US/conint_normmunumber.xml
A/tests/unit_tests/normvar.tst
A/license.txt
A/tests/unit_tests/bernoullipnumber.tst
A/macros/buildmacros.sce
A/tests/unit_tests/normmunumber.tst
A/tests
A/help/en_US/update_help.sce
A/help/en_US/build_help.sce
A/cleaner.sce
A/macros
A/macros/conint_bernoullip.sci
A/macros/conint_normmu.sci
A/help
A/builder.sce
A/macros/conint_bernoullipnumber.sci
A/etc/conint.start
A/help/en_US
A/macros/conint_normmunumber.sci
A/tests/unit_tests
A/changelog.txt
A/readme.txt
A/etc

File differences

cleaner.sce
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
// This file is released under the 3-clause BSD license. See COPYING-BSD.
// Generated by builder.sce: Please, do not edit this file
try
getversion("scilab");
catch
error("Scilab 5.0 or more is required.");
end
function perform_clean()
root_tlbx = get_absolute_file_path('cleaner.sce');
if isfile(root_tlbx + '/macros/cleanmacros.sce') then
exec(root_tlbx+'/macros/cleanmacros.sce');
end
if isfile(root_tlbx + '/src/cleaner_src.sce') then
exec(root_tlbx+'/src/cleaner_src.sce');
end
if isfile(root_tlbx + "/sci_gateway/cleaner_gateway.sce") then
exec(root_tlbx + "/sci_gateway/cleaner_gateway.sce");
mdelete(root_tlbx + "/sci_gateway/cleaner_gateway.sce");
end
if isfile(root_tlbx + "/help/cleaner_help.sce") then
exec(root_tlbx + "/help/cleaner_help.sce");
end
if isfile(root_tlbx + "/loader.sce") then
mdelete(root_tlbx + "/loader.sce");
end
endfunction
perform_clean();
clear perform_clean;
tests/unit_tests/bernoullip.tst
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
// Copyright (C) 2011 - Michael Baudin
//
// This file must be used under the terms of the GNU LGPL license.
// From Sheldon Ross, Example 7.5a, p.262
// We experimented with 100 transistors :
// 80 transistors works.
// Compute a 95% confidence interval for p,
// the probability that a transistor works.
pe = 80./100.;
n = 100;
[low, up] = conint_bernoullip ( pe , n , 1.-0.95 );
assert_checkalmostequal(low,0.7216014,1.e-5);
assert_checkalmostequal(up,0.8783986,1.e-5);
// Then (0.7216,0.8774) is a 95% C.I. for p.
// From Gilbert Saporta
// Section 13.5.4 Intervalle de confiance
// pour une proportion p
n = 400;
pe = 0.36;
[low, up] = conint_bernoullip ( pe , n , 1.-0.95 );
assert_checkalmostequal(low,0.3129609,1.e-5);
assert_checkalmostequal(up,0.4070391,1.e-5);
// Then (0.31,0.41) is a 95% C.I. for p.
// Check default value of level
n = 400;
pe = 0.36;
[low, up] = conint_bernoullip ( pe , n );
assert_checkalmostequal(low,0.3129609,1.e-5);
assert_checkalmostequal(up,0.4070391,1.e-5);
// Try various methods
// From Statistics: problems and solutions
// Edward Eryl Bassett,J. M. Bremner,B. J. T. Morgan
// "4C Binomial and Poisson distributions"
// "4C3 Women investors in building societies"
// p.148
pe = 22./80.;
n = 80;
level = 1.-0.90;
twosided = %t;
[low, up] = conint_bernoullip ( pe , n , level , twosided );
// Then (0.193,0.357) is a 90% C.I. for p.
assert_checkalmostequal(low,0.1928859,1.e-5);
assert_checkalmostequal(up,0.3571141,1.e-5);
// The quadratic or Wilson (1927) interval :
method = 2;
[low, up] = conint_bernoullip ( pe , n , level , twosided , method );
assert_checkalmostequal(low,0.2012659,1.e-5);
assert_checkalmostequal(up,0.3634549,1.e-5);
// Then (0.201,0.363) is an approximate 90% C.I. for p.
// Invert the Binomial: the Clopper-Pearson interval
method = 3;
[low, up] = conint_bernoullip ( pe , n , level , twosided , method );
assert_checkalmostequal(low,0.1884600,1.e-5);
assert_checkalmostequal(up,0.3522612,1.e-5);
// Then (0.188,0.352) is an approximate 90% C.I. for p.
tests/unit_tests/normmu.tst
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
// Copyright (C) 2011 - Michael Baudin
//
// This file must be used under the terms of the GNU LGPL license.
// Sheldon Ross, Example 7.3a, page 241
x = [5, 8.5, 12, 15, 7, 9, 7.5, 6.5, 10.5];
n = size(x,"*");
me = mean(x);
v = 4.;
[ low, up ] = conint_normmu ( n , me, v);
assert_checkalmostequal(me,9.,1.e-5);
assert_checkalmostequal(low,7.6933573,1.e-5);
assert_checkalmostequal(up,10.306643,1.e-5);
// Then P(7.69 <= mu <= 10.31) = 0.95
//
// Get one-sided intervals, Example 7.3b, p. 242
twosided = %f;
[ low, up ] = conint_normmu ( n , me, v, 1-0.95 , twosided );
assert_checkalmostequal(me,9.,1.e-5);
assert_checkalmostequal(low,7.9034309,1.e-5);
assert_checkalmostequal(up,10.096569,1.e-5);
// Then : P(7.903 <= mu) = 0.95
// P(mu <= 10.097) = 0.95
// Example 7.3c, p. 244
// Two-sided 99% interval, expected = (7.28, 10.72)
[ low, up ] = conint_normmu ( n , me, v, 1.-0.99 );
assert_checkalmostequal(me,9.,1.e-5);
assert_checkalmostequal(low,7.2827805,1.e-5);
assert_checkalmostequal(up,10.71722,1.e-5);
// One-sided 99% intervals are (7.447,inf), (−inf, 10.553)
[ low, up ] = conint_normmu ( n , me, v, 1.-0.99 , %f );
assert_checkalmostequal(me,9.,1.e-5);
assert_checkalmostequal(low,7.4491014,1.e-5);
assert_checkalmostequal(up,10.550899,1.e-5);
// EXAMPLE 7.3e, p. 247
// The variance is unknown: use Student's T random variable
v = variance(x);
[ low, up ] = conint_normmu ( n , me, v, 1.-0.95 , [] , %f );
assert_checkalmostequal(me,9.,1.e-5);
assert_checkalmostequal(low,6.630806,1.e-5);
assert_checkalmostequal(up,11.369194,1.e-5);
// Then P(6.63<= mu <= 11.37) = 0.95
// EXAMPLE 7.3f, p.249
x = [54, 63, 58, 72, 49, 92, 70, 73, ..
69, 104, 48, 66, 80, 64, 77];
n = size(x,"*");
me = mean(x);
v = variance(x);
// Two-sided :
[ low, up ] = conint_normmu ( n , me, v, 1.-0.95 , [] , %f );
assert_checkalmostequal(me,69.266667,1.e-5);
assert_checkalmostequal(low,60.866937,1.e-5);
assert_checkalmostequal(up,77.666397,1.e-5);
// Then P(60.865<= mu <=77.6683) = 0.95
// One-sided :
[ low, up ] = conint_normmu ( n , me, v, 1.-0.95 , %f , %f );
assert_checkalmostequal(me,69.266667,1.e-5);
assert_checkalmostequal(low,62.368764,1.e-5);
assert_checkalmostequal(up,76.164569,1.e-5);
// Then : P (mu <= 76.16) = 0.95
// Then : P (62.368 <= mu) = 0.95
// EXAMPLE 7.3g, Monte-Carlo simulation
grand("setsd",0);
U = grand(100,1,"def");
x = sqrt(1-U.^2);
n = size(x,"*");
me = mean(x);
v = variance(x);
[ low, up ] = conint_normmu ( n , me, v, 1.-0.95 , [] , %f );
// Exact integral is :
assert_checkalmostequal(me,%pi/4,1.e-1);
assert_checktrue(me>low);
assert_checktrue(me<up);
// Check default value of level
x = [5, 8.5, 12, 15, 7, 9, 7.5, 6.5, 10.5];
n = size(x,"*");
me = mean(x);
v = variance(x);
[ low, up ] = conint_normmu ( n , me, v , [] , [] , %f );
assert_checkalmostequal(me,9.,1.e-5);
assert_checkalmostequal(low,6.630806,1.e-5);
assert_checkalmostequal(up,11.369194,1.e-5);
tests/unit_tests/bernoullipnumber.tst
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
// Copyright (C) 2011 - Michael Baudin
//
// This file must be used under the terms of the GNU LGPL license.
// From Gilbert Saporta,
// Section 13.5.4 Intervalle de confiance
// pour une proportion p
// p.313
len = 0.05*2;
level = 1.-0.90;
n = conint_bernoullipnumber ( len , level );
assert_checkequal(n,271);
// Reproduce the table p 313.
mprintf("level = %9s %9s %9s %9s", ..
" ", "1.-0.90", "1.-0.95", "1.-0.98")
for len = 2*[0.01 0.02 0.05]
mprintf("len = %9s, n = ", string(len))
for level = 1.-[0.90 0.95 0.98]
n = conint_bernoullipnumber ( len , level );
mprintf("%9s ", string(n));
end
mprintf("\n");
end
// From Sheldon Ross,
// Example 7.5c
len = 0.05;
// Compute n for 95% C.I.
n = conint_bernoullipnumber ( len );
assert_checkequal(n,1537);
// Compute n for 99% C.I.
level = 1.-0.99;
n = conint_bernoullipnumber ( len , level );
assert_checkequal(n,2654);
// Compute n for 99% C.I. and probability estimate
pe = 26. / 30.;
n = conint_bernoullipnumber ( len , level , pe );
assert_checkequal(n,1227);
tests/unit_tests/normvar.tst
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
// Copyright (C) 2011 - Michael Baudin
//
// This file must be used under the terms of the GNU LGPL license.
// Sheldon Ross, Example 7.3h, page 252
// Compute a 90% confidence interval for the
// variance of a normal variable.
x = [0.123 0.124 0.126 0.120 0.130 ...
0.133 0.125 0.128 0.124 0.126];
n = size(x,"*");
v = variance(x);
[ low, up ] = conint_normvar ( n , v , 1.-0.90 );
assert_checkalmostequal(v,1.36555556D-05,1.e-5);
assert_checkalmostequal(low,7.26403231D-06,1.e-5);
assert_checkalmostequal(up,3.69611516D-05,1.e-5);
// Then P(7.26403231D-06 <= var(X) <= 36.9611516D-06) = 0.90
// Example from Gilbert Saporta,
// Section 13.5.3, "Variance d'une loi normale"
// Section 13.5.3.2, "m est inconnu"
x = [6;2;2;8;12;6;6;13;..
8;2;10;0;7;5;4;4;3;..
7;4;9;6;2;-3;5;7;2;4;5;2;2];
n = size(x,"*");
v = variance(x);
[ low, up ] = conint_normvar ( n , v , 1.-0.90 )
assert_checkalmostequal(v,12.,1.e-5);
assert_checkalmostequal(low,8.1772743,1.e-5);
assert_checkalmostequal(up,19.651728,1.e-5);
// Check default value of level
x = [0.123 0.124 0.126 0.120 0.130 ...
0.133 0.125 0.128 0.124 0.126];
n = size(x,"*");
v = variance(x);
[ low, up ] = conint_normvar ( n , v );
assert_checkalmostequal(v,1.36555556D-05,1.e-5);
assert_checkalmostequal(low,6.46067919D-06,1.e-5);
assert_checkalmostequal(up,4.55119530D-05,1.e-5);
tests/unit_tests/normmunumber.tst
1
2
3
4
5
6
7
8
9
10
11
// Copyright (C) 2011 - Michael Baudin
//
// This file must be used under the terms of the GNU LGPL license.
// Sheldon Ross, Example 7.3a, page 241
n = conint_normmunumber ( 0.2 , 0.3^2 , 1.-0.95 );
assert_checkalmostequal(n,35);
// Check default value of level
n = conint_normmunumber ( 0.2 , 0.3^2 );
assert_checkalmostequal(n,35);
macros/buildmacros.sce
1
2
3
4
5
6
7
8
9
10
11
// Copyright (C) 2010 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
tbx_build_macros("accsum", get_absolute_file_path('buildmacros.sce'));
clear tbx_build_macros;
macros/conint_normmunumber.sci
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
// Copyright (C) 2011 - Michael Baudin
//
// This file must be used under the terms of the GNU LGPL license.
function n = conint_normmunumber ( varargin )
// Number of experiments for a Normal variable.
//
// Calling Sequence
// n = conint_normmunumber ( len , v )
// n = conint_normmunumber ( len , v , level )
//
// Parameters
// len : a 1-by-1 matrix of doubles, positive, the length of the confidence interval.
// v : a 1-by-1 matrix of doubles, the variance
// level : a 1-by-1 matrix of doubles, the confidence level (default level = 1.-0.95=0.05). level is expected to be in the range [0.,0.5]
// n : a 1-by-1 matrix of doubles, integer value, the number of experiments
//
// Description
// Returns the number of experiments which guarantees that the
// length of the confidence interval is less than <literal>len</literal>.
// This two-sided confidence interval is for the expectation of the
// random variable and we make the assumption that the random
// variable has a Normal distribution.
//
// In other words, if <literal>x</literal> is a vector
// with <literal>n</literal> entries, then
//
// <screen>
// [ m, low, up ] = conint_normmu ( x , level , v )
// </screen>
//
// is so that
//
// <literal>
// abs(up-low) <= len
// </literal>
//
// To compute <literal>n</literal>, we assume that the random variable
// has Normal distribution with known variance <literal>v</literal>.
//
// Examples
// // Sheldon Ross, EXAMPLE 7.3d, p245
// // Compute the number of experiments for a random variable
// // with standard deviation 0.3 so that
// // the expectation has a two-sided 95% confidence interval
// // with length no greater than 0.2.
// n = conint_normmunumber ( 0.2 , 0.3^2 , 1.-0.95 )
// // Expected: n=35
//
// Authors
// Copyright (C) 2011 - Michael Baudin
//
// Bibliography
// http://en.wikipedia.org/wiki/Confidence_interval
// "Introduction to probability and statistics for engineers and scientists", Sheldon Ross, Third Edition, 2004
[lhs, rhs] = argn();
apifun_checkrhs ( "conint_normmunumber" , rhs , 2 : 3 )
apifun_checklhs ( "conint_normmunumber" , lhs , 0 : 1 )
//
// Get arguments
len = varargin(1)
v = varargin(2)
level = apifun_argindefault ( varargin , 3 , 1.-0.95 )
//
// Check type
apifun_checktype ( "conint_normmunumber" , len , "len" , 1 , "constant" )
apifun_checktype ( "conint_normmunumber" , v, "v" , 2 , "constant" )
apifun_checktype ( "conint_normmunumber" , level, "level" , 3 , "constant" )
//
// Check size
apifun_checkscalar ( "conint_normmunumber" , len , "len" , 1 )
apifun_checkscalar ( "conint_normmunumber" , v , "v" , 2 )
apifun_checkscalar ( "conint_normmunumber" , level , "level" , 3 )
//
// Check content
apifun_checkgreq ( "conint_normmunumber" , len , "len" , 1 , 0. )
apifun_checkgreq ( "conint_normmunumber" , v , "v" , 2 , 0. )
apifun_checkrange ( "conint_normmunumber" , level , "level" , 3 , 0. , 0.5 )
//
p = (2-level)/2
q = level/2
f = cdfnor("X",0.,1.,p,q)
n = v * (2.*f/len)^2
n = ceil(n)
endfunction
macros/cleanmacros.sce
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
// ====================================================================
// Allan CORNET
// DIGITEO 2009
// This file is released into the public domain
// ====================================================================
libpath = get_absolute_file_path('cleanmacros.sce');
binfiles = ls(libpath+'/*.bin');
for i = 1:size(binfiles,'*')
mdelete(binfiles(i));
end
mdelete(libpath+'/names');
mdelete(libpath+'/lib');
// ====================================================================
macros/conint_bernoullip.sci
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
// Copyright (C) 2011 - Michael Baudin
//
// This file must be used under the terms of the GNU LGPL license.
function [ low, up ] = conint_bernoullip ( varargin )
// C.I. of the probability of a Bernoulli variable.
//
// Calling Sequence
// [low, up] = conint_bernoullip ( pe , n )
// [low, up] = conint_bernoullip ( pe , n , level )
// [low, up] = conint_bernoullip ( pe , n , level , twosided )
// [low, up] = conint_bernoullip ( pe , n , level , twosided , method )
//
// Parameters
// pe : a 1-by-1 matrix of doubles, an estimate of the probability
// n : a 1-by-1 matrix of doubles, integer value, positive, the number of trials
// level : a 1-by-1 matrix of doubles, the confidence level (default level = 1.-0.95=0.05). level is expected to be in the range [0.,0.5]
// twosided : a 1-by-1 matrix of booleans, the side of the interval (default twosided = %t). If twosided is true, then low and up is a two-sided interval [low,up]. If twosided is false, then the intervals are [-inf,up] and [low,inf].
// method : a 1-by-1 matrix of doubles, integer value, the method to estimate the C.I. Available values are method=1,2,3 (default method = 1). If method=1, then approximates the Binomial distribution with the Normal distribution. If method=2, then
// low : a 1-by-1 matrix of doubles, the estimated lower bound
// up : a 1-by-1 matrix of doubles, the estimated upper bound
//
// Description
// Computes a confidence interval of the probability
// <literal>p</literal> of a Bernoulli variable.
//
// This function makes the assumption that the data
// in <literal>x</literal> has a binomial distribution with
// parameters <literal>p</literal> and
// <literal>n</literal>.
// In other words, we consider that the random variable is a "success" with probability
// <literal>p</literal> and is a failure with probability
// <literal>1-p</literal>.
// We assume that we have performed <literal>n</literal> trials
// of this experiment.
// We assume that <literal>pe</literal>
// is the fraction of the <literal>n</literal> trials which are successes.
// This function computes confidence intervals for <literal>p</literal>.
//
// If <literal>twosided</literal> is true, then
// <literal>low</literal> and <literal>up</literal> are such that
// the probability <literal>p</literal> satisfies the equality
//
// <literal>
// P(low < p < up) = 1-level
// </literal>
//
// If <literal>twosided</literal> is false, then
// <literal>low</literal> and <literal>up</literal> are such that
// the probability <literal>p</literal> satisfies the equality
//
// <literal>
// P(p < up) = 1-level
// </literal>
//
// and
//
// <literal>
// P(low < p) = 1-level
// </literal>
//
// The confidence interval that we compute is approximate.
// It is based on the Central Limit Theorem.
// We use a Normal distribution to compute the confidence interval.
//
// If method=1, then approximates the Binomial distribution
// with the Normal distribution.
// Then computes approximate roots of the quadratic equation.
// This is the "Textbook" formula.
//
// If method=2, then approximates the Binomial distribution with the Normal distribution.
// Then computes exact roots of the quadratic equation.
// This is the Wilson (1927).
//
// If method=3, then invert the quantiles of the
// Binomial distribution with the Normal distribution.
// This is the Clopper-Pearson interval.
//
// Examples
// // From Sheldon Ross, Example 7.5a, p.262
// // We experimented with 100 transistors :
// // 80 transistors works.
// // Compute a 95% confidence interval for p,
// // the probability that a transistor works.
// pe = 80./100.;
// n = 100;
// [low, up] = conint_bernoullip ( pe , n , 1.-0.95 )
// // Then (0.7216,0.8774) is an approximate 95% C.I. for p.
//
// // From Gilbert Saporta
// // Section 13.5.4 Intervalle de confiance
// // pour une proportion p
// n = 400;
// pe = 0.36;
// [low, up] = conint_bernoullip ( pe , n , 1.-0.95 )
// // Then (0.31,0.41) is an approximate 95% C.I. for p.
//
// // Try various methods
// // From Statistics: problems and solutions
// // Edward Eryl Bassett,J. M. Bremner,B. J. T. Morgan
// // "4C Binomial and Poisson distributions"
// // "4C3 Women investors in building societies"
// // p.148
// pe = 22./80.;
// n = 80;
// level = 1.-0.90
// twosided = %t
// [low, up] = conint_bernoullip ( pe , n , level , twosided )
// // Then (0.193,0.357) is a 90% C.I. for p.
// // The quadratic or Wilson (1927) interval :
// method = 2;
// [low, up] = conint_bernoullip ( pe , n , level , twosided , method )
// // Then (0.201,0.363) is an approximate 90% C.I. for p.
// // Invert the Binomial: the Clopper-Pearson interval
// method = 3
// [low, up] = conint_bernoullip ( pe , n , level , twosided , method )
// // Then (0.188,0.352) is an approximate 90% C.I. for p.
//
// Authors
// Copyright (C) 2011 - Michael Baudin
//
// Bibliography
// http://en.wikipedia.org/wiki/Confidence_interval
// "Introduction to probability and statistics for engineers and scientists", Sheldon Ross, Third Edition, 2004
// "Probabilités, Analyse de Données et Statistique", Gilbert Saporta, 2nd Ed., 2006
// http://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval
[lhs, rhs] = argn();
apifun_checkrhs ( "conint_bernoullip" , rhs , 2 : 5 )
apifun_checklhs ( "conint_bernoullip" , lhs , 2 : 2 )
//
// Get arguments
pe = varargin(1)
n = varargin(2)
level = apifun_argindefault ( varargin , 3 , 1.-0.95 )
twosided = apifun_argindefault ( varargin , 4 , %t )
method = apifun_argindefault ( varargin , 5 , 1 )
//
// Check type
apifun_checktype ( "conint_bernoullip" , pe , "pe" , 1 , "constant" )
apifun_checktype ( "conint_bernoullip" , n , "n" , 2 , "constant" )
apifun_checktype ( "conint_bernoullip" , level, "level" , 3 , "constant" )
apifun_checktype ( "conint_bernoullip" , twosided, "twosided" , 4 , "boolean" )
apifun_checktype ( "conint_bernoullip" , method, "method" , 5 , "constant" )
//
// Check size
apifun_checkscalar ( "conint_bernoullip" , pe , "pe" , 1 )
apifun_checkscalar ( "conint_bernoullip" , n , "n" , 2 )
apifun_checkscalar ( "conint_bernoullip" , level , "level" , 3 )
apifun_checkscalar ( "conint_bernoullip" , twosided , "twosided" , 4 )
apifun_checkscalar ( "conint_bernoullip" , method , "method" , 5 )
//
// Check content
apifun_checkrange ( "conint_bernoullip" , pe , "pe" , 1 , 0. , 1. )
apifun_checkgreq ( "conint_bernoullip" , n , "n" , 2 , 1 )
apifun_checkflint ( "conint_bernoullip" , n , "n" , 2 )
apifun_checkrange ( "conint_bernoullip" , level , "level" , 3 , 0. , 0.5 )
apifun_checkoption ( "conint_bernoullip" , method , "method" , 5 , [1,2,3] )
//
select method
case 1 then
// The Normal approximation to the Binomial.
// Moreover, computes approximate roots of the quadratic
// polynomial.
if ( twosided ) then
q = level/2.
p = 1.-q
else
q = level
p = 1.-q
end
f = cdfnor("X",0.,1.,p,q)
low = pe - f * sqrt(pe*(1.-pe)/n)
up = pe + f * sqrt(pe*(1.-pe)/n)
case 2 then
// The quadratic or Wilson (1927) interval :
// computes the exact roots of the quadratic polynomial
// which appears when approximating the Binomial distribution
// with a Normal distribution.
if ( twosided ) then
q = level/2.
p = 1.-q
else
q = level
p = 1.-q
end
f = cdfnor("X",0.,1.,p,q)
c = []
c(3) = pe^2
c(2) = -(f^2/n + 2*pe)
c(1) = 1+f^2/n
proot = roots(c)
if (or(imag(proot)<>0.)) then
error(msprintf("%s: Complex roots (%s,%s).",..
"conint_bernoullip",string(proot(1)),,string(proot(2))))
end
proot = real(proot)
low = min(proot)
up = max(proot)
case 3 then
// Invert the Binomial: the Clopper-Pearson interval
if ( twosided ) then
q = level/2.
p = 1.-q
else
q = level
p = 1.-q
end
S=cdfbin("S",n,pe,1.-pe,p,q)
up = S/n
if ( twosided ) then
p = level/2.
q = 1.-p
else
p = level
q = 1.-p
end
S=cdfbin("S",n,pe,1.-pe,p,q)
low = S/n
else
error(msprintf("%s: Unknown method %s.","conint_bernoullip",string(method)))
end
endfunction
macros/names
1
2
3
4
5
conint_normvar
conint_normmunumber
conint_normmu
conint_bernoullipnumber
conint_bernoullip
macros/conint_normmu.sci
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
// Copyright (C) 2011 - Michael Baudin
//
// This file must be used under the terms of the GNU LGPL license.
function [ low, up ] = conint_normmu ( varargin )
// C.I. of the mean of a normal variable.
//
// Calling Sequence
// [ low, up ] = conint_normmu ( n , me , v )
// [ low, up ] = conint_normmu ( n , me , v , level )
// [ low, up ] = conint_normmu ( n , me , v , level , twosided )
// [ low, up ] = conint_normmu ( n , me , v , level , twosided , known )
//
// Parameters
// n : a n-by-1 or 1-by-n matrix of doubles, integer value, positive, the number of samples
// me : a 1-by-1 matrix of doubles, the mean estimate
// level : a 1-by-1 matrix of doubles, the confidence level (default level = 1.-0.95=0.05). level is expected to be in the range <literal>[0.,0.5]</literal>
// v : a 1-by-1 matrix of doubles, the variance. By default, the variance is assumed to be exact (i.e. by default known=%t). If the variance <literal>v</literal> is estimated, then set the <literal>known</literal> option to %f.
// twosided : a 1-by-1 matrix of booleans, the side of the interval (default twosided = %t). If <literal>twosided</literal> is %t, then <literal>low</literal> and <literal>up</literal> make a two-sided interval <literal>[low,up]</literal>. If <literal>twosided</literal> is %f, then the intervals are <literal>[-inf,up]</literal> and <literal>[low,inf]</literal>.
// known : a 1-by-1 matrix of booleans, %t if the variance <literal>v</literal> is exact (default known = %t). If <literal>known</literal> is %t, then <literal>v</literal> is the exact variance. If known is %f, then <literal>v</literal> is the estimate of the variance.
// low : a 1-by-1 matrix of doubles, the estimated lower bound
// up : a 1-by-1 matrix of doubles, the estimated upper bound
//
// Description
// Computes a confidence interval of the mean of a normal variable, that is,
// computes <literal>m</literal> as the mean of <literal>x</literal> and
// computes confidence intervals for <literal>mu(X)</literal>.
//
// If <literal>twosided</literal> is true, then
// <literal>low</literal> and <literal>up</literal> are such that
// the expectation <literal>mu(X)</literal> is such that
//
// <literal>
// P(low < mu(X) < up) = 1-level
// </literal>
//
// If <literal>twosided</literal> is false, then
// <literal>low</literal> and <literal>up</literal> are such that
// the expectation <literal>mu</literal> is such that
//
// <literal>
// P(mu(X) < up) = 1-level
// </literal>
//
// and
//
// <literal>
// P(low < mu(X)) = 1-level
// </literal>
//
// This function makes the assumption that the
// data in x has a normal distribution with expectation <literal>mu(X)</literal> and
// variance <literal>var(X)</literal>.
//
// If the variance is exact (i.e. <literal>known=%t</literal>)
// then we use a Normal distribution to compute the confidence interval.
//
// If the variance is estimated (i.e. if <literal>known=%f</literal>),
// then we use Student's T distribution to compute the confidence interval.
//
// Examples
// // Sheldon Ross, Example 7.3a, page 241
// x = [5, 8.5, 12, 15, 7, 9, 7.5, 6.5, 10.5]
// n = size(x,"*")
// me = mean(x)
// // By default, we compute a two-sided 95% C.I.
// // with exact variance:
// v = 4.
// [ low, up ] = conint_normmu ( n , me , v )
// low_expected = 7.69
// up_expected = 10.31
// // Then P(7.69 <= mu <= 10.31) = 0.95
//
// // Get one-sided 95% C.I.
// // Example 7.3b, p. 242
// low_expected = 7.903
// up_expected = 10.097
// twosided = %f
// [ low, up ] = conint_normmu ( n , me , v , [] , twosided )
// // Then : P(7.903 <= mu) = 0.95
// // P(mu <= 10.097) = 0.95
//
// // Example 7.3c, p. 244
// // Two-sided 99% interval : (7.28, 10.72)
// [ low, up ] = conint_normmu ( n , me , v , 1.-0.99 )
// // One-sided 99% intervals : (7.447,inf), (−inf, 10.553)
// [ low, up ] = conint_normmu ( n , me , v , 1.-0.99 , %f )
//
// // EXAMPLE 7.3e, p. 247
// // The variance is unknown: use Student's T random variable
// v = variance(x)
// [ low, up ] = conint_normmu ( n , me , v , [] , [] , %f )
// // Then P(6.63<= mu <= 11.37) = 0.95
//
// // EXAMPLE 7.3f, p.249
// x = [54, 63, 58, 72, 49, 92, 70, 73, ..
// 69, 104, 48, 66, 80, 64, 77];
// n = size(x,"*")
// me = mean(x)
// v = variance(x)
// // Two-sided :
// [ low, up ] = conint_normmu ( n , me , v , 1.-0.95 , [] , %f )
// // Then P(60.865<= mu <=77.6683) = 0.95
// // One-sided :
// [ low, up ] = conint_normmu ( n , me , v , 1.-0.95 , %f , %f )
// // Then : P (mu <= 76.16) = 0.95
// // Then : P (62.368 <= mu) = 0.95
//
// // EXAMPLE 7.3g, Monte-Carlo simulation
// U = grand(100,1,"def");
// x = sqrt(1-U.^2);
// n = size(x,"*")
// me = mean(x)
// v = variance(x)
// [ low, up ] = conint_normmu ( n , me , v , [] , [] , %f )
// // Exact integral is :
// %pi/4
//
// Authors
// Copyright (C) 2011 - Michael Baudin
//
// Bibliography
// http://en.wikipedia.org/wiki/Confidence_interval
// "Introduction to probability and statistics for engineers and scientists", Sheldon Ross, Third Edition, 2004
[lhs, rhs] = argn();
apifun_checkrhs ( "conint_normmu" , rhs , 3 : 6 )
apifun_checklhs ( "conint_normmu" , lhs , 2 : 2 )
//
// Get arguments
n = varargin(1)
me = varargin(2)
v = varargin(3)
level = apifun_argindefault ( varargin , 4 , 1.-0.95 )
twosided = apifun_argindefault ( varargin , 5 , %t )
known = apifun_argindefault ( varargin , 6 , %t )
//
// Check type
apifun_checktype ( "conint_normmu" , n , "n" , 1 , "constant" )
apifun_checktype ( "conint_normmu" , me , "me" , 2 , "constant" )
apifun_checktype ( "conint_normmu" , v, "v" , 3 , "constant" )
apifun_checktype ( "conint_normmu" , level, "level" , 4 , "constant" )
apifun_checktype ( "conint_normmu" , twosided, "twosided" , 5 , "boolean" )
apifun_checktype ( "conint_normmu" , known, "known" , 6 , "boolean" )
//
// Check size
apifun_checkscalar ( "conint_normmu" , n , "n" , 1 )
apifun_checkscalar ( "conint_normmu" , me , "me" , 2 )
apifun_checkscalar ( "conint_normmu" , v , "v" , 3 )
apifun_checkscalar ( "conint_normmu" , level , "level" , 4 )
apifun_checkscalar ( "conint_normmu" , twosided , "twosided" , 5 )
apifun_checkscalar ( "conint_normmu" , known , "known" , 6 )
//
// Check content
apifun_checkflint ( "conint_normmu" , n , "n" , 1 )
apifun_checkgreq ( "conint_normmu" , n , "n" , 1 , 1. )
tiniest = number_properties("tiniest")
apifun_checkgreq ( "conint_normmu" , v , "v" , 3 , tiniest )
apifun_checkrange ( "conint_normmu" , level , "level" , 4 , 0. , 0.5 )
//
if ( twosided ) then
q = level/2.
p = 1.-q
else
q = level
p = 1.-q
end
s = sqrt(v)
if ( known ) then
// The variance is known:
// use a Standard Normal variable
f = cdfnor("X",0.,1.,p,q)
else
// The variance is unknown:
// use Student's T random variable.
f = cdft("T",n-1,p,q)
end
low = me - f * s/sqrt(n)
up = me + f * s/sqrt(n)
endfunction
macros/lib
macros/conint_bernoullipnumber.sci
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
// Copyright (C) 2011 - Michael Baudin
//
// This file must be used under the terms of the GNU LGPL license.
function n = conint_bernoullipnumber ( varargin )
// Number of experiments for a Bernoulli variable.
//
// Calling Sequence
// n = conint_bernoullipnumber ( len )
// n = conint_bernoullipnumber ( len , level )
// n = conint_bernoullipnumber ( len , level , pe )
//
// Parameters
// len : a 1-by-1 matrix of doubles, positive, the length of the required C.I.
// level : a 1-by-1 matrix of doubles, the confidence level (default level = 1.-0.95=0.05). level is expected to be in the range [0.,0.5]
// pe : a 1-by-1 matrix of doubles, an estimate of the probability of success (default pe = 0.5)
// n : a 1-by-1 matrix of doubles, the number of trials
//
// Description
// Returns the number of experiments which guarantees that the
// length of the confidence interval for the
// probability of success of a Bernoulli variable
// is less than <literal>len</literal>.
//
// In other words, if <literal>pe</literal> is the estimate of
// the probability of success
// with <literal>n</literal> trials, then
//
// <screen>
// [low, up] = conint_bernoullip ( pe , n )
// </screen>
//
// is so that
//
// <literal>
// abs(up-low) <= len
// </literal>
//
// To compute <literal>n</literal>, we assume that the random variable
// has Binomial distribution.
//
// This function makes the assumption that the data
// in <literal>x</literal> has a binomial distribution with
// parameters <literal>p</literal> and
// <literal>n</literal>.
//
// The default value <literal>pe=0.5</literal>
// guarantees that the number of experiments is large enough,
// no matter the actual value of the probability <literal>p</literal>.
// In other words, if <literal>pe</literal> is not provided
// be the user, the value of the output variable <literal>n</literal>
// is guaranteed, but might be much too large than required.
//
// Examples
// // From Gilbert Saporta,
// // Section 13.5.4 Intervalle de confiance
// // pour une proportion p
// // p.313
// len = 0.05*2
// level = 1.-0.90
// n = conint_bernoullipnumber ( len , level )
// // Therefore, it requires 271 experiments
// // to make sure that the 90% C.I. for p
// // has a length no greater than 0.1.
//
// // Reproduce the table p 313.
// mprintf("level = %9s %9s %9s %9s", ..
// " ", "1.-0.90", "1.-0.95", "1.-0.98")
// for len = 2*[0.01 0.02 0.05]
// mprintf("len = %9s, n = ", string(len))
// for level = 1.-[0.90 0.95 0.98]
// n = conint_bernoullipnumber ( len , level );
// mprintf("%9s ", string(n));
// end
// mprintf("\n");
// end
//
// // From Sheldon Ross,
// // Example 7.5c
// // Estimate the failure rate of computer chips.
// len = 0.05
// // Compute n for 95% C.I.
// n = conint_bernoullipnumber ( len )
// // Compute n for 99% C.I.
// level = 1.-0.99
// n = conint_bernoullipnumber ( len , level )
// // Compute n for 99% C.I. and probability estimate
// // From 30 sample computer chips, there are
// // 26. chips which works.
// pe = 26. / 30.
// n = conint_bernoullipnumber ( len , level , pe )
// // It may require 1227 experiments to get a
// // 99 C.I. with length less than 0.05.
//
// // How many experiments to estimate accurately a
// // small probability with 95% C.I. ?
// level = 1.-0.95;
// for pe = logspace(-1,-7,7)
// len = pe/10.;
// n = conint_bernoullipnumber ( len , level , pe );
// mprintf("pe = %-9s, len = %-9s, n = %-9s\n", ..
// string(pe), string(len), string(n))
// end
//
// Authors
// Copyright (C) 2011 - Michael Baudin
//
// Bibliography
// http://en.wikipedia.org/wiki/Confidence_interval
// "Introduction to probability and statistics for engineers and scientists", Sheldon Ross, Third Edition, 2004
// "Probabilités, Analyse de Données et Statistique", Gilbert Saporta, 2nd Ed., 2006
[lhs, rhs] = argn();
apifun_checkrhs ( "conint_bernoullipnumber" , rhs , 1 : 3 )
apifun_checklhs ( "conint_bernoullipnumber" , lhs , 0 : 3 )
//
// Get arguments
len = varargin(1)
level = apifun_argindefault ( varargin , 2 , 1.-0.95 )
pe = apifun_argindefault ( varargin , 3 , 0.5 )
//
// Check type
apifun_checktype ( "conint_bernoullipnumber" , len , "len" , 1 , "constant" )
apifun_checktype ( "conint_bernoullipnumber" , level, "level" , 2 , "constant" )
apifun_checktype ( "conint_bernoullipnumber" , pe, "pe" , 3 , "constant" )
//
// Check size
apifun_checkscalar ( "conint_bernoullipnumber" , len , "len" , 1 )
apifun_checkscalar ( "conint_bernoullipnumber" , level , "level" , 2 )
apifun_checkscalar ( "conint_bernoullipnumber" , pe , "pe" , 3 )
//
// Check content
apifun_checkgreq ( "conint_bernoullipnumber" , len , "len" , 1 , 0. )
apifun_checkrange ( "conint_bernoullipnumber" , level , "level" , 2 , 0. , 0.5 )
apifun_checkrange ( "conint_bernoullipnumber" , pe , "pe" , 1 , 0. , 1. )
//
q = level/2.
p = 1.-q
f = cdfnor("X",0.,1.,p,q)
n = (2*f)^2*pe*(1.-pe) / (len^2)
n = ceil(n)
endfunction
macros/conint_normvar.sci
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
// Copyright (C) 2011 - Michael Baudin
//
// This file must be used under the terms of the GNU LGPL license.
function [ low, up ] = conint_normvar ( varargin )
// C.I. of the variance of a normal variable.
//
// Calling Sequence
// [ low, up ] = conint_normvar ( n , v )
// [ low, up ] = conint_normvar ( n , v , level )
// [ low, up ] = conint_normvar ( n , v , level , twosided )
//
// Parameters
// n : a 1-by-1 matrix of doubles, positive, integer value, the number of samples
// v : a 1-by-1 matrix of doubles, positive, the variance estimate
// level : a 1-by-1 matrix of doubles, the confidence level (default level = 1.-0.95=0.05). level is expected to be in the range [0.,0.5]
// twosided : a 1-by-1 matrix of booleans, the side of the interval (default twosided = %t). If twosided is true, then low and up is a two-sided interval [low,up]. If twosided is false, then the intervals are [-inf,up] and [low,inf].
// low : a 1-by-1 matrix of doubles, the estimated lower bound
// up : a 1-by-1 matrix of doubles, the estimated upper bound
//
// Description
// Computes a confidence interval of the variance of a normal variable, that is,
// computes <literal>v</literal> as the variance of <literal>x</literal> and
// computes confidence intervals for <literal>v</literal>.
//
// If <literal>twosided</literal> is true, then
// <literal>low</literal> and <literal>up</literal> are such that
// the variance <literal>var(X)</literal> is such that
//
// <literal>
// P(low < var < up) = 1-level
// </literal>
//
// If <literal>twosided</literal> is false, then
// <literal>low</literal> and <literal>up</literal> are such that
// the variance <literal>var(X)</literal> is such that
//
// <literal>
// P(var(X) < up) = 1-level
// </literal>
//
// and
//
// <literal>
// P(low < var(X)) = 1-level
// </literal>
//
// This function makes the assumption that the
// data in <literal>x</literal> has a normal distribution with
// expectation <literal>mu(X)</literal> and
// variance <literal>var(X)</literal>.
//
// We use a Chi-square distribution to compute the confidence interval.
//
// Examples
// // Sheldon Ross, Example 7.3h, page 252
// // Compute a 90% confidence interval for the
// // variance of a normal variable.
// x = [0.123 0.124 0.126 0.120 0.130 ...
// 0.133 0.125 0.128 0.124 0.126]
// n = size(x,"*")
// v = variance(x)
// [ low, up ] = conint_normvar ( n , v , 1.-0.90 )
// // Then P(7.264D-06 <= var(X) <= 36.96D-06) = 0.90
// // Give the confidence interval for the standard
// // deviation
// [sqrt(low) sqrt(v) sqrt(up)]
//
// // Example from Gilbert Saporta,
// // Section 13.5.3, "Variance d'une loi normale"
// // Section 13.5.3.2, "m est inconnu"
// // Get the C.I. for the variance of a normal
// // variable with n=30 and var(X)=12.
// x = [6;2;2;8;12;6;6;13;..
// 8;2;10;0;7;5;4;4;3;..
// 7;4;9;6;2;-3;5;7;2;4;5;2;2];
// n = size(x,"*")
// v = variance(x)
// [ low, up ] = conint_normvar ( n , v , 1.-0.90 )
// // Then P(8.177 < var(X) < 19.65) = 0.90
// // The 90% C.I. is (8.177,19.65)
// // Notice that Saporta gets (8.46,20.33).
// // This is because Saporta uses n as the denominator
// // in the variance estimate, while we use (n-1).
//
// // Example from
// // "Probability and Statistics for Engineers and Scientists",
// // Walpole, Myers, Myers, Ye
// // Example 9.17:
// x = [46.4, 46.1, 45.8, 47.0, 46.1, ..
// 45.9, 45.8, 46.9, 45.2, 46.0];
// n = size(x,"*")
// v = variance(x)
// [ low, up ] = conint_normvar ( n , v , 1.-0.95 )
// // Then (0.135,0.953) is a 95% C.I. for the variance.
//
// Authors
// Copyright (C) 2011 - Michael Baudin
//
// Bibliography
// http://en.wikipedia.org/wiki/Confidence_interval
// "Introduction to probability and statistics for engineers and scientists", Sheldon Ross, Third Edition, 2004
// "Probabilités, Analyse de Données et Statistique", Gilbert Saporta, 2nd Ed., 2006
// // "Probability and Statistics for Engineers and Scientists", Ronald Walpole, Raymond Myers, Sharon Myers, Keying Ye, Eight Edition, 2006
[lhs, rhs] = argn();
apifun_checkrhs ( "conint_normvar" , rhs , 2 : 4 )
apifun_checklhs ( "conint_normvar" , lhs , 2 : 2 )
//
// Get arguments
n = varargin(1)
v = varargin(2)
level = apifun_argindefault ( varargin , 3 , 1.-0.95 )
twosided = apifun_argindefault ( varargin , 4 , %t )
//
// Check type
apifun_checktype ( "conint_normvar" , n , "n" , 1 , "constant" )
apifun_checktype ( "conint_normvar" , v , "v" , 2 , "constant" )
apifun_checktype ( "conint_normvar" , level, "level" , 3 , "constant" )
apifun_checktype ( "conint_normvar" , twosided, "twosided" , 4 , "boolean" )
//
// Check size
apifun_checkscalar ( "conint_normvar" , n , "n" , 1 )
apifun_checkscalar ( "conint_normvar" , v , "v" , 2 )
apifun_checkscalar ( "conint_normvar" , level , "level" , 3 )
apifun_checkscalar ( "conint_normvar" , twosided , "twosided" , 4 )
//
// Check content
apifun_checkgreq ( "conint_normvar" , n , "n" , 1 , 1. )
apifun_checkflint ( "conint_normvar" , n , "n" , 1 )
apifun_checkgreq ( "conint_normvar" , v , "v" , 2 , 0. )
apifun_checkrange ( "conint_normvar" , level , "level" , 3 , 0. , 0.5 )
//
if ( twosided ) then
q = level/2.
p = 1.-q
else
q = level
p = 1.-q
end
f = cdfchi("X",n-1,p,q)
low = (n-1)*v/f
if ( twosided ) then
p = level/2.
q = 1.-p
else
p = level
q = 1.-p
end
f = cdfchi("X",n-1,p,q)
up = (n-1)*v/f
endfunction
changelog.txt
1
2
3
4
5
changelog of the Conint Toolbox
conint (0.1)
* Initial import.
help/builder_help.sce
1
2
3
4
5
6
7
8
9
// ====================================================================
// Copyright INRIA 2008
// This file is released into the public domain
// ====================================================================
help_dir = get_absolute_file_path('builder_help.sce');
tbx_builder_help_lang("en_US", help_dir);
clear help_dir;
help/en_US/conint_bernoullipnumber.xml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from conint_bernoullipnumber.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="conint_bernoullipnumber" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>conint_bernoullipnumber</refname><refpurpose>Number of experiments for a Bernoulli variable.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
n = conint_bernoullipnumber ( len )
n = conint_bernoullipnumber ( len , level )
n = conint_bernoullipnumber ( len , level , pe )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>len :</term>
<listitem><para> a 1-by-1 matrix of doubles, positive, the length of the required C.I.</para></listitem></varlistentry>
<varlistentry><term>level :</term>
<listitem><para> a 1-by-1 matrix of doubles, the confidence level (default level = 1.-0.95=0.05). level is expected to be in the range [0.,0.5]</para></listitem></varlistentry>
<varlistentry><term>pe :</term>
<listitem><para> a 1-by-1 matrix of doubles, an estimate of the probability of success (default pe = 0.5)</para></listitem></varlistentry>
<varlistentry><term>n :</term>
<listitem><para> a 1-by-1 matrix of doubles, the number of trials</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
Returns the number of experiments which guarantees that the
length of the confidence interval for the
probability of success of a Bernoulli variable
is less than <literal>len</literal>.
</para>
<para>
In other words, if <literal>pe</literal> is the estimate of
the probability of success
with <literal>n</literal> trials, then
</para>
<para>
<screen>
[low, up] = conint_bernoullip ( pe , n )
</screen>
</para>
<para>
is so that
</para>
<para>
<literal>
abs(up-low) &lt;= len
</literal>
</para>
<para>
To compute <literal>n</literal>, we assume that the random variable
has Binomial distribution.
</para>
<para>
This function makes the assumption that the data
in <literal>x</literal> has a binomial distribution with
parameters <literal>p</literal> and
<literal>n</literal>.
</para>
<para>
The default value <literal>pe=0.5</literal>
guarantees that the number of experiments is large enough,
no matter the actual value of the probability <literal>p</literal>.
In other words, if <literal>pe</literal> is not provided
be the user, the value of the output variable <literal>n</literal>
is guaranteed, but might be much too large than required.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
// From Gilbert Saporta,
// Section 13.5.4 Intervalle de confiance
// pour une proportion p
// p.313
len = 0.05*2
level = 1.-0.90
n = conint_bernoullipnumber ( len , level )
// Therefore, it requires 271 experiments
// to make sure that the 90% C.I. for p
// has a length no greater than 0.1.
// Reproduce the table p 313.
mprintf("level = %9s %9s %9s %9s", ..
" ", "1.-0.90", "1.-0.95", "1.-0.98")
for len = 2*[0.01 0.02 0.05]
mprintf("len = %9s, n = ", string(len))
for level = 1.-[0.90 0.95 0.98]
n = conint_bernoullipnumber ( len , level );
mprintf("%9s ", string(n));
end
mprintf("\n");
end
// From Sheldon Ross,
// Example 7.5c
// Estimate the failure rate of computer chips.
len = 0.05
// Compute n for 95% C.I.
n = conint_bernoullipnumber ( len )
// Compute n for 99% C.I.
level = 1.-0.99
n = conint_bernoullipnumber ( len , level )
// Compute n for 99% C.I. and probability estimate
// From 30 sample computer chips, there are
// 26. chips which works.
pe = 26. / 30.
n = conint_bernoullipnumber ( len , level , pe )
// It may require 1227 experiments to get a
// 99 C.I. with length less than 0.05.
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Copyright (C) 2011 - Michael Baudin</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>http://en.wikipedia.org/wiki/Confidence_interval</para>
<para>"Introduction to probability and statistics for engineers and scientists", Sheldon Ross, Third Edition, 2004</para>
<para>"Probabilités, Analyse de Données et Statistique", Gilbert Saporta, 2nd Ed., 2006</para>
</refsection>
</refentry>
help/en_US/conint_normvar.xml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from conint_normvar.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="conint_normvar" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>conint_normvar</refname><refpurpose>C.I. of the variance of a normal variable.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
[ low, up ] = conint_normvar ( n , v )
[ low, up ] = conint_normvar ( n , v , level )
[ low, up ] = conint_normvar ( n , v , level , twosided )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>n :</term>
<listitem><para> a 1-by-1 matrix of doubles, positive, integer value, the number of samples</para></listitem></varlistentry>
<varlistentry><term>v :</term>
<listitem><para> a 1-by-1 matrix of doubles, positive, the variance estimate</para></listitem></varlistentry>
<varlistentry><term>level :</term>
<listitem><para> a 1-by-1 matrix of doubles, the confidence level (default level = 1.-0.95=0.05). level is expected to be in the range [0.,0.5]</para></listitem></varlistentry>
<varlistentry><term>twosided :</term>
<listitem><para> a 1-by-1 matrix of booleans, the side of the interval (default twosided = %t). If twosided is true, then low and up is a two-sided interval [low,up]. If twosided is false, then the intervals are [-inf,up] and [low,inf].</para></listitem></varlistentry>
<varlistentry><term>low :</term>
<listitem><para> a 1-by-1 matrix of doubles, the estimated lower bound</para></listitem></varlistentry>
<varlistentry><term>up :</term>
<listitem><para> a 1-by-1 matrix of doubles, the estimated upper bound</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
Computes a confidence interval of the variance of a normal variable, that is,
computes <literal>v</literal> as the variance of <literal>x</literal> and
computes confidence intervals for <literal>v</literal>.
</para>
<para>
If <literal>twosided</literal> is true, then
<literal>low</literal> and <literal>up</literal> are such that
the variance <literal>var(X)</literal> is such that
</para>
<para>
<literal>
P(low &lt; var &lt; up) = 1-level
</literal>
</para>
<para>
If <literal>twosided</literal> is false, then
<literal>low</literal> and <literal>up</literal> are such that
the variance <literal>var(X)</literal> is such that
</para>
<para>
<literal>
P(var(X) &lt; up) = 1-level
</literal>
</para>
<para>
and
</para>
<para>
<literal>
P(low &lt; var(X)) = 1-level
</literal>
</para>
<para>
This function makes the assumption that the
data in <literal>x</literal> has a normal distribution with
expectation <literal>mu(X)</literal> and
variance <literal>var(X)</literal>.
</para>
<para>
We use a Chi-square distribution to compute the confidence interval.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
// Sheldon Ross, Example 7.3h, page 252
// Compute a 90% confidence interval for the
// variance of a normal variable.
x = [0.123 0.124 0.126 0.120 0.130 ...
0.133 0.125 0.128 0.124 0.126]
n = size(x,"*")
v = variance(x)
[ low, up ] = conint_normvar ( n , v , 1.-0.90 )
// Then P(7.264D-06 <= var(X) <= 36.96D-06) = 0.90
// Give the confidence interval for the standard
// deviation
[sqrt(low) sqrt(v) sqrt(up)]
// Example from Gilbert Saporta,
// Section 13.5.3, "Variance d'une loi normale"
// Section 13.5.3.2, "m est inconnu"
// Get the C.I. for the variance of a normal
// variable with n=30 and var(X)=12.
x = [6;2;2;8;12;6;6;13;..
8;2;10;0;7;5;4;4;3;..
7;4;9;6;2;-3;5;7;2;4;5;2;2];
n = size(x,"*")
v = variance(x)
[ low, up ] = conint_normvar ( n , v , 1.-0.90 )
// Then P(8.177 < var(X) < 19.65) = 0.90
// The 90% C.I. is (8.177,19.65)
// Notice that Saporta gets (8.46,20.33).
// This is because Saporta uses n as the denominator
// in the variance estimate, while we use (n-1).
// Example from
// "Probability and Statistics for Engineers and Scientists",
// Walpole, Myers, Myers, Ye
// Example 9.17:
x = [46.4, 46.1, 45.8, 47.0, 46.1, ..
45.9, 45.8, 46.9, 45.2, 46.0];
n = size(x,"*")
v = variance(x)
[ low, up ] = conint_normvar ( n , v , 1.-0.95 )
// Then (0.135,0.953) is a 95% C.I. for the variance.
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Copyright (C) 2011 - Michael Baudin</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>http://en.wikipedia.org/wiki/Confidence_interval</para>
<para>"Introduction to probability and statistics for engineers and scientists", Sheldon Ross, Third Edition, 2004</para>
<para>"Probabilités, Analyse de Données et Statistique", Gilbert Saporta, 2nd Ed., 2006</para>
<para>// "Probability and Statistics for Engineers and Scientists", Ronald Walpole, Raymond Myers, Sharon Myers, Keying Ye, Eight Edition, 2006</para>
</refsection>
</refentry>
help/en_US/conint_normmunumber.xml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from conint_normmunumber.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="conint_normmunumber" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>conint_normmunumber</refname><refpurpose>Number of experiments for a Normal variable.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
n = conint_normmunumber ( len , v )
n = conint_normmunumber ( len , v , level )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>len :</term>
<listitem><para> a 1-by-1 matrix of doubles, positive, the length of the confidence interval.</para></listitem></varlistentry>
<varlistentry><term>v :</term>
<listitem><para> a 1-by-1 matrix of doubles, the variance</para></listitem></varlistentry>
<varlistentry><term>level :</term>
<listitem><para> a 1-by-1 matrix of doubles, the confidence level (default level = 1.-0.95=0.05). level is expected to be in the range [0.,0.5]</para></listitem></varlistentry>
<varlistentry><term>n :</term>
<listitem><para> a 1-by-1 matrix of doubles, integer value, the number of experiments</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
Returns the number of experiments which guarantees that the
length of the confidence interval is less than <literal>len</literal>.
This two-sided confidence interval is for the expectation of the
random variable and we make the assumption that the random
variable has a Normal distribution.
</para>
<para>
In other words, if <literal>x</literal> is a vector
with <literal>n</literal> entries, then
</para>
<para>
<screen>
[ m, low, up ] = conint_normmu ( x , level , v )
</screen>
</para>
<para>
is so that
</para>
<para>
<literal>
abs(up-low) &lt;= len
</literal>
</para>
<para>
To compute <literal>n</literal>, we assume that the random variable
has Normal distribution with known variance <literal>v</literal>.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
// Sheldon Ross, EXAMPLE 7.3d, p245
// Compute the number of experiments for a random variable
// with standard deviation 0.3 so that
// the expectation has a two-sided 95% confidence interval
// with length no greater than 0.2.
n = conint_normmunumber ( 0.2 , 0.3^2 , 1.-0.95 )
// Expected: n=35
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Copyright (C) 2011 - Michael Baudin</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>http://en.wikipedia.org/wiki/Confidence_interval</para>
<para>"Introduction to probability and statistics for engineers and scientists", Sheldon Ross, Third Edition, 2004</para>
</refsection>
</refentry>
help/en_US/build_help.sce
1
2
3
4
5
6
7
8
9
10
11
12
// ====================================================================
// Copyright INRIA 2008
// Allan CORNET
// Simon LIPP
// This file is released into the public domain
// ====================================================================
help_lang_dir = get_absolute_file_path('build_help.sce');
tbx_build_help(TOOLBOX_TITLE, help_lang_dir);
clear help_lang_dir;
help/en_US/update_help.sce
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
// Copyright (C) 2010-2011 - Michael Baudin
// Updates the .xml files by deleting existing files and
// creating them again from the .sci with help_from_sci.
//
cwd = get_absolute_file_path("update_help.sce");
mprintf("Working dir = %s\n",cwd);
//
// Generate the library help
mprintf("Updating conint\n");
helpdir = cwd;
funmat = [
"conint_normmu"
"conint_normmunumber"
"conint_normvar"
"conint_bernoullip"
"conint_bernoullipnumber"
];
macrosdir = fullfile(cwd,"..","..","macros");
demosdir = [];
modulename = "conint";
helptbx_helpupdate ( funmat , helpdir , macrosdir , demosdir , modulename , %t );
help/en_US/conint_bernoullip.xml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from conint_bernoullip.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="conint_bernoullip" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>conint_bernoullip</refname><refpurpose>C.I. of the probability of a Bernoulli variable.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
[low, up] = conint_bernoullip ( pe , n )
[low, up] = conint_bernoullip ( pe , n , level )
[low, up] = conint_bernoullip ( pe , n , level , twosided )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>pe :</term>
<listitem><para> a 1-by-1 matrix of doubles, an estimate of the probability</para></listitem></varlistentry>
<varlistentry><term>n :</term>
<listitem><para> a 1-by-1 matrix of doubles, the number of trials</para></listitem></varlistentry>
<varlistentry><term>level :</term>
<listitem><para> a 1-by-1 matrix of doubles, the confidence level (default level = 1.-0.95=0.05). level is expected to be in the range [0.,0.5]</para></listitem></varlistentry>
<varlistentry><term>twosided :</term>
<listitem><para> a 1-by-1 matrix of booleans, the side of the interval (default twosided = %t). If twosided is true, then low and up is a two-sided interval [low,up]. If twosided is false, then the intervals are [-inf,up] and [low,inf].</para></listitem></varlistentry>
<varlistentry><term>low :</term>
<listitem><para> a 1-by-1 matrix of doubles, the estimated lower bound</para></listitem></varlistentry>
<varlistentry><term>up :</term>
<listitem><para> a 1-by-1 matrix of doubles, the estimated upper bound</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
Computes a confidence interval of the probability
<literal>p</literal> of a Bernoulli variable.
This random variable is a "success" with probability
<literal>p</literal> and is a failure with probability
<literal>1-p</literal>.
We assume that we have performed <literal>n</literal> trials
of this experiment and we assume that <literal>pe</literal>
is the fraction of the <literal>n</literal> trials which are successes.
This function computes confidence intervals for <literal>p</literal>.
</para>
<para>
If <literal>twosided</literal> is true, then
<literal>low</literal> and <literal>up</literal> are such that
the probability <literal>p</literal> satisfies the equality
</para>
<para>
<literal>
P(low &lt; p &lt; up) = 1-level
</literal>
</para>
<para>
If <literal>twosided</literal> is false, then
<literal>low</literal> and <literal>up</literal> are such that
the probability <literal>var(X)</literal> satisfies the equality
</para>
<para>
<literal>
P(p &lt; up) = 1-level
</literal>
</para>
<para>
and
</para>
<para>
<literal>
P(low &lt; p) = 1-level
</literal>
</para>
<para>
This function makes the assumption that the data
in <literal>x</literal> has a binomial distribution with
parameters <literal>p</literal> and
<literal>n</literal>.
</para>
<para>
The confidence interval that we compute is approximate.
It is based on the Central Limit Theorem.
We use a Normal distribution to compute the confidence interval.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
// From Sheldon Ross, Example 7.5a, p.262
// We experimented with 100 transistors :
// 80 transistors works.
// Compute a 95% confidence interval for p,
// the probability that a transistor works.
pe = 80./100.;
n = 100;
[low, up] = conint_bernoullip ( pe , n , 1.-0.95 )
// Then (0.7216,0.8774) is a 95% C.I. for p.
// From Gilbert Saporta
// Section 13.5.4 Intervalle de confiance
// pour une proportion p
n = 400;
pe = 0.36;
[low, up] = conint_bernoullip ( pe , n , 1.-0.95 )
// Then (0.31,0.41) is a 95% C.I. for p.
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Copyright (C) 2011 - Michael Baudin</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>http://en.wikipedia.org/wiki/Confidence_interval</para>
<para>"Introduction to probability and statistics for engineers and scientists", Sheldon Ross, Third Edition, 2004</para>
<para>"Probabilités, Analyse de Données et Statistique", Gilbert Saporta, 2nd Ed., 2006</para>
</refsection>
</refentry>
help/en_US/conint_normmu.xml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
<?xml version="1.0" encoding="UTF-8"?>
<!--
*
* This help file was generated from conint_normmu.sci using help_from_sci().
*
-->
<refentry version="5.0-subset Scilab" xml:id="conint_normmu" xml:lang="en"
xmlns="http://docbook.org/ns/docbook"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns:ns3="http://www.w3.org/1999/xhtml"
xmlns:mml="http://www.w3.org/1998/Math/MathML"
xmlns:db="http://docbook.org/ns/docbook">
<refnamediv>
<refname>conint_normmu</refname><refpurpose>C.I. of the mean of a normal variable.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<title>Calling Sequence</title>
<synopsis>
[ low, up ] = conint_normmu ( n , me , v )
[ low, up ] = conint_normmu ( n , me , v , level )
[ low, up ] = conint_normmu ( n , me , v , level , twosided )
[ low, up ] = conint_normmu ( n , me , v , level , twosided , known )
</synopsis>
</refsynopsisdiv>
<refsection>
<title>Parameters</title>
<variablelist>
<varlistentry><term>n :</term>
<listitem><para> a n-by-1 or 1-by-n matrix of doubles, integer value, positive, the number of samples</para></listitem></varlistentry>
<varlistentry><term>me :</term>
<listitem><para> a 1-by-1 matrix of doubles, the mean estimate</para></listitem></varlistentry>
<varlistentry><term>level :</term>
<listitem><para> a 1-by-1 matrix of doubles, the confidence level (default level = 1.-0.95=0.05). level is expected to be in the range <literal>[0.,0.5]</literal></para></listitem></varlistentry>
<varlistentry><term>v :</term>
<listitem><para> a 1-by-1 matrix of doubles, the variance. By default, the variance is assumed to be exact (i.e. by default known=%t). If the variance <literal>v</literal> is estimated, then set the <literal>known</literal> option to %f.</para></listitem></varlistentry>
<varlistentry><term>twosided :</term>
<listitem><para> a 1-by-1 matrix of booleans, the side of the interval (default twosided = %t). If <literal>twosided</literal> is %t, then <literal>low</literal> and <literal>up</literal> make a two-sided interval <literal>[low,up]</literal>. If <literal>twosided</literal> is %f, then the intervals are <literal>[-inf,up]</literal> and <literal>[low,inf]</literal>.</para></listitem></varlistentry>
<varlistentry><term>known :</term>
<listitem><para> a 1-by-1 matrix of booleans, %t if the variance <literal>v</literal> is exact (default known = %t). If <literal>known</literal> is %t, then <literal>v</literal> is the exact variance. If known is %f, then <literal>v</literal> is the estimate of the variance.</para></listitem></varlistentry>
<varlistentry><term>low :</term>
<listitem><para> a 1-by-1 matrix of doubles, the estimated lower bound</para></listitem></varlistentry>
<varlistentry><term>up :</term>
<listitem><para> a 1-by-1 matrix of doubles, the estimated upper bound</para></listitem></varlistentry>
</variablelist>
</refsection>
<refsection>
<title>Description</title>
<para>
Computes a confidence interval of the mean of a normal variable, that is,
computes <literal>m</literal> as the mean of <literal>x</literal> and
computes confidence intervals for <literal>mu(X)</literal>.
</para>
<para>
If <literal>twosided</literal> is true, then
<literal>low</literal> and <literal>up</literal> are such that
the expectation <literal>mu(X)</literal> is such that
</para>
<para>
<literal>
P(low &lt; mu(X) &lt; up) = 1-level
</literal>
</para>
<para>
If <literal>twosided</literal> is false, then
<literal>low</literal> and <literal>up</literal> are such that
the expectation <literal>mu</literal> is such that
</para>
<para>
<literal>
P(mu(X) &lt; up) = 1-level
</literal>
</para>
<para>
and
</para>
<para>
<literal>
P(low &lt; mu(X)) = 1-level
</literal>
</para>
<para>
This function makes the assumption that the
data in x has a normal distribution with expectation <literal>mu(X)</literal> and
variance <literal>var(X)</literal>.
</para>
<para>
If the variance is exact (i.e. <literal>known=%t</literal>)
then we use a Normal distribution to compute the confidence interval.
</para>
<para>
If the variance is estimated (i.e. if <literal>known=%f</literal>),
then we use Student's T distribution to compute the confidence interval.
</para>
<para>
</para>
</refsection>
<refsection>
<title>Examples</title>
<programlisting role="example"><![CDATA[
// Sheldon Ross, Example 7.3a, page 241
x = [5, 8.5, 12, 15, 7, 9, 7.5, 6.5, 10.5]
n = size(x,"*")
me = mean(x)
// By default, we compute a two-sided 95% C.I.
// with exact variance:
v = 4.
[ low, up ] = conint_normmu ( n , me , v )
low_expected = 7.69
up_expected = 10.31
// Then P(7.69 <= mu <= 10.31) = 0.95
// Get one-sided 95% C.I.
// Example 7.3b, p. 242
low_expected = 7.903
up_expected = 10.097
twosided = %f
[ low, up ] = conint_normmu ( n , me , v , [] , twosided )
// Then : P(7.903 <= mu) = 0.95
// P(mu <= 10.097) = 0.95
// Example 7.3c, p. 244
// Two-sided 99% interval : (7.28, 10.72)
[ low, up ] = conint_normmu ( n , me , v , 1.-0.99 )
// One-sided 99% intervals : (7.447,inf), (−inf, 10.553)
[ low, up ] = conint_normmu ( n , me , v , 1.-0.99 , %f )
// EXAMPLE 7.3e, p. 247
// The variance is unknown: use Student's T random variable
v = variance(x)
[ low, up ] = conint_normmu ( n , me , v , [] , [] , %f )
// Then P(6.63<= mu <= 11.37) = 0.95
// EXAMPLE 7.3f, p.249
x = [54, 63, 58, 72, 49, 92, 70, 73, ..
69, 104, 48, 66, 80, 64, 77];
n = size(x,"*")
me = mean(x)
v = variance(x)
// Two-sided :
[ low, up ] = conint_normmu ( n , me , v , 1.-0.95 , [] , %f )
// Then P(60.865<= mu <=77.6683) = 0.95
// One-sided :
[ low, up ] = conint_normmu ( n , me , v , 1.-0.95 , %f , %f )
// Then : P (mu <= 76.16) = 0.95
// Then : P (62.368 <= mu) = 0.95
// EXAMPLE 7.3g, Monte-Carlo simulation
U = grand(100,1,"def");
x = sqrt(1-U.^2);
n = size(x,"*")
me = mean(x)
v = variance(x)
[ low, up ] = conint_normmu ( n , me , v , [] , [] , %f )
// Exact integral is :
%pi/4
]]></programlisting>
</refsection>
<refsection>
<title>Authors</title>
<simplelist type="vert">
<member>Copyright (C) 2011 - Michael Baudin</member>
</simplelist>
</refsection>
<refsection>
<title>Bibliography</title>
<para>http://en.wikipedia.org/wiki/Confidence_interval</para>
<para>"Introduction to probability and statistics for engineers and scientists", Sheldon Ross, Third Edition, 2004</para>
</refsection>
</refentry>
etc/conint.start
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) 2008 - INRIA - Michael Baudin
// Copyright (C) 2008 - INRIA - Allan CORNET
// Copyright (C) 2010 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function conintlib = loadconintlib ()
TOOLBOX_NAME = "conint"
TOOLBOX_TITLE = "Conint"
mprintf("Start %s\n",TOOLBOX_TITLE);
etc_tlbx = get_absolute_file_path(TOOLBOX_NAME+".start");
etc_tlbx = getshortpathname(etc_tlbx);
root_tlbx = strncpy( etc_tlbx, length(etc_tlbx)-length("\etc\") );
//Load functions library
// =============================================================================
if ( %t ) then
mprintf("\tLoad macros\n");
pathmacros = pathconvert( root_tlbx ) + "macros" + filesep();
conintlib = lib(pathmacros);
end
// load gateways
// =============================================================================
if ( %f ) then
mprintf("\tLoad gateways\n");
ilib_verbose(0);
exec( pathconvert(root_tlbx+"/sci_gateway/loader_gateway.sce",%f));
end
// Load and add help chapter
// =============================================================================
if ( %t ) then
if or(getscilabmode() == ["NW";"STD"]) then
mprintf("\tLoad help\n");
path_addchapter = pathconvert(root_tlbx+"/jar");
if ( isdir(path_addchapter) <> [] ) then
add_help_chapter(TOOLBOX_TITLE, path_addchapter, %F);
end
end
end
// add demos
// =============================================================================
if ( %F ) then
if or(getscilabmode() == ["NW";"STD"]) then
mprintf("\tLoad demos\n");
demoscript = TOOLBOX_NAME + ".dem.gateway.sce"
pathdemos = pathconvert(fullfile(root_tlbx,"demos",demoscript),%f,%t);
add_demo(TOOLBOX_TITLE,pathdemos);
end
end
// ====================================================================
// A Welcome message.
//mprintf("\tType ""help intprb_overview"" for quick start.\n");
//mprintf("\tType ""demo_gui()"" and search for "+TOOLBOX_TITLE+" for Demonstrations.\n");
endfunction
if ( isdef("conintlib") ) then
warning("Library is already loaded (""ulink(); clear conintlib;"" to unload.)");
return;
end
conintlib = loadconintlib ();
clear loadconintlib;
etc/conint.quit
1
readme.txt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
Conint toolbox
Purpose
----
The goal of this toolbox is to provide functions to
compute confidence intervals for statistical estimators.
The toolbox is based on macros.
Features
--------
The following is a list of the current accsum functions :
* accsum_dcs : A Doubly Self Compensated Sum algorithm
* accsum_scs : A Self Compensated Sum algorithm
* accsum_compsum : The compensated sum of a matrix.
* accsum_dblcompsum : The doubly compensated sum of a matrix.
* accsum_fasttwosum : The fast2sum sum of a and b.
* accsum_orderdynamic : Returns the sum with a dynamic re-ordering.
* accsum_straight : The straightforward sum of a matrix.
* accsum_twosum : The twosum sum of a and b.
* accsum_fcompsum : The compensated sum of a matrix.
* accsum_fdcs : A Doubly Self Compensated Sum algorithm
* accsum_fscs : A Self Compensated Sum algorithm
and support functions:
* accsum_getpath : Returns the path to the current module.
* accsum_order : Re-order the matrix.
* accsum_priestx : A difficult example for SCS by Priest.
* accsum_shuffle : Randomly shuffles the input.
* accsum_sumcond : Condition number of the sum function.
* accsum_wilkinson : A test vector by Wilkinson.
* accsum_higham : Returns an example designed by Higham.
The accsum_fcompsum, accsum_fdcs and accsum_fscs functions
are based on compiled source code and are faster than the other.
Dependencies
----
* This modules depends on the assert module.
* This modules depends on the helptbx module.
* This modules depends on the apifun module.
Author
------
Copyright (C) 2011 - Michael Baudin
Forge
-----
http://forge.scilab.org/index.php/p/accsum/
ATOMS
-----
Not yet
Licence
----
This toolbox is released under the CeCILL_V2 licence :
http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
Bibliography
------------
* "Stability and numerical accuracy of algorithms", Nicolas Higham
* "Handbook of Floating Point Computations", Muller et al
* https://hpcrd.lbl.gov/SCG/ocean/NRS/ECMWF/img14.htm
* https://hpcrd.lbl.gov/SCG/ocean/NRS/SCSsum.F
* "On properties of floating point arithmetics: numerical stability and the cost of accurate computations", Douglas Priest, 1992
* "Using Accurate Arithmetics to Improve Numerical Reproducibility and Stability in Parallel Applications". Yun He and Chris H.Q. Ding. Journal of Supercomputing, Vol.18, Issue 3, 259-277, March 2001. Also Proceedings of International Conference on Supercomputing (ICS'00), May 2000, 225-234.
builder.sce
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
//
// builder.sce --
// Builder for the Accsum Scilab Toolbox
//
mode(-1);
lines(0);
//setenv("DEBUG_SCILAB_DYNAMIC_LINK","YES")
try
getversion("scilab");
catch
error(gettext("Scilab 5.0 or more is required."));
end;
// ====================================================================
if ~with_module("development_tools") then
error(msprintf(gettext("%s module not installed."),"development_tools"));
end
// ====================================================================
TOOLBOX_NAME = "conint";
TOOLBOX_TITLE = "Conint";
// ====================================================================
toolbox_dir = get_absolute_file_path("builder.sce");
//tbx_builder_src(toolbox_dir);
//tbx_builder_gateway(toolbox_dir);
tbx_builder_macros(toolbox_dir);
tbx_builder_help(toolbox_dir);
tbx_build_loader(TOOLBOX_NAME, toolbox_dir);
tbx_build_cleaner(TOOLBOX_NAME, toolbox_dir);
clear toolbox_dir TOOLBOX_NAME TOOLBOX_TITLE;
// ====================================================================
license.txt
1
2
3
4
5
Accsum Scilab Toolbox
This toolbox is released under the terms of the CeCILL license :
http://www.cecill.info/index.en.html

Archive Download the corresponding diff file

Revision: 1