Different result between Octave and GSL (Gnu Scientific Library) C function call

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

Different result between Octave and GSL (Gnu Scientific Library) C function call

Tallis Huther da Costa
I wanted to call the function "gsl_stats_wvariance" within Octave.

The documentation for the function is here (https://www.gnu.org/software/gsl/doc/html/statistics.html).

I made a C++ code to obtain the result of the function and I made a corresponding Octave script to confirm the result of the function, but the results didn't match.

At first, the result from the function gsl_stats_wvariance was different than my Octave code:

I wrote a file called "gsl_stats_wvariance_test_c.cc".The contents of the file is this:
// --------------------------------------------------------------------------------
#include <octave/oct.h>
#include <octave/parse.h>
#include <gsl/gsl_statistics.h>
// in shell:
// sudo apt-get install libgsl-dev
// in Octave:
// filename = "gsl_stats_wvariance_test_c.cc"
// mkoctfile ("--verbose", "-lgsl", filename);

DEFUN_DLD (gsl_stats_wvariance_test_c, args, nargout,
           "Hello World Help String")
{
  octave_stdout << "Hello World has "
                << args.length () << " input arguments and "
                << nargout << " output arguments.\n";
int nargin = args.length ();
double data[] = {0, 2, 2, 3, 1, 5, 2.5, 7, 3, 3.9, 1.1};
size_t stride = 1;
double w[] = {2, 20, 1.8, 3, 4, 5.4, 4.5, 1.2, 30, 1.9, 10.1};
size_t wstride = 1;
size_t n = sizeof (data)/sizeof (double);
octave_stdout << "n: " << n << "\n";

double result = gsl_stats_wvariance(w, wstride, data, stride, n);
return octave_value (result);
}
// --------------------------------------------------------------------------------

I compiled this code by typing in the Octave prompt:
>> filename = "gsl_stats_wvariance_test_c.cc"
filename = gsl_stats_wvariance_test_c.cc

>> mkoctfile ("--verbose", "-lgsl", filename);
g++ -c -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC -I/usr/include/octave-5.2.0/octave/.. -I/usr/include/octave-5.2.0/octave  -pthread -fopenmp -g -O2 -fdebug-prefix-map=/build/octave-rvRilm/octave-5.2.0=. -fstack-protector-strong -Wformat -Werror=format-security    gsl_stats_wvariance_test_c.cc -o /tmp/oct-EDOZzk.o
g++ -I/usr/include/octave-5.2.0/octave/.. -I/usr/include/octave-5.2.0/octave  -pthread -fopenmp -g -O2 -fdebug-prefix-map=/build/octave-rvRilm/octave-5.2.0=. -fstack-protector-strong -Wformat -Werror=format-security-shared -Wl,-Bsymbolic -Wl,-Bsymbolic-functions -Wl,-z,relro  -o gsl_stats_wvariance_test_c.oct  /tmp/oct-EDOZzk.o   -lgsl  -L/usr/lib/x86_64-linux-gnu  -Wl,-Bsymbolic-functions -Wl,-z,relro
>> gsl_stats_wvariance_test_c
Hello World has 0 input arguments and 0 output arguments.
n: 11
ans =  1.7643

As you can see, the answer is 1.7643.

Now, from the function's documentation:
image.png

So, I wrote this little Octave script to confirm that the call was correct:
data = [0 2 2 3 1 5 2.5 7 3 3.9 1.1];
_mean = mean (data);
stride = 1;
w = [2 20 1.8 3 4 5.4 4.5 1.2 30 1.9 10.1];
wstride = 1;
n = numel (data);
expected = (sum (w) / ((sum (w))^2 - sum (w .^ 2))) * sum (w .* ((data - _mean) .^ 2));

The result was:
expected =  1.8427

As you can see, the result from the C function call is 1.7643 while the result from Octave's calculation is 1.8427.

This difference seems too much compared to previous function calls of the same statistics functions category (data not shown).

I would like some help to figure out where the wrong part is.

The code was compiled in a Ubuntu virtual machine in a Windows host. Octave was running in the virtual machine.
If the same Octave script is run in Octave in Windows, it gives the result:
expected = 1.8427

Here is the Ubuntu version:
$ uname -a
Linux ubuntuuservm 5.4.0-58-generic #64-Ubuntu SMP Wed Dec 9 08:16:25 UTC 2020 x86_64 x86_64 x86_64 GNU/Linux

Here is the gsl version:
$ dpkg-query -l libgsl-dev
Desired=Unknown/Install/Remove/Purge/Hold
| Status=Not/Inst/Conf-files/Unpacked/halF-conf/Half-inst/trig-aWait/Trig-pend
|/ Err?=(none)/Reinst-required (Status,Err: uppercase=bad)
||/ Name           Version          Architecture Description
+++-==============-================-============-===================================================
ii  libgsl-dev     2.5+dfsg-6build1 amd64        GNU Scientific Library (GSL) -- development package

Here is the Octave version:
$ octave --version
GNU Octave, version 5.2.0

Thank you,
Tallis Huther da Costa







Reply | Threaded
Open this post in threaded view
|

Re: Different result between Octave and GSL (Gnu Scientific Library) C function call

mmuetzel
Am 01. Januar 2021 um 07:24 Uhr schrieb "Tallis Huther da Costa":

> >> gsl_stats_wvariance_test_c
> Hello World has 0 input arguments and 0 output arguments.
> n: 11
> ans =  1.7643

> As you can see, the answer is 1.7643.

> Now, from the function's documentation:

> So, I wrote this little Octave script to confirm that the call was correct:
> data = [0 2 2 3 1 5 2.5 7 3 3.9 1.1];
> _mean = mean (data);
> stride = 1;
> w = [2 20 1.8 3 4 5.4 4.5 1.2 30 1.9 10.1];
> wstride = 1;
> n = numel (data);
> expected = (sum (w) / ((sum (w))^2 - sum (w .^ 2))) * sum (w .* ((data - _mean) .^ 2));
>
> The result was:
> expected =  1.8427
 
Looking at the definition of the weighted variance as used by gsl, the *weighted* mean has to be used in its calculation.
Correcting for that gives the same result in Octave as the one you observed when using the gsl function:
data = [0 2 2 3 1 5 2.5 7 3 3.9 1.1];
w = [2 20 1.8 3 4 5.4 4.5 1.2 30 1.9 10.1];
weighted_mean = sum (w .* data) / sum (w);
expected = ((sum (w) / ((sum (w))^2 - sum (w .^ 2))) * sum (w .* ((data - weighted_mean) .^ 2));

>> expected
expected =  1.7643

HTH,
Markus