* Bug #7858 fixed - Statistics: a priori mean for variance and variancef 09/13009/8
Paul BIGNIER [Fri, 25 Oct 2013 16:01:26 +0000 (18:01 +0200)]
Completing commit https://codereview.scilab.org/#/c/12855/ ,
variance and variancef can now take the a priori mean as input, and it can be a scalar.

Change-Id: I4ea4f041dfd5d27cbc31964ca32f6aa1e8c6ed7d

14 files changed:
SEP/INDEX
SEP/SEP_112_variance_variancef.odt
scilab/CHANGES_5.5.X
scilab/modules/statistics/help/en_US/descriptive_statistics/variance.xml
scilab/modules/statistics/help/en_US/descriptive_statistics/variancef.xml
scilab/modules/statistics/help/fr_FR/descriptive_statistics/variance.xml
scilab/modules/statistics/macros/variance.sci
scilab/modules/statistics/macros/variancef.sci
scilab/modules/statistics/tests/nonreg_tests/bug_7858.dia.ref [new file with mode: 0644]
scilab/modules/statistics/tests/nonreg_tests/bug_7858.tst [new file with mode: 0644]
scilab/modules/statistics/tests/unit_tests/variance.dia.ref
scilab/modules/statistics/tests/unit_tests/variance.tst
scilab/modules/statistics/tests/unit_tests/variancef.dia.ref
scilab/modules/statistics/tests/unit_tests/variancef.tst

index d018488..e6b5b9c 100644 (file)
--- a/SEP/INDEX
+++ b/SEP/INDEX
@@ -107,5 +107,5 @@ SEP #108: 3-D arrows with xarrows
 SEP #109: New input argument for neldermead_search(), to control warnings.
 SEP #110: /* RESERVED */
 SEP #111: /* RESERVED */
-SEP #112: variance and variancef can now return the mean of the input in a new output argument.
+SEP #112: variance and variancef can now return the mean of the input in a new output argument and take the a priori mean as input.
 SEP #113: resize_matrix hypermatrices support.
index ad7d13f..c5865e7 100644 (file)
Binary files a/SEP/SEP_112_variance_variancef.odt and b/SEP/SEP_112_variance_variancef.odt differ
index e7af7f8..e76d916 100644 (file)
@@ -35,7 +35,7 @@ Scilab Bug Fixes
                     "dimension", "minbounds" and "maxbounds" fields.
 
 * Bug #7858 fixed - variance and variancef can now return the mean of the input
-                    in a new output argument.
+                    in a new output argument and take the a priori mean as input.
 
 * Bug #7879 fixed - string now accepts plist type, and printing a plist displays that string.
 
index 248d139..2260c2b 100644 (file)
@@ -3,11 +3,11 @@
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) 2000 - INRIA - Carlos Klimann
  * Copyright (C) 2013 - Samuel GOUGEON
- * 
+ *
  * 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    
+ * are also available at
  * http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
  *
  -->
     </refnamediv>
     <refsynopsisdiv>
         <title>Calling Sequence</title>
-        <synopsis>s = variance(x [,orien [,w]])
-            [s,m] = variance(x,'r') or s = variance(x,1)
-            [s,m] = variance(x,'c') or s = variance(x,2)
-            [s,m] = variance(x,'*',1)
-            [s,m] = variance(x,'r',1)
-            [s,m] = variance(x,'c',1)
+        <synopsis>
+            s = variance(x [,orien [,m]])
+            [s, m] = variance(x, 'r') or s = variance(x, 1)
+            [s, m] = variance(x, 'c') or s = variance(x, 2)
+            [s, m] = variance(x, '*', m)
+            [s, m] = variance(x, 'r', m)
+            [s, m] = variance(x, 'c', m)
         </synopsis>
     </refsynopsisdiv>
     <refsection>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>w</term>
+                <term>m</term>
                 <listitem>
-                    <para>w : type of normalization to use. Valid values are 
-                        <itemizedlist>
-                            <listitem>0 : provides the best unbiased estimator of the variance when the mean is not known a priori. This is the default.
-                                The sum is normalized by 1/(nE-1), where nE is the number of summed elements 
-                                (the number of rows nR if "r" is used; of columns nC if "c" is used; or by default nC*nR)
-                            </listitem>
-                            <listitem>1 : provides the second moment around the mean. The sum is normalized by 1/nE.
-                                This estimation is unbiased only if the mean is known a priori.
-                            </listitem>
-                        </itemizedlist>
+                    <para>
+                        m: the a priori mean. If it is passed to <literal>variance</literal>, then the sum is normalized by
+                        <literal>n</literal>, the function returns the second moment around the mean.
+                        Otherwise, it is normalized by <literal>n-1</literal> and <literal>variance</literal> provides the best unbiased estimator of the variance.
                     </para>
                 </listitem>
             </varlistentry>
         <para>
             The second output argument <literal>m</literal> is the mean of the input, with respect to the <literal>orien</literal> parameter.
         </para>
+        <para>
+            <note>
+                If <code>m=%nan</code>, then the biased variance is returned:
+                that is, the chosen a priori mean is that of <varname>x</varname> and the denominator is <code>n</code>.
+            </note>
+        </para>
     </refsection>
     <refsection>
         <title>Examples</title>
-        <programlisting role="example"><![CDATA[ 
-x = [ 0.2113249 0.0002211 0.6653811;0.7560439 0.4453586 0.6283918 ]
+        <programlisting role="example"><![CDATA[
+x = [ 0.2113249 0.0002211 0.6653811; 0.7560439 0.4453586 0.6283918 ]
 s = variance(x)
 s = variance(x, "r")
 s = variance(x, "c")
+// If the a priori mean is known
+m_r = mean(x, "r");
+m_c = mean(x, "c");
+s = variance(x, "r", m_r)
+s = variance(x, "c", m_c)
 
 // with complex numbers
-x = rand(4,3) + rand(4,3)*%i
+x = rand(4, 3) + rand(4, 3)*%i
 s = variance(x)
-s = variance(x, "*", 1)        // should be smaller than the previous one
 s = variance(x, "r")
-s = variance(x, "r", 1)
 s = variance(x, "c")
 
 // with an hypermatrix
-x = rand(3,2,2)
+x = rand(3, 2, 2)
 s = variance(x)
-s = variance(x, "*", 1)        // should be smaller than the previous one
 // s = variance(x, "r")  // is not supported for an hypermatrix
 // s = variance(x, "c")  // is not supported for an hypermatrix
  ]]></programlisting>
@@ -135,6 +138,9 @@ s = variance(x, "*", 1)     // should be smaller than the previous one
         <title>History</title>
         <revhistory>
             <revision>
+                <revnumber>5.5.0</revnumber>
+                <revdescription>variance(x,orien,***) introduced. Variance can now take the a priori mean.
+                </revdescription>
                 <revnumber>5.4.1</revnumber>
                 <revdescription>variance(complexes) fixed. variance(x,"*",1) introduced. Vectorization of the code for directional usages variance(x,"r"|"c"). Full revision of the help page
                 </revdescription>
index 22f208a..1daf7e8 100644 (file)
@@ -2,11 +2,11 @@
 <!--
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) 2000 - INRIA - Carlos Klimann
- * 
+ *
  * 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    
+ * are also available at
  * http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
  *
  -->
     </refnamediv>
     <refsynopsisdiv>
         <title>Calling Sequence</title>
-        <synopsis>s=variancef(x,fre)
-            [s,m]=variancef(x,fre,'r') or s=variancef(x,fre,1)
-            [s,m]=variancef(x,fre,'c') or s=variancef(x,fre,2)
+        <synopsis>s = variancef(x, fre [,orien [,m]])
+            [s, m] = variancef(x, fre, 'r') or s = variancef(x, fre, 1)
+            [s, m] = variancef(x, fre, 'c') or s = variancef(x, fre, 2)
+            [s, m] = variancef(x, fre, 'r', m) or s = variancef(x, fre, m)
+            [s, m] = variancef(x, fre, 'c', m) or s = variancef(x, fre, m)
         </synopsis>
     </refsynopsisdiv>
     <refsection>
                     <para>real or complex vector or matrix</para>
                 </listitem>
             </varlistentry>
+            <varlistentry>
+                <term>m</term>
+                <listitem>
+                    <para>real or complex vector or matrix. The a priori weighted mean, to simplify computations</para>
+                </listitem>
+            </varlistentry>
         </variablelist>
     </refsection>
     <refsection>
             The second output argument <literal>m</literal> is the mean of the input,
             with respect to the <literal>orien</literal> parameter.
         </para>
+        <para>
+            If the weighted mean of <literal>x</literal> is known, then it can be passed to <literal>variancef</literal> to simplify its computations.
+        </para>
+        <para>
+            <note>
+                If <code>m=%nan</code>, then the biased variance is returned:
+                that is, the chosen a priori mean is that of <varname>x</varname> and the denominator is <code>size(x,"*")</code>.
+            </note>
+        </para>
     </refsection>
     <refsection>
         <title>Examples</title>
-        <programlisting role="example"><![CDATA[ 
+        <programlisting role="example"><![CDATA[
 x = [0.2113249 0.0002211 0.6653811; 0.7560439 0.9546254 0.6283918]
 fre = [1 2 3; 3 4 3]
 [s, m] = variancef(x, fre)
 [s, m] = variancef(x, fre, "r")
 [s, m] = variancef(x, fre, "c")
+
+// If the weighted mean is known
+m_r = meanf(x, fre, "r");
+m_c = meanf(x, fre, "c");
+[s, m] = variancef(x, fre, "r", mean_r)
+[s, m] = variancef(x, fre, "c", mean_c)
  ]]></programlisting>
     </refsection>
     <refsection role="see also">
@@ -97,4 +120,14 @@ fre = [1 2 3; 3 4 3]
             Wonacott, T.H. &amp; Wonacott, R.J.; Introductory Statistics, fifth edition, J.Wiley &amp; Sons, 1990.
         </para>
     </refsection>
+    <refsection>
+        <title>History</title>
+        <revhistory>
+            <revision>
+                <revnumber>5.5.0</revnumber>
+                <revdescription>variancef(x,fre,orien,***) introduced. Variancef can now take the a priori mean if it is known, not to compute it internally.
+                </revdescription>
+            </revision>
+        </revhistory>
+    </refsection>
 </refentry>
index 883f704..da37752 100644 (file)
@@ -3,11 +3,11 @@
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) 2000 - INRIA - Carlos Klimann
  * Copyright (C) 2013 - Samuel GOUGEON
- * 
+ *
  * 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    
+ * are also available at
  * http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
  *
  -->
     </refnamediv>
     <refsynopsisdiv>
         <title>Séquence d'appel</title>
-        <synopsis>s = variance(x [,orien [,w]])
-            s = variance(x,'r') or s = variance(x,1)
-            s = variance(x,'c') or s = variance(x,2)
-            s = variance(x,'*',1)
-            s = variance(x,'r',1)
-            s = variance(x,'c',1)
+        <synopsis>s = variance(x [,orien [,m]])
+            s = variance(x, 'r') or s = variance(x, 1)
+            s = variance(x, 'c') or s = variance(x, 2)
+            s = variance(x, '*', m)
+            s = variance(x, 'r', m)
+            s = variance(x, 'c', m)
         </synopsis>
     </refsynopsisdiv>
     <refsection>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>w</term>
+                <term>m</term>
                 <listitem>
-                    <para>w : indique le type de normalisation de la variance : 
-                        <itemizedlist>
-                            <listitem>
-                                0 : la somme des écarts quadratiques est divisée par (nE-1), où nE est le nombre d'écarts sommés (mode de calcul utilisé a priori). La variance ainsi estimée est non biaisée (l'espérance de la variable aléatoire n'est pas connue a priori ; elle est estimée par la moyenne des valeurs de x). nE est le nombre nR de lignes de <literal>x</literal> si "r" est utilisée ; ou le nombre nC de colonnes si "c" est utilisée ; ou le nombre total nR*nC d'éléments, sinon.
-                            </listitem>
-                            <listitem>1 : la somme des écarts quadratiques est divisée par nE. La variance résultante est non biaisée uniquement si l'espérance de la variable aléatoire est connue (Scilab ne permet actuellement pas de l'indiquer en paramètre).
-                            </listitem>
-                        </itemizedlist>
+                    <para>
+                        m: la moyenne a priori. Si cet argument est présent, alors la somme est normalisée par
+                        <literal>n</literal>.
+                        Sinon, elle est normalisée par <literal>n-1</literal> et <literal>variance</literal> retourne le mielleur estimateur non-biaisé de la variance.
                     </para>
                 </listitem>
             </varlistentry>
             retourne  dans le scalaire <literal>s</literal> la variance de tous les éléments de <literal>x</literal>.
         </para>
         <para>
-            <literal>s = variance(x,'c')</literal> (ou indifféremment <literal>s = variance(x, 2)</literal>) calcule la variance de chaque ligne. 
+            <literal>s = variance(x,'c')</literal> (ou indifféremment <literal>s = variance(x, 2)</literal>) calcule la variance de chaque ligne.
             Le vecteur colonne <literal>s</literal> est retourné, avec <literal>s(j) = variance(x(j,:),..)</literal>.
         </para>
         <para>
-            <literal>s = variance(x,'r')</literal> (ou indifféremment <literal>s = variance(x,1)</literal>) calcule la variance de chaque colonne. 
+            <literal>s = variance(x,'r')</literal> (ou indifféremment <literal>s = variance(x,1)</literal>) calcule la variance de chaque colonne.
             Le vecteur ligne <literal>s</literal> est retourné, avec <literal>s(i) = variance(x(:,i),..)</literal>.
         </para>
+        <para>
+            <note>
+                Si <code>m=%nan</code>, alors la variance biaisée est retournée;
+                c'est-à-dire que la moyenne a priori choisie est celle de <varname>x</varname> et le dénominateur est <code>n</code>.
+            </note>
+        </para>
     </refsection>
     <refsection>
         <title>Exemples</title>
-        <programlisting role="example"><![CDATA[ 
+        <programlisting role="example"><![CDATA[
 x = [ 0.2113249 0.0002211 0.6653811;0.7560439 0.4453586 0.6283918 ]
 s = variance(x)
 s = variance(x, "r")
 s = variance(x, "c")
 
-// with complex numbers
+// Si la moyenne est connue a priori
+m_r = mean(x, "r");
+m_c = mean(x, "c");
+s = variance(x, "r", m_r)
+s = variance(x, "c", m_c)
+
+// avec des nombres complexes
 x = rand(4,3) + rand(4,3)*%i
 s = variance(x)
-s = variance(x, "*", 1)        // doit être inférieure à la précédente valeur
 s = variance(x, "r")
-s = variance(x, "r", 1)
 s = variance(x, "c")
 
-// with an hypermatrix
-x = rand(3,2,2)
+// avec une hypermatrice
+x = rand(3, 2, 2)
 s = variance(x)
-s = variance(x, "*", 1)        // doit être inférieure à la précédente valeur
 // s = variance(x, "r")  // utilisation non admise pour une hypermatrice
 // s = variance(x, "c")  // utilisation non admise pour une hypermatrice
  ]]></programlisting>
index 6a189a1..3b4ae37 100644 (file)
@@ -9,7 +9,7 @@
 // http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
 //
 
-function [s,m] = variance(x, orien, w)
+function [s, m] = variance(x, orien, m)
     //
     //This function computes  the  variance  of  the values of  a  vector or
     //matrix x.
@@ -25,6 +25,8 @@ function [s,m] = variance(x, orien, w)
     //deviation.   It  returns in  each  entry  of the   column vector y the
     //variance of each row of x.
     //
+    //The input argument m represents the a priori mean. If it is present, then the sum is
+    //divided by n. Otherwise ("sample variance"), it is divided by n-1.
     //
 
     // Checking and normalizing input arguments:
@@ -44,6 +46,28 @@ function [s,m] = variance(x, orien, w)
         orien = "*"
     end
 
+    if rhs==3 then
+        if typeof(m)~="constant" then
+            tmp = gettext("%s: Wrong value of m : a priori mean expected.\n")
+            error(msprintf(tmp, "variance"))
+        elseif orien=="*" then
+            if ~isscalar(m) then
+                tmp = gettext("%s: Wrong value of m : a priori mean expected.\n")
+                error(msprintf(tmp, "variance"))
+            end
+        elseif orien=="r" | orien==1 then
+            if size(m)~=[1 size(x,"c")] & ~isscalar(m) then
+                tmp = gettext("%s: Wrong value of m : a priori mean expected.\n")
+                error(msprintf(tmp, "variance"))
+            end
+        elseif orien=="c" | orien==2 then
+            if size(m)~=[size(x,"r") 1] & ~isscalar(m) then
+                tmp = gettext("%s: Wrong value of m : a priori mean expected.\n")
+                error(msprintf(tmp, "variance"))
+            end
+        end
+    end
+
     transposed = %f // to refer and process as in "r", we priorly transpose any "c" request
     if orien=="r" | orien==1 | orien=="c" | orien==2 | orien=="*"
         if orien=="c" | orien==2 then
@@ -56,21 +80,39 @@ function [s,m] = variance(x, orien, w)
         error(msprintf(tmp, "variance", "c", "r", 1, 2))
     end
 
-    if rhs<3 then
-        w = 0
-    elseif typeof(w)~="constant" | and(w~=[0 1]) then
-        tmp = gettext("%s: Wrong value of w : %s; %d or %d expected.\n")
-        error(msprintf(tmp, "variance", string(w), 0, 1))
-    end
-
     // Calculations
     // ------------
-    d = size(x, orien) - 1 + w    // denominator
+
+    d = size(x, orien) - 1 + exists("m","local") // Denominator. If m is given, then the a priori mean is known and we divide by size(n,orien)
+    if rhs == 3 & isnan(m) then
+        // This will compute the "biased variance": the denominator is size(x,orien) but the a priori mean is not considered as provided.
+        rhs = 2
+    end
     if orien=="*" then
-        m = mean(x)
+        if rhs < 3 then
+            m = mean(x)
+        end
     else
-        m = mean(x, orien).*.ones(size(x,1),1)
+        if rhs < 3 then
+            m = mean(x, orien).*.ones(size(x,1),1)
+        else
+            if isscalar(m) then
+                if or(m==[0 1]) then
+                    tmp = _("%s: The use of input argument ''%s'' is now obsolete, please use ''%s'' instead.\n")
+                    warning(msprintf(tmp, "variance", "w", "m"))
+                end
+                // If m is a scalar, extend it to the size of x.
+                // If lhs==1, we do not need to perform this operation, because in the following 'x - m', m can be a scalar
+                m = m*ones(x)
+            else
+                if transposed == %t then
+                    m = m';
+                end
+                m = m.*.ones(size(x,1),1)
+            end
+        end
     end
+
     s = sum(abs(x - m).^2, orien) / d
 
     m = m(1, :);
index 931b9b0..7d30847 100644 (file)
@@ -9,7 +9,7 @@
 // http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
 //
 
-function [s,m]=variancef(x,fre,orien)
+function [s,m]=variancef(x,fre,orien,m)
     //
     //This function  computes  the variance  of the values  of   a vector or
     //matrix x, each  of  them  counted with  a  frequency signaled   by the
@@ -30,11 +30,14 @@ function [s,m]=variancef(x,fre,orien)
     //each row of  x, each value counted with  the multiplicity indicated by
     //the corresponding value of fre.
     //
+    //The input argument m represents the a priori mean. If it is present, then the sum is
+    //divided by n. Otherwise ("sample variance"), it is divided by n-1.
+    //
     //
     if x==[] then s=%nan, return, end
     [lhs,rhs]=argn(0)
-    if rhs<2|rhs>3 then
-        error(msprintf(gettext("%s: Wrong number of input arguments: %d to %d expected.\n"),"variancef",2,3)),
+    if rhs<2|rhs>4 then
+        error(msprintf(gettext("%s: Wrong number of input arguments: %d to %d expected.\n"),"variancef",2,4)),
     end
     if x==[]|fre==[]|fre==0, s=%nan;return,end
     if rhs==2 then
@@ -44,23 +47,71 @@ function [s,m]=variancef(x,fre,orien)
         s=(sum(((x-m).^2).*fre))/(sumfre-1),
         return,
     end
+    biased = %f
+    if rhs==4 then
+        if typeof(m)~="constant" then
+            tmp = gettext("%s: Wrong value of m : a priori mean expected.\n")
+            error(msprintf(tmp, "variance", ))
+        elseif orien=="*" then
+            if ~isscalar(m) then
+                tmp = gettext("%s: Wrong value of m : a priori mean expected.\n")
+                error(msprintf(tmp, "variance", ))
+            end
+        elseif orien=="r" | orien==1 then
+            if size(m)~=[1 size(x,"c")] & ~isscalar(m) then
+                tmp = gettext("%s: Wrong value of m : a priori mean expected.\n")
+                error(msprintf(tmp, "variance", ))
+            end
+        elseif orien=="c" | orien==2 then
+            if size(m)~=[size(x,"r") 1] & ~isscalar(m) then
+                tmp = gettext("%s: Wrong value of m : a priori mean expected.\n")
+                error(msprintf(tmp, "variance", ))
+            end
+        end
+        if isnan(m) then
+            biased = %t; // Compute the biased variance
+        end
+    end
     if orien=="*",
         sumfre=sum(fre)
         if sumfre <= 1 then error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be > %d.\n"),"variancef", 2, 1)),end
-        m = meanf(x,fre)
-        s=(sum(((x-m).^2).*fre))/(sumfre-1),
+        if rhs<4 then
+            m = meanf(x,fre)
+            s=(sum(((x-m).^2).*fre))/(sumfre-1),
+        elseif biased == %t
+            m = meanf(x,fre)
+            s=(sum(((x-m).^2).*fre))/sumfre,
+        else
+            s=(sum(((x-m).^2).*fre))/sumfre,
+        end
     elseif orien=="r"|orien==1,
         sumfre=sum(fre,"r")
         if or(sumfre==0) then error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be > %d.\n"),"variancef",2,1)),end
-        m = meanf(x,fre,"r")
+        if rhs<4 | biased == %t then
+            m = meanf(x,fre,"r")
+        elseif isscalar(m) then
+            m = m*ones(1, size(x,"c"));
+        end
         m2 = ones(size(x,"r"),1)*m
-        s=(sum(((x-m2).^2).*fre))./(sumfre-1)
+        if rhs<4 then
+            s=(sum(((x-m2).^2).*fre))./(sumfre-1)
+        else
+            s=(sum(((x-m2).^2).*fre))./sumfre
+        end
     elseif orien=="c"|orien==2,
         sumfre=sum(fre,"c")
         if or(sumfre==0) then error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be > %d.\n"),"variancef",2,1)),end
-        m = meanf(x,fre,"c")
+        if rhs<4 | biased == %t then
+            m = meanf(x,fre,"c")
+        elseif isscalar(m) then
+            m = m*ones(size(x,"r"), 1);
+        end
         m2 = m*ones(1,size(x,"c"))
-        s=(sum((x-m2).^2,"c"))./(sumfre-1)
+        if rhs<4 then
+            s=(sum((x-m2).^2,"c"))./(sumfre-1)
+        else
+            s=(sum((x-m2).^2,"c"))./sumfre
+        end
     else error(msprintf(gettext("%s: Wrong value for input argument #%d: ''%s'', ''%s'', ''%s'', %d or %d expected.\n"),"variancef",3,"*","c","r",1,2))
     end
 
diff --git a/scilab/modules/statistics/tests/nonreg_tests/bug_7858.dia.ref b/scilab/modules/statistics/tests/nonreg_tests/bug_7858.dia.ref
new file mode 100644 (file)
index 0000000..e3a1a1e
--- /dev/null
@@ -0,0 +1,61 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// <-- Non-regression test for bug 7858 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=7858
+//
+// <-- Short Description -->
+// variance and variancef can take now take the a priori mean as input, and it can be a scalar.
+// =============================================================================
+//--------------------------------------------------------------
+// Variance
+x = [
+0.9    0.7
+0.1    0.1
+0.5    0.4 ];
+orien = "r";
+expectedM = [1 1];
+expectedV = [1.07 1.26]/3;
+// Voluntarily passing a bad a priori mean, to check that we get the same results with scalar
+[Variance, Mean] = variance( x, orien, ones(mean(x,orien)) );
+assert_checkalmostequal([Variance Mean], [expectedV expectedM]);
+// With a scalar
+[Variance, Mean] = variance( x, orien, 1 );
+WARNING: variance: The use of input argument 'w' is now obsolete, please use 'm' instead.
+assert_checkalmostequal([Variance Mean], [expectedV expectedM]);
+orien = "c";
+expectedM = [1; 1; 1];
+expectedV = [0.05; 0.81; 0.305];
+// Voluntarily passing a bad a priori mean, to check that we get the same results with scalar
+[Variance, Mean] = variance( x, orien, ones(mean(x,orien)) );
+assert_checkalmostequal([Variance Mean], [expectedV expectedM]);
+// With a scalar
+[Variance, Mean] = variance( x, orien, 1 );
+WARNING: variance: The use of input argument 'w' is now obsolete, please use 'm' instead.
+assert_checkalmostequal([Variance Mean], [expectedV expectedM]);
+//--------------------------------------------------------------
+// Variancef
+x = [0.2113249 0.0002211 0.6653811; 0.7560439 0.9546254 0.6283918];
+fre = [1 2 3; 3 4 3];
+orien = "r";
+refM = [1 1 1];
+refV = [0.889522663062 0.593015108708 0.593015108708];
+[v, m] = variancef( x, fre, orien, ones(meanf(x,fre,orien)) );
+assert_checkalmostequal([v m], [refV refM]);
+[v, m] = variancef( x, fre, orien, 1 );
+assert_checkalmostequal([v m], [refV refM]);
+orien = "c";
+refM = [1; 1];
+refV = [0.288922678414; 0.019966608736];
+[v, m] = variancef( x, fre, orien, ones(meanf(x,fre,orien)) );
+assert_checkalmostequal([v m], [refV refM]);
+[v, m] = variancef( x, fre, orien, 1 );
+assert_checkalmostequal([v m], [refV refM]);
diff --git a/scilab/modules/statistics/tests/nonreg_tests/bug_7858.tst b/scilab/modules/statistics/tests/nonreg_tests/bug_7858.tst
new file mode 100644 (file)
index 0000000..4da83e2
--- /dev/null
@@ -0,0 +1,71 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// <-- Non-regression test for bug 7858 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=7858
+//
+// <-- Short Description -->
+// variance and variancef can take now take the a priori mean as input, and it can be a scalar.
+// =============================================================================
+
+//--------------------------------------------------------------
+// Variance
+
+x = [
+0.9    0.7
+0.1    0.1
+0.5    0.4 ];
+orien = "r";
+expectedM = [1 1];
+expectedV = [1.07 1.26]/3;
+
+// Voluntarily passing a bad a priori mean, to check that we get the same results with scalar
+[Variance, Mean] = variance( x, orien, ones(mean(x,orien)) );
+assert_checkalmostequal([Variance Mean], [expectedV expectedM]);
+// With a scalar
+[Variance, Mean] = variance( x, orien, 1 );
+assert_checkalmostequal([Variance Mean], [expectedV expectedM]);
+
+orien = "c";
+expectedM = [1; 1; 1];
+expectedV = [0.05; 0.81; 0.305];
+
+// Voluntarily passing a bad a priori mean, to check that we get the same results with scalar
+[Variance, Mean] = variance( x, orien, ones(mean(x,orien)) );
+assert_checkalmostequal([Variance Mean], [expectedV expectedM]);
+// With a scalar
+[Variance, Mean] = variance( x, orien, 1 );
+assert_checkalmostequal([Variance Mean], [expectedV expectedM]);
+
+
+//--------------------------------------------------------------
+// Variancef
+
+x = [0.2113249 0.0002211 0.6653811; 0.7560439 0.9546254 0.6283918];
+fre = [1 2 3; 3 4 3];
+orien = "r";
+refM = [1 1 1];
+refV = [0.889522663062 0.593015108708 0.593015108708];
+
+[v, m] = variancef( x, fre, orien, ones(meanf(x,fre,orien)) );
+assert_checkalmostequal([v m], [refV refM]);
+
+[v, m] = variancef( x, fre, orien, 1 );
+assert_checkalmostequal([v m], [refV refM]);
+
+orien = "c";
+refM = [1; 1];
+refV = [0.288922678414; 0.019966608736];
+[v, m] = variancef( x, fre, orien, ones(meanf(x,fre,orien)) );
+assert_checkalmostequal([v m], [refV refM]);
+
+[v, m] = variancef( x, fre, orien, 1 );
+assert_checkalmostequal([v m], [refV refM]);
index 5c80c8f..dd46fde 100644 (file)
@@ -26,31 +26,28 @@ expectedM = 46.25;
 expectedV = 712.5;
 assert_checkequal([Variance Mean], [expectedV expectedM]);
 // With x as a row vector and specified orient
-data = [10 20 30 40; 50 60 70 90];
 [Variance, Mean] = variance(data, 1);
 expectedM = [30 40 50 65];
 expectedV = [800 800 800 1250];
 assert_checkequal([Variance Mean], [expectedV expectedM]);
 // With x as a row vector and specified orient
-data = [10 20 30 40; 50 60 70 90];
 [Variance, Mean] = variance(data, "r");
 expectedM = [30 40 50 65];
 expectedV = [800 800 800 1250];
 assert_checkequal([Variance Mean], [expectedV expectedM]);
 // With x as a row vector and specified orient
-data = [10 20 30 40; 50 60 70 90];
 [Variance, Mean] = variance(data, 2);
 expectedM = [50; 135]/2;
 expectedV = [500; 875]/3;
 assert_checkequal([Variance Mean], [expectedV expectedM]);
 // With x as a row vector and specified orient
-data = [10 20 30 40; 50 60 70 90];
 [Variance, Mean] = variance(data, "c");
 expectedM = [50; 135]/2;
 expectedV = [500; 875]/3;
 assert_checkequal([Variance Mean], [expectedV expectedM]);
 // With x as a complex row vector and 1 argument
-a = [ 0.9, 0.7;
+a = [
+0.9, 0.7;
 0.1, 0.1;
 0.5, 0.4 ];
 data = a + a * 2 * %i;
@@ -59,40 +56,34 @@ expectedM = 0.45+0.9*%i;
 expectedV = 0.515;
 assert_checkequal([Variance Mean], [expectedV expectedM]); // Must be variance(real(data)) + variance(imag(data))
 // With x as a complex row vector and computation by column
-a = [ 0.9, 0.7;
-0.1, 0.1;
-0.5, 0.4 ];
 data = a + a * 2 * %i;
 [Variance, Mean] = variance(data, 1);
 expectedM = [0.5+%i 0.4+0.8*%i];
 expectedV = [0.8 0.45];
 assert_checkalmostequal([Variance Mean], [expectedV expectedM]);
 // With x as a complex row vector and computation by row
-a = [ 0.9, 0.7;
-0.1, 0.1;
-0.5, 0.4 ];
 data = a + a * 2 * %i;
 [Variance, Mean] = variance(data, 2);
 expectedM = [0.8+1.6*%i; 0.1+0.2*%i; 0.45+0.9*%i];
 expectedV = [0.1; 0; 0.025];
 assert_checkalmostequal([Variance Mean], [expectedV expectedM]);
-// Normalization with N-1
-x = [0.9    0.7
+// Normalization with N-1, no a priori mean
+x = [
+0.9    0.7
 0.1    0.1
-0.5    0.4];
+0.5    0.4 ];
 orien = 1;
-w = 0;
-[Variance, Mean] = variance(x,orien,w);
+[Variance, Mean] = variance(x,orien);
 expectedM = [0.5 0.4];
 expectedV = [0.16 0.09];
 assert_checkalmostequal([Variance Mean], [expectedV expectedM]);
-// Normalization with N
-x = [0.9    0.7
-0.1    0.1
-0.5    0.4];
-orien = 1;
-w = 1;
-[Variance, Mean] = variance(x,orien,w);
+// Normalization with N, the a priori can be known
+[Variance, Mean] = variance(x,orien,mean(x,orien));
+expectedM = [0.5 0.4];
+expectedV = [0.32 0.18]/3;
+assert_checkalmostequal([Variance Mean], [expectedV expectedM]);
+// Biased variance
+[Variance, Mean] = variance(x,orien,%nan);
 expectedM = [0.5 0.4];
 expectedV = [0.32 0.18]/3;
 assert_checkalmostequal([Variance Mean], [expectedV expectedM]);
index d7909c6..31594da 100644 (file)
@@ -27,32 +27,29 @@ expectedM = 46.25;
 expectedV = 712.5;
 assert_checkequal([Variance Mean], [expectedV expectedM]);
 // With x as a row vector and specified orient
-data = [10 20 30 40; 50 60 70 90];
 [Variance, Mean] = variance(data, 1);
 expectedM = [30 40 50 65];
 expectedV = [800 800 800 1250];
 assert_checkequal([Variance Mean], [expectedV expectedM]);
 // With x as a row vector and specified orient
-data = [10 20 30 40; 50 60 70 90];
 [Variance, Mean] = variance(data, "r");
 expectedM = [30 40 50 65];
 expectedV = [800 800 800 1250];
 assert_checkequal([Variance Mean], [expectedV expectedM]);
 // With x as a row vector and specified orient
-data = [10 20 30 40; 50 60 70 90];
 [Variance, Mean] = variance(data, 2);
 expectedM = [50; 135]/2;
 expectedV = [500; 875]/3;
 assert_checkequal([Variance Mean], [expectedV expectedM]);
 // With x as a row vector and specified orient
-data = [10 20 30 40; 50 60 70 90];
 [Variance, Mean] = variance(data, "c");
 expectedM = [50; 135]/2;
 expectedV = [500; 875]/3;
 assert_checkequal([Variance Mean], [expectedV expectedM]);
 
 // With x as a complex row vector and 1 argument
-a = [ 0.9, 0.7;
+a = [
+0.9, 0.7;
 0.1, 0.1;
 0.5, 0.4 ];
 data = a + a * 2 * %i;
@@ -62,41 +59,35 @@ expectedV = 0.515;
 assert_checkequal([Variance Mean], [expectedV expectedM]); // Must be variance(real(data)) + variance(imag(data))
 
 // With x as a complex row vector and computation by column
-a = [ 0.9, 0.7;
-0.1, 0.1;
-0.5, 0.4 ];
 data = a + a * 2 * %i;
 [Variance, Mean] = variance(data, 1);
 expectedM = [0.5+%i 0.4+0.8*%i];
 expectedV = [0.8 0.45];
 assert_checkalmostequal([Variance Mean], [expectedV expectedM]);
 // With x as a complex row vector and computation by row
-a = [ 0.9, 0.7;
-0.1, 0.1;
-0.5, 0.4 ];
 data = a + a * 2 * %i;
 [Variance, Mean] = variance(data, 2);
 expectedM = [0.8+1.6*%i; 0.1+0.2*%i; 0.45+0.9*%i];
 expectedV = [0.1; 0; 0.025];
 assert_checkalmostequal([Variance Mean], [expectedV expectedM]);
 
-// Normalization with N-1
-x = [0.9    0.7
+// Normalization with N-1, no a priori mean
+x = [
+0.9    0.7
 0.1    0.1
-0.5    0.4];
+0.5    0.4 ];
 orien = 1;
-w = 0;
-[Variance, Mean] = variance(x,orien,w);
+[Variance, Mean] = variance(x,orien);
 expectedM = [0.5 0.4];
 expectedV = [0.16 0.09];
 assert_checkalmostequal([Variance Mean], [expectedV expectedM]);
-// Normalization with N
-x = [0.9    0.7
-0.1    0.1
-0.5    0.4];
-orien = 1;
-w = 1;
-[Variance, Mean] = variance(x,orien,w);
+// Normalization with N, the a priori can be known
+[Variance, Mean] = variance(x,orien,mean(x,orien));
+expectedM = [0.5 0.4];
+expectedV = [0.32 0.18]/3;
+assert_checkalmostequal([Variance Mean], [expectedV expectedM]);
+// Biased variance
+[Variance, Mean] = variance(x,orien,%nan);
 expectedM = [0.5 0.4];
 expectedV = [0.32 0.18]/3;
 assert_checkalmostequal([Variance Mean], [expectedV expectedM]);
index b954e09..3cb5eaf 100644 (file)
@@ -20,3 +20,21 @@ refM = [0.367985066667; 0.79718087];
 refV = [0.049647428728; 0.006107864498];
 [v, m] = variancef(x, fre, "c");
 assert_checkalmostequal([v m], [refV refM]);
+// With the a priori mean
+refM = 0.63623244375;
+refV = 0.090053830785;
+[v, m] = variancef(x, fre, "*", meanf(x,fre));
+assert_checkalmostequal([v m], [refV refM]);
+refM = [0.61986415 0.636490633333 0.64688645];
+refV = [0.35977704 0.23985136 0.23985136];
+[v, m] = variancef(x, fre, "r", meanf(x,fre,"r"));
+assert_checkalmostequal([v m], [refV refM]);
+refM = [0.367985066667; 0.79718087];
+refV = [0.041372857273; 0.005497078047];
+[v, m] = variancef(x, fre, "c", meanf(x,fre,"c"));
+assert_checkalmostequal([v m], [refV refM]);
+// Biased variance
+refM = 0.63623244375;
+refV = 0.090053830785;
+[v, m] = variancef(x, fre, "*", %nan);
+assert_checkalmostequal([v m], [refV refM]);
index f2d42ee..15d3c91 100644 (file)
@@ -24,3 +24,25 @@ refM = [0.367985066667; 0.79718087];
 refV = [0.049647428728; 0.006107864498];
 [v, m] = variancef(x, fre, "c");
 assert_checkalmostequal([v m], [refV refM]);
+
+// With the a priori mean
+refM = 0.63623244375;
+refV = 0.090053830785;
+[v, m] = variancef(x, fre, "*", meanf(x,fre));
+assert_checkalmostequal([v m], [refV refM]);
+
+refM = [0.61986415 0.636490633333 0.64688645];
+refV = [0.35977704 0.23985136 0.23985136];
+[v, m] = variancef(x, fre, "r", meanf(x,fre,"r"));
+assert_checkalmostequal([v m], [refV refM]);
+
+refM = [0.367985066667; 0.79718087];
+refV = [0.041372857273; 0.005497078047];
+[v, m] = variancef(x, fre, "c", meanf(x,fre,"c"));
+assert_checkalmostequal([v m], [refV refM]);
+
+// Biased variance
+refM = 0.63623244375;
+refV = 0.090053830785;
+[v, m] = variancef(x, fre, "*", %nan);
+assert_checkalmostequal([v m], [refV refM]);