Turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Intel Community
- Software
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library
- MKL Summary Statistics computation of skewness and kurtosis seems to be broken when the standard deviation is more than four orders of magnitude small...

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

tennican

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-04-2011
03:32 PM

163 Views

Can anyone tell me how to fix this, please?

For example, given normal random samples around 2, I should get skewness and excess kurtosis close to zero. However, at a standard deviation of 1.0e-4, Excess kurtosis jumps to -311. At 1.0e-5, skewness jumps to 11. At 1.0e-6, the absolute values of both skewness and kurtosis are gigantic.

In R, I get the correct values much longer.

I have confirmed that I can workaround this problem by copying the data, subtracting the mean and computing the moments with another MKL task.

But, man, this is ugly.

Here is my C++ code for computing summary stats with MKL

------------------------------------------------------------------

// Setup a Summary Statistic task pointing at the data

int status;

VSLSSTaskPtr task;

MKL_INT numRows = m_count;

MKL_INT numCols = 1;

MKL_INT xstorage = VSL_SS_MATRIX_STORAGE_ROWS;

int indices[1] = {1};

status = vsldSSNewTask( &task, &numCols, &numRows, &xstorage, &(m_data[0]), 0, indices );

// Setup 1st, 2nd and 3rd quartile computation

double quants[3];

double q_order[3] = {.25,.5,.75};

MKL_INT q_order_n = 3;

status = vsldSSEditQuantiles( task, &q_order_n, q_order, quants, NULL, NULL );

// Setup max and min computation

double min_est = m_data[0];

double max_est = m_data[0];

status = vsldSSEditTask(task, VSL_SS_ED_MIN, &min_est);

status = vsldSSEditTask(task, VSL_SS_ED_MAX, &max_est);

// Setup moment computation. Computing kurtosis requires all the other moments.

double mean_est, skewness_est, kurtosis_est;

double raw2mom, raw3mom, raw4mom;

double cen2mom, cen3mom, cen4mom;

status = vsldSSEditTask(task, VSL_SS_ED_SKEWNESS, &skewness_est);

status = vsldSSEditTask(task, VSL_SS_ED_KURTOSIS, &kurtosis_est);

status = vsldSSEditMoments( task, &mean_est, &raw2mom, &raw3mom, &raw4mom, &cen2mom, &cen3mom, &cen4mom );

// Request five number summary and moment computation

MKL_UINT64 mask = VSL_SS_MEAN | VSL_SS_KURTOSIS | VSL_SS_SKEWNESS |

VSL_SS_2R_MOM | VSL_SS_3R_MOM | VSL_SS_4R_MOM |

VSL_SS_2C_MOM | VSL_SS_3C_MOM | VSL_SS_4C_MOM |

VSL_SS_QUANTS | VSL_SS_MIN | VSL_SS_MAX;

status = vsldSSCompute(task, mask, VSL_SS_METHOD_FAST);

--------------------------------------------------------------------------

Here is the R code that I used to generate the test data as well and check skewness and kurtosis:

-----------------------------------------------------------------------

library(moments)

setwd("C:/Users/stennican/Documents/Rcode")

N = 500

normDist = rnorm(N)

trial = seq(1,N)

twos = rep(2,N)

twoSd1 = twos + 1e-1*normDist

twoSd2 = twos + 1e-2*normDist

twoSd3 = twos + 1e-3*normDist

twoSd4 = twos + 1e-4*normDist

twoSd5 = twos + 1e-5*normDist

twoSd7 = twos + 1e-7*normDist

twoSd10 = twos + 1e-10*normDist

twoSd13 = twos + 1e-13*normDist

zeros = rep(0,N)

zeroSd1 = 1e-1*normDist

zeroSd2 = 1e-2*normDist

zeroSd3 = 1e-3*normDist

zeroSd4 = 1e-4*normDist

zeroSd5 = 1e-5*normDist

zeroSd7 = 1e-7*normDist

zeroSd10 = 1e-10*normDist

zeroSd13 = 1e-13*normDist

out = data.frame(trial,zeros,twos,

zeroSd1,zeroSd2,zeroSd3,zeroSd4,zeroSd5,

zeroSd7,zeroSd10,zeroSd13,

twoSd1,twoSd2,twoSd3,twoSd4,twoSd5,

twoSd7,twoSd10,twoSd13)

apply(out,2,sd)

apply(out,2,skewness)

apply(out,2,kurtosis)

options(digits = 22)

write.csv(out,"testB42382.csv",row.names=FALSE)

Link Copied

13 Replies

Andrey_N_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-07-2011
01:13 AM

163 Views

Thanks for the post and the detailed description of the issue. We will have a look atthe dataset and the test case on our side. Can you please clarifywhat the value of the m_count variable in your code is? Is it equal to 500?

Thanks,

Andrey

tennican

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-07-2011
08:58 AM

163 Views

Yes, m_count is 500 because I just imported the data generated in the R code where N = 500.

I believe there is some change in the behavior if you vary sample size from 50 to 500 to 5000.

But, I wasn't testing in an organized fashion at the time. So, I'd have to retest to be sure.

thanks for you help, Scott

Sergey_M_Intel2

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-07-2011
09:02 PM

163 Views

May I clarify one more thing? Is this behavior reproducible on different random sequences? In other words, isn't it because of large statistical error due to relatively small sample size?

Sergey

tennican

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-08-2011
09:15 AM

163 Views

Yes, it is reproducible with various random samples of various sizes.

Most would call 500 a large sample since the CLT is working for you big time at that point.

I haven't tried any small samples of less than 30 yet.

Let me know if you need me to do additional testing and send you the results.

thanks, Scott

Andrey_N_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-08-2011
09:33 AM

163 Views

We did the analysis of the datasetsyou work with,computed confidence intervals for sample kurtosis for the dataset from Gaussian distribution, and the output of Summary Stats routines for kurtosis/skewness estimation.

The current conlusion is that it makes sensefor us tounderstand the opportunities for having the algorithms for computation of the moments/covariance which would bemore tolerant to numerical errors to address the case of the datasets with variance several orders smaller than its expectation.

I wonder if the work-around you mention earlier and which sounds like a two pass algorithm would currently work for your needs?

Also, canyouplease specify howoftenyouget the datasets with such structure in your work?

Iwould appreciate if you could also clarify if you work with in-memory datasets only, or data that can have such structure and which does not fit into memory of the computer is also the case?

Thanks,

Andrey

tennican

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-08-2011
09:59 AM

163 Views

Our tool makes it very easy for people with no statistical background to access any sort of database, transform columns of data from that database and then request summary stats for the column. So, while it doesn't make sense to care about the kurtosis of a variable with almost no variance, we will have to compute this at some point. For a constant, kurtosis is undefined and it is correct to say so. But, for small variance, kurtosis and skewness are well defined and we should report them correctly.

The data on which MKL operates is in an in-memory database.

As I said before, I can work around the problem by first computing the mean of the original data, subtracting the mean from a copy of the data and computing the higher moments using the demeaned data.

The extra cost in time and memory are acceptable in the short term.

How long do you think it will take for MKL to be improved to match the precision of R in this respect?

thanks, Scott

Andrey_N_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-08-2011
11:21 PM

163 Views

We would do the deeper analysis of your requestwithin the next several weeks. Once the analysisis completedwe would understand the next steps.

Does it sound reasonable for you?

Thanks,

Andrey

tennican

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-09-2011
11:12 AM

163 Views

Sure, I'm not in a huge rush.

The workaround only triggers in this very special case.

So, it won't impact performance.

thanks, Scott

Sergey_M_Intel2

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-09-2011
09:27 PM

163 Views

One more clarification to be 100% sure that we address your needs.

You indicated that majority of your customers are not experts in statistics so you would like to provide basic statistics capabilities as simple (to use) as possible. Does it mean that you want to avoid giving customers some controls to let your software know about properties of incoming dataset?

For example, Andrey gave you an example of the one-pass method for the mean estimate, which is rather expensive (compared to naive implementation) but it is more numerically robust. At the same time for datasets with certain known a priori properties the naive method may work quite well either (e.g. relatively small size dataset, mean isn't too far from 0, etc).

Would you like us to consider such cases either?

The opposite direction is an option either. One-pass method described by Andrey has its own limitations. There are more robust but at the same time more computationally intensive methods for moment estimates. Should we be looking into that either?

Many thanks!

Sergey

tennican

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-10-2011
10:50 AM

163 Views

Actually, yes, if you could.

If MKL could switch from naive to one-pass based on parameters of the data, it would be awesome.

In fact, the current vsldSSNewTask function requires me to specify number of rows.

So, if you switched to the naive method whenever, the number of rows is not large that would be great.

I could code that myself. But, it would be cleaner if it were all build in to MKL.

Also, I don't need this yet.

But, for future reference, how would I set up the computation if I don't yet know how many rows I will be feeding to the vsldSS task?

thanks, scott

Sergey_M_Intel2

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-10-2011
09:28 PM

163 Views

Hi Scott,

Perfect. That certainly helps us to define the scope of experiments.

One of Summary Statistics features is the possibility to feed data in blocks, so-called streaming data processing. That may help. Andrey may have some comments on that either.But, for future reference, how would I set up the computation if I don't yet know how many rows I will be feeding to the vsldSS task?

Regards,

Sergey

Andrey_N_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-11-2011
06:10 AM

163 Views

You might want to have a look atSummary Statistics Application Notes, http://software.intel.com/sites/products/documentation/hpc/mkl/sslnotes/index.htm, chapter "Processing Data in Blocks" which describes the usage model for streaming data processing. In Intel MKL Manual, Chapter StatisticalFunctions/Summary Statistics,Section Task Computation Routine you would find the tablelisting estimates that support streaming processing. It's worth noting thatbesidesbasic estimates like variance-covariance or meanthis feature is available for quantiles.

Thanks,

Andrey

tennican

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-11-2011
03:04 PM

163 Views

thanks, Scott

Topic Options

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page

For more complete information about compiler optimizations, see our Optimization Notice.