ks.c

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

  
  d = 2.0;
  for (j = 1; j <= n[0]; j++)
  {
    xmu[0] = x[j - 1] - mu[0];
    xmu[1] = y[j - 1] - mu[1];
    dens[j - 1] = pow( 1 + mult(xmu, viSigma, xmu, d)/df[0], -0.5*(-d +df[0]));
  } 
}

/*************************************************************************
* 3-dim normal density 
*************************************************************************/

void dmvnorm_3d(double *x1, double *x2, double *x3, double *mu, 
		double *visigma, double *detsigma, int *n, double *dens)
{
  double norm;
  double xmu[3];
  int j, d;
  
  d = 3;
  norm = 1/sqrt(pow(2*M_PI, d) * detsigma[0]);
  for (j = 1; j <= n[0]; j++)
  {
    xmu[0] = x1[j - 1] - mu[0];
    xmu[1] = x2[j - 1] - mu[1];
    xmu[2] = x3[j - 1] - mu[2]; 
    dens[j - 1] = norm * exp(-0.5 * mult(xmu, visigma, xmu, d));
  }
}

/*************************************************************************
* 3-dim normal density - double sum 
*************************************************************************/

void dmvnorm_3d_sum(double *x1, double *x2, double *x3, double *visigma, 
		    double *detsigma, int *n, double *sum)
{
  int i, j, k, N;
  double mu[] = {0.0, 0.0, 0.0};
  double sum1;
  double *y1, *y2, *y3, *dens;
  
  N = n[0];
  y1 = malloc(sizeof(double)*N);
  y2 = malloc(sizeof(double)*N);
  y3 = malloc(sizeof(double)*N);
  dens = malloc(sizeof(double)*N);
  sum1 = 0.0;
 
  for (i = 1; i <= N; i++)
  {
    for (j = i; j <= N; j++)
    {
      y1[j - 1] = x1[i - 1] - x1[j - 1];
      y2[j - 1] = x2[i - 1] - x2[j - 1];
      y3[j - 1] = x3[i - 1] - x3[j - 1];
    }
    dmvnorm_3d(y1, y2, y3, mu, visigma, detsigma, n, dens);
    for (k = i; k <= N; k++)
      sum1 = sum1 + dens[k - 1];
  }

  sum[0] = sum1;
  free(y1);
  free(y2);
  free(y3);
  free(dens);
}



/*************************************************************************
* 4-dim normal density 
*************************************************************************/

void dmvnorm_4d(double *x1, double *x2, double *x3, double *x4, double *mu, 
		double *visigma, double *detsigma, int *n, double *dens)
{
  double norm;
  double xmu[4];
  int j, d;
  
  d = 4;
  norm = 1/sqrt(pow(2*M_PI, d) * detsigma[0]);
  for (j = 1; j <= n[0]; j++)
  {
    xmu[0] = x1[j - 1] - mu[0];
    xmu[1] = x2[j - 1] - mu[1];
    xmu[2] = x3[j - 1] - mu[2]; 
    xmu[3] = x4[j - 1] - mu[3];
    dens[j - 1] = norm * exp(-0.5 * mult(xmu, visigma, xmu, d));
  }
}

/*************************************************************************
* 4-dim normal density - double sum 
*************************************************************************/

void dmvnorm_4d_sum(double *x1, double *x2, double *x3, double *x4, double *visigma, 
		    double *detsigma, int *n, double *sum)
{
  int i, j, k, N;
  double mu[] = {0.0, 0.0, 0.0, 0.0};
  double sum1;
  double *y1, *y2, *y3, *y4, *dens;
  
  N = n[0];
  y1 = malloc(sizeof(double)*N);
  y2 = malloc(sizeof(double)*N);
  y3 = malloc(sizeof(double)*N);
  y4 = malloc(sizeof(double)*N);
  dens = malloc(sizeof(double)*N);
  sum1 = 0.0;
 
  for (i = 1; i <= N; i++)
  {
    for (j = i; j <= N; j++)
    {
      y1[j - 1] = x1[i - 1] - x1[j - 1];
      y2[j - 1] = x2[i - 1] - x2[j - 1];
      y3[j - 1] = x3[i - 1] - x3[j - 1];
      y4[j - 1] = x4[i - 1] - x4[j - 1];
    }
    dmvnorm_4d(y1, y2, y3, y4, mu, visigma, detsigma, n, dens);
    for (k = i; k <= N; k++)
      sum1 = sum1 + dens[k - 1];
  }

  sum[0] = sum1;
  free(y1);
  free(y2);
  free(y3);
  free(y4);
  free(dens);
}


/*************************************************************************
* 2nd order partial derivative of 4-dim normal density 
*************************************************************************/

void dmvnormd2_4d(double *x1, double *x2, double *x3, double *x4, double *vsigma, 
		  int *r, int *n, double *derivt)
{
  double sigma11, sigma12, sigma13, sigma14, sigma22, sigma23, sigma24, 
    sigma33, sigma34, sigma44, y1, y2, y3, y4;
  int i;
    
  sigma11 = sqrt(vsigma[0]);
  sigma12 = vsigma[1];
  sigma13 = vsigma[2];
  sigma14 = vsigma[3];
  sigma22 = sqrt(vsigma[4]);
  sigma23 = vsigma[5];
  sigma24 = vsigma[6];
  sigma33 = sqrt(vsigma[7]);
  sigma34 = vsigma[8];
  sigma44 = sqrt(vsigma[9]);


  /* second order derivatives */
    
  for (i = 1; i <= n[0]; i++)
  {
    y1 = x1[i - 1];
    y2 = x2[i - 1];
    y3 = x3[i - 1];
    y4 = x4[i - 1];

  if ((r[0] == 2) && (r[1] == 0) && (r[2] == 0) && (r[3] == 0)) 	  
  {  
    derivt[i - 1] =
(-(1/(pow(M_E,(pow(y1,2)/pow(sigma11,2) + pow(y2,2)/pow(sigma22,2) + 
              pow(y3,2)/pow(sigma33,2) + pow(y4,2)/pow(sigma44,2))/2.)*
          pow(sigma11,2))) + pow(y1,2)/
      (pow(M_E,(pow(y1,2)/pow(sigma11,2) + pow(y2,2)/pow(sigma22,2) + 
            pow(y3,2)/pow(sigma33,2) + pow(y4,2)/pow(sigma44,2))/2.)*pow(sigma11,4)
       ))/(4.*pow(M_PI,2)*sigma11*sigma22*sigma33*sigma44);
  }
  else if ((r[0] == 0) && (r[1] == 2) && (r[2] == 0) && (r[3] == 0))
  {
    derivt[i - 1] =
      (-(1/(pow(M_E,(pow(y1,2)/pow(sigma11,2) + pow(y2,2)/pow(sigma22,2) + 
              pow(y3,2)/pow(sigma33,2) + pow(y4,2)/pow(sigma44,2))/2.)*
          pow(sigma22,2))) + pow(y2,2)/
      (pow(M_E,(pow(y1,2)/pow(sigma11,2) + pow(y2,2)/pow(sigma22,2) + 
            pow(y3,2)/pow(sigma33,2) + pow(y4,2)/pow(sigma44,2))/2.)*pow(sigma22,4)
       ))/(4.*pow(M_PI,2)*sigma11*sigma22*sigma33*sigma44);
  }
  else if ((r[0] == 0) && (r[1] == 0) && (r[2] == 2) && (r[3] == 0))
  {
    derivt[i - 1] = 
      (-(1/(pow(M_E,(pow(y1,2)/pow(sigma11,2) + pow(y2,2)/pow(sigma22,2) + 
              pow(y3,2)/pow(sigma33,2) + pow(y4,2)/pow(sigma44,2))/2.)*
          pow(sigma33,2))) + pow(y3,2)/
      (pow(M_E,(pow(y1,2)/pow(sigma11,2) + pow(y2,2)/pow(sigma22,2) + 
            pow(y3,2)/pow(sigma33,2) + pow(y4,2)/pow(sigma44,2))/2.)*pow(sigma33,4)
       ))/(4.*pow(M_PI,2)*sigma11*sigma22*sigma33*sigma44);
  }
  else if ((r[0] == 0) && (r[1] == 0) && (r[2] == 0) && (r[3] == 2))
  {    
    derivt[i - 1] =
	(-(1/(pow(M_E,(pow(y1,2)/pow(sigma11,2) + pow(y2,2)/pow(sigma22,2) + 
              pow(y3,2)/pow(sigma33,2) + pow(y4,2)/pow(sigma44,2))/2.)*
          pow(sigma44,2))) + pow(y4,2)/
      (pow(M_E,(pow(y1,2)/pow(sigma11,2) + pow(y2,2)/pow(sigma22,2) + 
            pow(y3,2)/pow(sigma33,2) + pow(y4,2)/pow(sigma44,2))/2.)*pow(sigma44,4)
       ))/(4.*pow(M_PI,2)*sigma11*sigma22*sigma33*sigma44);
  }
  else if ((r[0] == 1) && (r[1] == 1) && (r[2] == 0) && (r[3] == 0))
  {
    derivt[i - 1] = (y1*y2)/(4.*pow(M_E,(pow(y1,2)/pow(sigma11,2) + 
         pow(y2,2)/pow(sigma22,2) + pow(y3,2)/pow(sigma33,2) + 
         pow(y4,2)/pow(sigma44,2))/2.)*pow(M_PI,2)*pow(sigma11,3)*
			       pow(sigma22,3)*sigma33*sigma44);
  }
  else if ((r[0] == 1) && (r[1] == 0) && (r[2] == 1) && (r[3] == 0))
  {
    derivt[i - 1] = (y1*y3)/(4.*pow(M_E,(pow(y1,2)/pow(sigma11,2) + 
         pow(y2,2)/pow(sigma22,2) + pow(y3,2)/pow(sigma33,2) + 
         pow(y4,2)/pow(sigma44,2))/2.)*pow(M_PI,2)*pow(sigma11,3)*sigma22*
			     pow(sigma33,3)*sigma44);
  }
  else if ((r[0] == 1) && (r[1] == 0) && (r[2] == 0) && (r[3] == 1))
  {
    derivt[i - 1] = (y1*y4)/(4.*pow(M_E,(pow(y1,2)/pow(sigma11,2) + 
         pow(y2,2)/pow(sigma22,2) + pow(y3,2)/pow(sigma33,2) + 
         pow(y4,2)/pow(sigma44,2))/2.)*pow(M_PI,2)*pow(sigma11,3)*sigma22*
			     sigma33*pow(sigma44,3));
  }
  else if ((r[0] == 0) && (r[1] == 1) && (r[2] == 1) && (r[3] == 0))
  {    
    derivt[i - 1] = (y2*y3)/(4.*pow(M_E,(pow(y1,2)/pow(sigma11,2) + 
         pow(y2,2)/pow(sigma22,2) + pow(y3,2)/pow(sigma33,2) + 
         pow(y4,2)/pow(sigma44,2))/2.)*pow(M_PI,2)*sigma11*pow(sigma22,3)*
			     pow(sigma33,3)*sigma44);
  }
  else if ((r[0] == 0) && (r[1] == 1) && (r[2] == 0) && (r[3] == 1))		     
  {    
    derivt[i - 1] = (y2*y4)/(4.*pow(M_E,(pow(y1,2)/pow(sigma11,2) + 
         pow(y2,2)/pow(sigma22,2) + pow(y3,2)/pow(sigma33,2) + 
         pow(y4,2)/pow(sigma44,2))/2.)*pow(M_PI,2)*sigma11*pow(sigma22,3)*
			     sigma33*pow(sigma44,3));
  }
  else if ((r[0] == 0) && (r[1] == 0) && (r[2] == 1) && (r[3] == 1))
  {    
    derivt[i - 1] = (y3*y4)/(4.*pow(M_E,(pow(y1,2)/pow(sigma11,2) + 
         pow(y2,2)/pow(sigma22,2) + pow(y3,2)/pow(sigma33,2) + 
         pow(y4,2)/pow(sigma44,2))/2.)*pow(M_PI,2)*sigma11*sigma22*
			     pow(sigma33,3)*pow(sigma44,3));
  }
  }
}

/*************************************************************************
* 5-dim normal density 
*************************************************************************/

void dmvnorm_5d(double *x1, double *x2, double *x3, double *x4, double *x5, 
		double *mu, 
		double *visigma, double *detsigma, int *n, double *dens)
{
  double norm;
  double xmu[5];
  int j, d;
  
  d = 5;
  norm = 1/sqrt(pow(2*M_PI, d) * detsigma[0]);
  for (j = 1; j <= n[0]; j++)
  {
    xmu[0] = x1[j - 1] - mu[0];
    xmu[1] = x2[j - 1] - mu[1];
    xmu[2] = x3[j - 1] - mu[2]; 
    xmu[3] = x4[j - 1] - mu[3];
    xmu[4] = x5[j - 1] - mu[4]; 
    dens[j - 1] = norm * exp(-0.5 * mult(xmu, visigma, xmu, d));
  }
}

/*************************************************************************
* 5-dim normal density - double sum 
*************************************************************************/

void dmvnorm_5d_sum(double *x1, double *x2, double *x3, double *x4, double *x5, 
		    double *visigma, double *detsigma, int *n, double *sum)
{
  int i, j, k, N;
  double mu[] = {0.0, 0.0, 0.0, 0.0, 0.0};
  double sum1;
  double *y1, *y2, *y3, *y4, *y5, *dens;
  
  N = n[0];
  y1 = malloc(sizeof(double)*N);
  y2 = malloc(sizeof(double)*N);
  y3 = malloc(sizeof(double)*N);
  y4 = malloc(sizeof(double)*N);
  y5 = malloc(sizeof(double)*N);
  dens = malloc(sizeof(double)*N);
  sum1 = 0.0;
 
  for (i = 1; i <= N; i++)
  {
    for (j = i; j <= N; j++)
    {
      y1[j - 1] = x1[i - 1] - x1[j - 1];
      y2[j - 1] = x2[i - 1] - x2[j - 1];
      y3[j - 1] = x3[i - 1] - x3[j - 1];
      y4[j - 1] = x4[i - 1] - x4[j - 1];
      y5[j - 1] = x5[i - 1] - x5[j - 1];
    }
    dmvnorm_5d(y1, y2, y3, y4, y5, mu, visigma, detsigma, n, dens);
    for (k = i; k <= N; k++)
      sum1 = sum1 + dens[k - 1];
  }

  sum[0] = sum1;
  free(y1);
  free(y2);
  free(y3);
  free(y4);
  free(y5);
  free(dens);
}





/*************************************************************************
* 6-dim normal density 
*************************************************************************/

void dmvnorm_6d(double *x1, double *x2, double *x3, double *x4, double *x5, 
		double *x6, double *mu, 
		double *visigma, double *detsigma, int *n, double *dens)
{
  double norm;
  double xmu[6];
  int j, d;
  
  d = 6;
  norm = 1/sqrt(pow(2*M_PI, d) * detsigma[0]);
  for (j = 1; j <= n[0]; j++)
  {
    xmu[0] = x1[j - 1] - mu[0];
    xmu[1] = x2[j - 1] - mu[1];
    xmu[2] = x3[j - 1] - mu[2]; 
    xmu[3] = x4[j - 1] - mu[3];
    xmu[4] = x5[j - 1] - mu[4]; 
    xmu[5] = x6[j - 1] - mu[5];
    dens[j - 1] = norm * exp(-0.5 * mult(xmu, visigma, xmu, d));
  }
}

/*************************************************************************
* 6-dim normal density - double sum 
*************************************************************************/

void dmvnorm_6d_sum(double *x1, double *x2, double *x3, double *x4, double *x5, 
		    double *x6, double *visigma, 
		    double *detsigma, int *n, double *sum)
{
  int i, j, k, N;
  double mu[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
  double sum1;
  double *y1, *y2, *y3, *y4, *y5, *y6, *dens;
  
  N = n[0];
  y1 = malloc(sizeof(double)*N);
  y2 = malloc(sizeof(double)*N);
  y3 = malloc(sizeof(double)*N);
  y4 = malloc(sizeof(double)*N);
  y5 = malloc(sizeof(double)*N);
  y6 = malloc(sizeof(double)*N);
  dens = malloc(sizeof(double)*N);
  sum1 = 0.0;
 
  for (i = 1; i <= N; i++)
  {
    for (j = i; j <= N; j++)
    {
      y1[j - 1] = x1[i - 1] - x1[j - 1];
      y2[j - 1] = x2[i - 1] - x2[j - 1];
      y3[j - 1] = x3[i - 1] - x3[j - 1];
      y4[j - 1] = x4[i - 1] - x4[j - 1];
      y5[j - 1] = x5[i - 1] - x5[j - 1];
      y6[j - 1] = x6[i - 1] - x6[j - 1];
    }
    dmvnorm_6d(y1, y2, y3, y4, y5, y6, mu, visigma, detsigma, n, dens);
    for (k = i; k <= N; k++)
      sum1 = sum1 + dens[k - 1];
  }

  sum[0] = sum1;
  free(y1);
  free(y2);
  free(y3);
  free(y4);
  free(y5);
  free(y6);
  free(dens);
}