ks.c

r软件 另一款可以计算核估计的软件包 需安装r软件

  free(y);
  free(dens);

}

/*************************************************************************
* Double sum of normal density values * x * x^T - bivariate
*
* Parameters
* x1 - x values
* x2 - y values
* visigma - vec of inverse Sigma
* detsigma - det of Sigma
* n - number of values
* sum - contains sum of density values
*************************************************************************/

void dmvnorm_2d_xxt_sum(double *x1, double *x2, double *visigma, 
     double *detsigma, int *n, double *sum)
{
    int i, j, k, N;
    double mu[] = {0.0, 0.0};
    double *Xmat0, *Xmat1, *Xmat2, *Xmat3;
    double *x, *y, *dens;

    N = n[0];
    x = malloc(sizeof(double)*N);
    y = malloc(sizeof(double)*N);
    dens = malloc(sizeof(double)*N);
    Xmat0 = malloc(sizeof(double)*N);
    Xmat1 = malloc(sizeof(double)*N);
    Xmat2 = malloc(sizeof(double)*N);
    Xmat3 = malloc(sizeof(double)*N);
     
    for (i = 1; i <= N; i++)
    {
        for (j = i; j <= N; j++)
        {
            x[j - 1] = x1[i - 1] - x1[j - 1];
            y[j - 1] = x2[i - 1] - x2[j - 1];
            Xmat0[j - 1] = x[j - 1] * x[j - 1];
            Xmat2[j - 1] = Xmat1[j - 1] = x[j - 1] * y[j - 1];
            Xmat3[j - 1] = y[j - 1] * y[j - 1];
        }
        dmvnorm_2d(x, y, mu, visigma, detsigma, n, dens);
        for (k = i; k <= N; k++)
        {
            sum[0] += dens[k - 1]*Xmat0[k - 1];
            sum[1] += dens[k - 1]*Xmat1[k - 1];
            sum[2] += dens[k - 1]*Xmat2[k - 1];
            sum[3] += dens[k - 1]*Xmat3[k - 1];
        }
    }

    free(x);
    free(y);
    free(dens);
    free(Xmat0);   
    free(Xmat1);
    free(Xmat2);    
    free(Xmat3);
}

/*************************************************************************
* Double sum used in density derivative estimation - 1st order
*
* Parameters
* x1 - x values
* x2 - y values
* v - vec Sigma
* r - (r1, r2) partial derivative
* n - number of values
* sum - contains sum of density derivative values
*************************************************************************/

void dmvnormd1_2d_sum(double *x1, double *x2, double *vsigma, int *r, int *n,
    double *sum)
{
    int i, j, k, N;
    double sum1;
    double *x, *y, *derivt;

    N = n[0];
    x = malloc(sizeof(double)*N);
    y = malloc(sizeof(double)*N);
    derivt = malloc(sizeof(double)*N);
    sum1 = 0.0;

    for (i = 1; i <= N; i++)
    {
        for (j = i; j <= N; j++)
        {
            x[j - 1] = x1[i - 1] - x1[j - 1];
            y[j - 1] = x2[i - 1] - x2[j - 1];
        }
        dmvnormd1_2d(x, y, vsigma, r, n, derivt);

        for (k = i; k <= N; k++)
            sum1 = sum1 + derivt[k - 1];
    }

    sum[0] = sum1;
    free(x);
    free(y);
    free(derivt);
}


/*************************************************************************
* Double sum used in density derivative estimation - 2nd order
*
* Parameters
* x1 - x values
* x2 - y values
* v - vec Sigma
* r - (r1, r2) partial derivative
* n - number of values
* sum - contains sum of density derivative values
*************************************************************************/

void dmvnormd2_2d_sum(double *x1, double *x2, double *vsigma, int *r, int *n,
    double *sum)
{
    int i, j, k, N;
    double sum1;
    double *x, *y, *derivt;

    N = n[0];
    x = malloc(sizeof(double)*N);
    y = malloc(sizeof(double)*N);
    derivt = malloc(sizeof(double)*N);
    sum1 = 0.0;

    for (i = 1; i <= N; i++)
    {
        for (j = i; j <= N; j++)
        {
            x[j - 1] = x1[i - 1] - x1[j - 1];
            y[j - 1] = x2[i - 1] - x2[j - 1];
        }
        dmvnormd2_2d(x, y, vsigma, r, n, derivt);

        for (k = i; k <= N; k++)
            sum1 = sum1 + derivt[k - 1];
    }

    sum[0] = sum1;
    free(x);
    free(y);
    free(derivt);
}



/*************************************************************************
* Double sum used in density derivative estimation - 3rd order
*
* Parameters
* x1 - x values
* x2 - y values
* v - vec Sigma
* r - (r1, r2) partial derivative
* n - number of values
* sum - contains sum of density derivative values
*************************************************************************/

void dmvnormd3_2d_sum(double *x1, double *x2, double *vsigma, int *r, int *n,
    double *sum)
{
    int i, j, k, N;
    double sum1;
    double *x, *y, *derivt;

    N = n[0];
    x = malloc(sizeof(double)*N);
    y = malloc(sizeof(double)*N);
    derivt = malloc(sizeof(double)*N);
    sum1 = 0.0;

    for (i = 1; i <= N; i++)
    {
        for (j = i; j <= N; j++)
        {
            x[j - 1] = x1[i - 1] - x1[j - 1];
            y[j - 1] = x2[i - 1] - x2[j - 1];
        }
        dmvnormd3_2d(x, y, vsigma, r, n, derivt);  

        for (k = i; k <= N; k++)
            sum1 = sum1 + derivt[k - 1];
    }

    sum[0] = sum1;

    free(x);
    free(y);
    free(derivt);
}


/*************************************************************************
* Double sum used in density derivative estimation - 4th order
*
* Parameters
* x1 - x values
* x2 - y values
* v - vec Sigma
* r - (r1, r2) partial derivative
* n - number of values
* sum - contains sum of density derivative values
*************************************************************************/

void dmvnormd4_2d_sum(double *x1, double *x2, double *vsigma, int *r, int *n,
    double *sum)
{
    int i, j, k, N;
    double sum1;
    double *x, *y, *derivt;

    N = n[0];
    x = malloc(sizeof(double)*N);
    y = malloc(sizeof(double)*N);
    derivt = malloc(sizeof(double)*N);
    sum1 = 0.0;

    for (i = 1; i <= N; i++)
    {
        for (j = i; j <= N; j++)
        {
            x[j - 1] = x1[i - 1] - x1[j - 1];
            y[j - 1] = x2[i - 1] - x2[j - 1];
        }
        dmvnormd4_2d(x, y, vsigma, r, n, derivt);

        for (k = i; k <= N; k++)
            sum1 = sum1 + derivt[k - 1];
    }

    sum[0] = sum1;
    free(x);
    free(y);
    free(derivt);
}

/*************************************************************************
* Double sum used in density derivative estimation - 5th order
*
* Parameters
* x1 - x values
* x2 - y values
* v - vec Sigma
* r - (r1, r2) partial derivative
* n - number of values
* sum - contains sum of density derivative values
*************************************************************************/

void dmvnormd5_2d_sum(double *x1, double *x2, double *vsigma, int *r, int *n,
    double *sum)
{
    int i, j, k, N;
    double sum1;
    double *x, *y, *derivt;

    N = n[0];
    x = malloc(sizeof(double)*N);
    y = malloc(sizeof(double)*N);
    derivt = malloc(sizeof(double)*N);
    sum1 = 0.0;

    for (i = 1; i <= N; i++)
    {
        for (j = i; j <= N; j++)
        {
            x[j - 1] = x1[i - 1] - x1[j - 1];
            y[j - 1] = x2[i - 1] - x2[j - 1];
        }
        dmvnormd5_2d(x, y, vsigma, r, n, derivt);

        for (k = i; k <= N; k++)
            sum1 = sum1 + derivt[k - 1];
    }

    sum[0] = sum1;
    free(x);
    free(y);
    free(derivt);
}

/*************************************************************************
* Double sum used in density derivative estimation - 6th order
*
* Parameters
* x1 - x values
* x2 - y values
* v - vec Sigma
* r - (r1, r2) partial derivative
* n - number of values
* sum - contains sum of density derivative values
*************************************************************************/

void dmvnormd6_2d_sum(double *x1, double *x2, double *vsigma, int *r, int *n,
    double *sum)
{
    int i, j, k, N;
    double sum1;
    double *x, *y, *derivt;

    N = n[0];
    x = malloc(sizeof(double)*N);
    y = malloc(sizeof(double)*N);
    derivt = malloc(sizeof(double)*N);
    sum1 = 0.0;

    for (i = 1; i <= N; i++)
    {
        for (j = i; j <= N; j++)
        {
            x[j - 1] = x1[i - 1] - x1[j - 1];
            y[j - 1] = x2[i - 1] - x2[j - 1];
        }
        dmvnormd6_2d(x, y, vsigma, r, n, derivt);

        for (k = i; k <= N; k++)
            sum1 = sum1 + derivt[k - 1];
    }

    sum[0] = sum1;
    free(x);
    free(y);
    free(derivt);
}



/*************************************************************************
* Double sum of 4th order derivative normal density values * x x^T
*
* Parameters
* x1 - x values
* x2 - y values
* visigma - vec of inverse Sigma
* detsigma - det of Sigma
* n - number of values
* sum - contains sum of density values
*************************************************************************/

void dmvnormd4_2d_xxt_sum(double *x1, double *x2, double *visigma, 
			  int*r, int *n, double *sum)
{
    int i, j, k, N;
    double *Xmat0, *Xmat1, *Xmat2, *Xmat3;
    double *x, *y, *derivt;

    N = n[0];
    x = malloc(sizeof(double)*N);
    y = malloc(sizeof(double)*N);
    derivt = malloc(sizeof(double)*N);
    Xmat0 = malloc(sizeof(double)*N);
    Xmat1 = malloc(sizeof(double)*N);
    Xmat2 = malloc(sizeof(double)*N);
    Xmat3 = malloc(sizeof(double)*N);
     
    for (i = 1; i <= N; i++)
    {
        for (j = i; j <= N; j++)
        {
            x[j - 1] = x1[i - 1] - x1[j - 1];
            y[j - 1] = x2[i - 1] - x2[j - 1];
            Xmat0[j - 1] = x[j - 1] * x[j - 1];
            Xmat2[j - 1] = Xmat1[j - 1] = x[j - 1] * y[j - 1];
            Xmat3[j - 1] = y[j - 1] * y[j - 1];
        }
        dmvnormd4_2d(x, y, visigma, r, n, derivt);
        for (k = i; k <= N; k++)
        {
            sum[0] += derivt[k - 1]*Xmat0[k - 1];
            sum[1] += derivt[k - 1]*Xmat1[k - 1];
            sum[2] += derivt[k - 1]*Xmat2[k - 1];
            sum[3] += derivt[k - 1]*Xmat3[k - 1];
        }
    }

    free(x);
    free(y);
    free(derivt);
    free(Xmat0);   
    free(Xmat1);
    free(Xmat2);    
    free(Xmat3);
}

/*************************************************************************
 * Bivariate t density
 *
 * Parameters
 * x - x values
 * y - y values
 * mu - mean
 * viSigma - vector form of inverse of dispersion matrix
 * 
 * n - number of values 
 * dens - contains density values
 *************************************************************************/
void dmvt_2d(double *x, double *y, double *mu, double *viSigma,  
	     double *df, int *n, double *dens)
{
  double xmu[2], d;
  int j;