vslsrobustcovariancerow.c

使用INTEL矢量统计类库的程序,包括以下功能: &#61623 Raw and central moments up to 4th order &#61623 Kurtosis and

/*******************************************************************************
!                             INTEL CONFIDENTIAL
!  Copyright(C) 2007-2008 Intel Corporation. All Rights Reserved.
!  The source code contained  or  described herein and all documents related to
!  the source code ("Material") are owned by Intel Corporation or its suppliers
!  or licensors.  Title to the  Material remains with  Intel Corporation or its
!  suppliers and licensors. The Material contains trade secrets and proprietary
!  and  confidential  information of  Intel or its suppliers and licensors. The
!  Material  is  protected  by  worldwide  copyright  and trade secret laws and
!  treaty  provisions. No part of the Material may be used, copied, reproduced,
!  modified, published, uploaded, posted, transmitted, distributed or disclosed
!  in any way without Intel's prior express written permission.
!  No license  under any  patent, copyright, trade secret or other intellectual
!  property right is granted to or conferred upon you by disclosure or delivery
!  of the Materials,  either expressly, by implication, inducement, estoppel or
!  otherwise.  Any  license  under  such  intellectual property  rights must be
!  express and approved by Intel in writing.
!
!*******************************************************************************
!  Content:
!    Computation of robust covariance matrix and mean Example Program Text
!******************************************************************************/


#include <stdio.h>
#include "vsl_ss.h"
#include "mkl.h"


/* parameters of the task */
#define DIM                   5    /* dimension of task */
#define N                  10000   /* number of observations */

static float x[DIM*N];            /* matrix of observations */
static float t[DIM];              /* vector of means */
static float cov[DIM*DIM];        /* covariance matrix */ 

static float t_est[DIM];       /* vector of robust mean estimates */
static float cov_est[DIM*DIM]; /* matrix of robust covariances */

static float W[2];             /* array of weights used in covariance 
                                   estimation */

#define RETURN_ON_ERROR                 \
    if(errcode<0)                        \
    {                                   \
        printf("Error: %i\n", errcode);  \
        printf("\nTEST FAILED\n");		\
        return errcode;                  \
    }

int GenerateDataset( MKL_INT p, MKL_INT n, float t[], float c[], 
                     MKL_INT eps, float m, float sigma, float r[] );

int main()
{
    int errcode;
    MKL_INT i, j;
    MKL_INT p, n;
    MKL_INT robust_cov_storage, robust_params_n, eps, storage_format_x;
    MKL_INT cov_storage;
    float m, sigma;
    float breakdown_point, arp, method_accuracy, iter_num;
    float robust_method_params[VSL_SS_TBS_PARAMS_N];
    VSLSSTaskPtr task=0;

    p = DIM;
    n = N;
    storage_format_x   = VSL_SS_MATRIX_ROWS_STORAGE;
    cov_storage        = VSL_SS_MATRIX_FULL_STORAGE;
    robust_cov_storage = VSL_SS_MATRIX_FULL_STORAGE;
    robust_params_n    = VSL_SS_TBS_PARAMS_N;

    /* Generate covariance matrix and vector of means */
    for ( i = 0; i < p; i++ )
    {
        t[i] = 0.0;
        for ( j = 0;   j < i; j++ )   cov[i*p+j] = 1.0 / (2.0 * p);
        for ( j = i+1; j < p;   j++ ) cov[i*p+j] = 1.0 / (2.0 * p);
        cov[i*p+i] = 1.0;
    }

    printf("Dimension of the task: %d\n", p);
    printf("Number of observations: %d\n\n",n);

    printf("Original covariance matrix:\n");
    for ( i = 0; i < p; i++ )
    {
         for( j = 0; j < p; j++ ) printf("%f ", cov[i*DIM+j]);
         printf("\n");
    }
    printf("\n");
    printf("Original vector of means:\n");
    for ( i = 0; i < p; i++ )   printf("%f, ", t[i]);
    printf("\n\n");

    eps   = 20;    /* ratio of outliers in the dataset */
    m     = 5.0;    /* mean of the outliers */
    sigma = 0.1;    /* coefficient to compute covarince of outliers */
    errcode = GenerateDataset( p, n, t, cov, eps, m, sigma, x );
    if ( errcode < 0 )
    {
        printf("Dataset was not generated\n");
        return 0;
    }

    /* Compute "classical" covariance and mean estimates */
    /* Initialize arrays which will hold estimates */
    for ( i = 0; i < p*p; i++ ) cov_est[i]=0.0;
    for ( i = 0; i < p; i++ ) t_est[i]=0.0;
    W[0]=0.0;
    W[1]=0.0;

    /* Create task */
    errcode = vslsSSNewTask( &task, &p, &n, &storage_format_x, x, 0, 0 );
	RETURN_ON_ERROR;

    /* Register array of weights in the task */
    errcode = vslsSSEditTask( task, VSL_SS_ACCUMULATED_WEIGHT, W );
	RETURN_ON_ERROR;

    /* Register array for covariance estimate */  
    errcode = vslsSSEditCovCor( task, t_est, cov_est, &cov_storage, 0, 0 );
	RETURN_ON_ERROR;

    errcode = vslsSSCompute( task, 
                         VSL_SS_COVARIANCE_MATRIX, VSL_SS_FAST_METHOD );
	RETURN_ON_ERROR;

    printf("Classical covariance estimate:\n");
	RETURN_ON_ERROR;

    printf("\n");
    printf("Classical mean estimate:\n");
    for ( i = 0; i < p; i++ )   printf("%f, ", t_est[i]);
    printf("\n\n");


    for ( i = 0; i < p*p; i++ ) cov_est[i]=0.0f;
    for ( i = 0; i < p; i++ ) t_est[i]=0.0f;


    /* Compute robust  covariance and mean estimates */
    breakdown_point = 0.4;
    arp             = 0.001;
    method_accuracy = 0.001;
    iter_num        = 20;

    robust_method_params[0] = breakdown_point;
    robust_method_params[1] = arp;
    robust_method_params[2] = method_accuracy;
    robust_method_params[3] = iter_num;

    errcode =  vslsSSEditRobustCovariance( task, &robust_cov_storage, 
                              &robust_params_n, robust_method_params,
                                                     t_est, cov_est );
	RETURN_ON_ERROR;

    errcode = vslsSSCompute( task, VSL_SS_ROBUST_COVARIANCE, 
                                   VSL_SS_ROBUST_TBS_METHOD );
	RETURN_ON_ERROR;

    printf("Robust covariance estimate:\n");
    for ( i = 0; i < p; i++ )
    {
         for( j = 0; j < p; j++ ) printf("%f ", cov_est[i*DIM+j]);
         printf("\n");
    }
    printf("\n");
    printf("Robust mean estimate:\n");
    for ( i = 0; i < p; i++ )   printf("%f, ", t_est[i]);
    printf("\n\n");

    errcode = vslSSDeleteTask( &task );
	RETURN_ON_ERROR;

    printf("\nTEST PASSED\n");

    return 0;
}



static float r1[DIM*N], r2[DIM*N], mean[DIM], T[DIM*DIM];

int GenerateDataset( MKL_INT p, MKL_INT n, float t[], float c[], 
                     MKL_INT eps, float m, float sigma, float r[] )
{
    int errcode;
    VSLStreamStatePtr stream;
    int i,j, mc;

    char uplo;
    int lda;
    int info;


    errcode = 0;
    mc = (n * eps)/100; /* number of outliers */


    for ( i = 0; i < p; i++ ) 
    {
        for ( j = 0; j < p; j++ ) T[i*p+j] = c[i*p+j];
    }

    uplo = 'U';                
    lda=p;
    spotrf( &uplo, &p, T, &lda, &info );

    errcode = vslNewStream( &stream, VSL_BRNG_MT19937, 7777777 );
    if ( errcode < 0 ) return errcode;

    /* Generate "good" points of the matrix of observations */
    if ( n-mc )
    {
        errcode = vsRngGaussianMV( VSL_METHOD_DGAUSSIAN_ICDF, stream, n-mc,
                                   r1, p, VSL_MATRIX_STORAGE_FULL, t, T );
        if ( errcode < 0 ) return errcode;
    }

    /* Generate "bad" points of the matrix of observations */
    if ( mc )
    {
        for ( i = 0; i < p; i++ ) 
        {
            mean[i] = m;
            for ( j = 0; j < p; j++ ) 
            {
                T[i*p+j] = eps * c[i*p+j];
            }
        }
        
        spotrf( &uplo, &p, T, &lda, &info );
        errcode = vsRngGaussianMV( VSL_METHOD_DGAUSSIAN_ICDF, stream, mc, 
                                   r2, p, VSL_MATRIX_STORAGE_FULL, mean, T );
        if ( errcode < 0 ) return errcode;
    }

    /* Copy results to the matrix r */
    for ( i = 0; i < n-mc; i++ )
    {
       for ( j = 0; j < p; j++ ) r[j*n+i] = r1[i*p+j];
    }


    for ( i = n-mc; i < n; i++ )
    {
       for ( j = 0; j < p; j++ ) r[j*n+i] = r2[(i-(n-mc))*p+j];
    }

    errcode = vslDeleteStream( &stream );

    return errcode;
}