ks.c

来自「r软件 另一款可以计算核估计的软件包 需安装r软件」· C语言 代码 · 共 2,039 行 · 第 1/5 页

        }
    
    else if ((r[0] == 4) && (r[1] == 2))
      for (i = 1; i <= n[0]; i++)
        {
	  x = x1[i - 1];
	  y = x2[i - 1];
	  derivt[i - 1] =
	    (pow(M_E,(pow(sigma22,2)*pow(x,2) - 
		      2*rho*sigma11*sigma22*x*y + 
		      pow(sigma11,2)*pow(y,2))/
		 (2.*(-1 + pow(rho,2))*pow(sigma11,2)*
		  pow(sigma22,2)))*
	     (pow(rho,2)*pow(sigma22,6)*pow(x,6) - 
	      2*rho*(1 + 2*pow(rho,2))*sigma11*pow(sigma22,5)*
	      pow(x,5)*y - 4*rho*pow(sigma11,3)*
	      pow(sigma22,3)*pow(x,3)*y*
	      ((-6 - 3*pow(rho,2) + 9*pow(rho,4))*
	       pow(sigma22,2) + 
	       (1 + 3*pow(rho,2) + pow(rho,4))*pow(y,2)) + 
	      pow(sigma11,2)*pow(sigma22,4)*pow(x,4)*
	      ((-1 - 13*pow(rho,2) + 14*pow(rho,4))*
	       pow(sigma22,2) + 
	       (1 + 8*pow(rho,2) + 6*pow(rho,4))*pow(y,2))
	      - 2*rho*pow(sigma11,5)*sigma22*x*y*
	      (3*pow(-1 + pow(rho,2),2)*(7 + 8*pow(rho,2))*
	       pow(sigma22,4) + 
	       2*(-3 - 7*pow(rho,2) + 8*pow(rho,4) + 
		  2*pow(rho,6))*pow(sigma22,2)*pow(y,2) + 
	       pow(rho,2)*(2 + pow(rho,2))*pow(y,4)) + 
	      pow(sigma11,4)*pow(sigma22,2)*pow(x,2)*
	      (3*pow(-1 + pow(rho,2),2)*(2 + 13*pow(rho,2))*
	       pow(sigma22,4) + 
	       6*(-1 - 8*pow(rho,2) + 4*pow(rho,4) + 
		  5*pow(rho,6))*pow(sigma22,2)*pow(y,2) + 
	       pow(rho,2)*(6 + 8*pow(rho,2) + pow(rho,4))*
	       pow(y,4)) + 
	      pow(sigma11,6)*
	      (3*pow(-1 + pow(rho,2),3)*(1 + 4*pow(rho,2))*
	       pow(sigma22,6) + 
	       3*pow(-1 + pow(rho,2),2)*
	       (1 + 10*pow(rho,2) + 4*pow(rho,4))*
	       pow(sigma22,4)*pow(y,2) + 
	       3*pow(rho,2)*
	       (-2 - pow(rho,2) + 3*pow(rho,4))*
	       pow(sigma22,2)*pow(y,4) + 
	       pow(rho,4)*pow(y,6))))/
	    (2.*M_PI*pow(1 - pow(rho,2),6.5)*pow(sigma11,11)*
	     pow(sigma22,9));
	}
    
    else if ((r[0] == 3) && (r[1] == 3))
      for (i = 1; i <= n[0]; i++)
        {
	  x = x1[i - 1];
	  y = x2[i - 1];
	  derivt[i - 1] = 
	    (pow(M_E,(pow(sigma22,2)*pow(x,2) - 
		      2*rho*sigma11*sigma22*x*y + 
		      pow(sigma11,2)*pow(y,2))/
		 (2.*(-1 + pow(rho,2))*pow(sigma11,2)*
		  pow(sigma22,2)))*
	     (-6*pow(rho,9)*pow(sigma11,6)*pow(sigma22,6) + 
	      18*pow(rho,8)*pow(sigma11,5)*pow(sigma22,5)*x*
	      y + pow(sigma11,3)*pow(sigma22,3)*x*
	      (3*pow(sigma11,2) - pow(x,2))*y*
	      (3*pow(sigma22,2) - pow(y,2)) + 
	      9*pow(rho,7)*pow(sigma11,4)*pow(sigma22,4)*
	      (pow(sigma11,2)*
	       (pow(sigma22,2) - 3*pow(y,2)) - 
	       pow(x,2)*(3*pow(sigma22,2) + pow(y,2))) + 
	      pow(rho,6)*pow(sigma11,3)*pow(sigma22,3)*x*y*
	      (pow(x,2)*(21*pow(sigma22,2) + pow(y,2)) + 
	       3*pow(sigma11,2)*
	       (9*pow(sigma22,2) + 7*pow(y,2))) + 
	      3*pow(rho,5)*pow(sigma11,2)*pow(sigma22,2)*
        (-(pow(sigma22,2)*pow(x,4)*
	   (4*pow(sigma22,2) + pow(y,2))) + 
	 pow(sigma11,4)*
	 (3*pow(sigma22,4) + 
	  12*pow(sigma22,2)*pow(y,2) - 4*pow(y,4))
	 + pow(sigma11,2)*pow(x,2)*
	 (12*pow(sigma22,4) - 
	  15*pow(sigma22,2)*pow(y,2) - pow(y,4)))
	      + 3*pow(rho,2)*sigma11*sigma22*x*y*
	      (pow(sigma22,4)*pow(x,4) + 
	       pow(sigma11,2)*pow(sigma22,2)*pow(x,2)*
	       (-11*pow(sigma22,2) + 3*pow(y,2)) + 
	       pow(sigma11,4)*
	       (15*pow(sigma22,4) - 
		11*pow(sigma22,2)*pow(y,2) + pow(y,4)))
	      + 3*rho*pow(sigma11,2)*pow(sigma22,2)*
	      (pow(sigma22,2)*pow(x,4)*
	       (pow(sigma22,2) - pow(y,2)) - 
	       pow(sigma11,2)*pow(x,2)*
	       (6*pow(sigma22,4) - 
		9*pow(sigma22,2)*pow(y,2) + pow(y,4)) + 
	       pow(sigma11,4)*
	       (3*pow(sigma22,4) - 
		6*pow(sigma22,2)*pow(y,2) + pow(y,4))) + 
	      3*pow(rho,4)*sigma11*sigma22*x*y*
	      (pow(sigma22,4)*pow(x,4) + 
	       pow(sigma11,2)*pow(sigma22,2)*pow(x,2)*
	       (5*pow(sigma22,2) + 3*pow(y,2)) + 
	       pow(sigma11,4)*
	       (-33*pow(sigma22,4) + 
		5*pow(sigma22,2)*pow(y,2) + pow(y,4))) - 
	      pow(rho,3)*(pow(sigma22,6)*pow(x,6) + 
			  9*pow(sigma11,2)*pow(sigma22,4)*pow(x,4)*
			  (-pow(sigma22,2) + pow(y,2)) - 
			  9*pow(sigma11,4)*pow(sigma22,2)*pow(x,2)*
			  (pow(sigma22,4) + 
			   3*pow(sigma22,2)*pow(y,2) - pow(y,4)) + 
			  pow(sigma11,6)*
			  (21*pow(sigma22,6) - 
			   9*pow(sigma22,4)*pow(y,2) - 
			   9*pow(sigma22,2)*pow(y,4) + pow(y,6)))))/
	    (2.*M_PI*pow(1 - pow(rho,2),6.5)*pow(sigma11,10)*
	     pow(sigma22,10));
        }
    
    else if ((r[0] == 2) && (r[1] == 4))
      for (i = 1; i <= n[0]; i++)
        {
	  x = x1[i - 1];
	  y = x2[i - 1];  
	  derivt[i - 1] =
	    (pow(M_E,(pow(sigma22,2)*pow(x,2) - 
		      2*rho*sigma11*sigma22*x*y + 
		      pow(sigma11,2)*pow(y,2))/
		 (2.*(-1 + pow(rho,2))*pow(sigma11,2)*
		  pow(sigma22,2)))*
	     (pow(rho,4)*pow(sigma22,6)*pow(x,6) - 
	      2*pow(rho,3)*(2 + pow(rho,2))*sigma11*
	      pow(sigma22,5)*pow(x,5)*y - 
	      4*rho*pow(sigma11,3)*pow(sigma22,3)*pow(x,3)*y*
	      ((-3 - 7*pow(rho,2) + 8*pow(rho,4) + 
		2*pow(rho,6))*pow(sigma22,2) + 
	       (1 + 3*pow(rho,2) + pow(rho,4))*pow(y,2)) + 
	      pow(rho,2)*pow(sigma11,2)*pow(sigma22,4)*
	      pow(x,4)*((-6 - 3*pow(rho,2) + 9*pow(rho,4))*
			pow(sigma22,2) + 
			(6 + 8*pow(rho,2) + pow(rho,4))*pow(y,2)) - 
	      2*rho*pow(sigma11,5)*sigma22*x*y*
	      (3*pow(-1 + pow(rho,2),2)*(7 + 8*pow(rho,2))*
	       pow(sigma22,4) + 
	       6*(-2 - pow(rho,2) + 3*pow(rho,4))*
	       pow(sigma22,2)*pow(y,2) + 
	       (1 + 2*pow(rho,2))*pow(y,4)) + 
	      pow(sigma11,4)*pow(sigma22,2)*pow(x,2)*
	      (3*pow(-1 + pow(rho,2),2)*
	       (1 + 10*pow(rho,2) + 4*pow(rho,4))*
	       pow(sigma22,4) + 
	       6*(-1 - 8*pow(rho,2) + 4*pow(rho,4) + 
		  5*pow(rho,6))*pow(sigma22,2)*pow(y,2) + 
	       (1 + 8*pow(rho,2) + 6*pow(rho,4))*pow(y,4))
	      + pow(sigma11,6)*
	      (3*pow(-1 + pow(rho,2),3)*(1 + 4*pow(rho,2))*
	       pow(sigma22,6) + 
	       3*pow(-1 + pow(rho,2),2)*
	       (2 + 13*pow(rho,2))*pow(sigma22,4)*pow(y,2)
	       + (-1 - 13*pow(rho,2) + 14*pow(rho,4))*
	       pow(sigma22,2)*pow(y,4) + 
	       pow(rho,2)*pow(y,6))))/
	    (2.*M_PI*pow(1 - pow(rho,2),6.5)*pow(sigma11,9)*
	     pow(sigma22,11));
        }

    else if ((r[0] == 1) && (r[1] == 5))
        for (i = 1; i <= n[0]; i++)
        {
	  x = x1[i - 1];
	  y = x2[i - 1];
	  derivt[i - 1] =
	    (pow(M_E,(pow(sigma22,2)*pow(x,2) - 
		      2*rho*sigma11*sigma22*x*y + 
		      pow(sigma11,2)*pow(y,2))/
		 (2.*(-1 + pow(rho,2))*pow(sigma11,2)*
		  pow(sigma22,2)))*
	     (-5*pow(rho,7)*pow(sigma11,2)*pow(sigma22,6)*
        (3*pow(sigma11,4) + 
	 6*pow(sigma11,2)*pow(x,2) + pow(x,4)) + 
	      pow(rho,6)*sigma11*pow(sigma22,5)*x*
	      (75*pow(sigma11,4) + 
	       30*pow(sigma11,2)*pow(x,2) + pow(x,4))*y + 
	      pow(sigma11,5)*sigma22*x*y*
	      (15*pow(sigma22,4) - 
	       10*pow(sigma22,2)*pow(y,2) + pow(y,4)) + 
	      pow(rho,5)*pow(sigma22,4)*
	      (-(pow(sigma22,2)*pow(x,6)) + 
	       15*pow(sigma11,4)*pow(x,2)*
	       (3*pow(sigma22,2) - 4*pow(y,2)) + 
	       45*pow(sigma11,6)*
	       (pow(sigma22,2) - pow(y,2)) - 
	       5*pow(sigma11,2)*pow(x,4)*
	       (pow(sigma22,2) + pow(y,2))) + 
	      5*pow(rho,4)*sigma11*pow(sigma22,3)*x*y*
	      (pow(sigma22,2)*pow(x,4) + 
	       2*pow(sigma11,2)*pow(x,2)*pow(y,2) + 
	       pow(sigma11,4)*
	       (-27*pow(sigma22,2) + 10*pow(y,2))) + 
	      5*pow(rho,2)*pow(sigma11,3)*sigma22*x*y*
        (2*pow(sigma22,2)*pow(x,2)*
	 (-3*pow(sigma22,2) + pow(y,2)) + 
	 pow(sigma11,2)*
	 (9*pow(sigma22,4) - 
	  8*pow(sigma22,2)*pow(y,2) + pow(y,4))) - 
	      5*pow(rho,3)*pow(sigma11,2)*pow(sigma22,2)*
	      (2*pow(sigma11,2)*pow(x,2)*pow(y,2)*
	       (-3*pow(sigma22,2) + pow(y,2)) + 
	       2*pow(sigma22,2)*pow(x,4)*
	       (-pow(sigma22,2) + pow(y,2)) + 
	       3*pow(sigma11,4)*
	       (3*pow(sigma22,4) - 
		6*pow(sigma22,2)*pow(y,2) + pow(y,4))) + 
	      rho*pow(sigma11,4)*
	      (-5*pow(sigma22,2)*pow(x,2)*
	       (3*pow(sigma22,4) - 
		6*pow(sigma22,2)*pow(y,2) + pow(y,4)) + 
	       pow(sigma11,2)*
	       (15*pow(sigma22,6) - 
		45*pow(sigma22,4)*pow(y,2) + 
		15*pow(sigma22,2)*pow(y,4) - pow(y,6)))))
	    /(2.*M_PI*pow(1 - pow(rho,2),6.5)*pow(sigma11,8)*
	      pow(sigma22,12));
        }

    else if ((r[0] == 0) && (r[1] == 6))
        for (i = 1; i <= n[0]; i++)
        {
	  x = x1[i - 1];
	  y = x2[i - 1];
	  
	  derivt[i - 1] =
	    (pow(M_E,(pow(sigma22,2)*pow(x,2) - 
		      2*rho*sigma11*sigma22*x*y + 
		      pow(sigma11,2)*pow(y,2))/
     (2.*(-1 + pow(rho,2))*pow(sigma11,2)*
      pow(sigma22,2)))*
	     (pow(rho,6)*pow(sigma22,6)*pow(x,6) - 
	      6*pow(rho,5)*sigma11*pow(sigma22,5)*pow(x,5)*
	      y + 15*pow(rho,4)*pow(sigma11,2)*
	      pow(sigma22,4)*pow(x,4)*
	      ((-1 + pow(rho,2))*pow(sigma22,2) + pow(y,2))
	      - 20*pow(rho,3)*pow(sigma11,3)*pow(sigma22,3)*
	      pow(x,3)*y*(3*(-1 + pow(rho,2))*
			  pow(sigma22,2) + pow(y,2)) + 
	      15*pow(rho,2)*pow(sigma11,4)*pow(sigma22,2)*
	      pow(x,2)*(3*pow(-1 + pow(rho,2),2)*
			pow(sigma22,4) + 
			6*(-1 + pow(rho,2))*pow(sigma22,2)*
			pow(y,2) + pow(y,4)) - 
	      6*rho*pow(sigma11,5)*sigma22*x*y*
	      (15*pow(-1 + pow(rho,2),2)*pow(sigma22,4) + 
	       10*(-1 + pow(rho,2))*pow(sigma22,2)*
	       pow(y,2) + pow(y,4)) +pow(sigma11,6)*
	      (15*pow(-1 + pow(rho,2),3)*pow(sigma22,6) + 
	       45*pow(-1 + pow(rho,2),2)*pow(sigma22,4)*pow(y,2) + 
	       15*(-1 + pow(rho,2))*pow(sigma22,2)*
	       pow(y,4) + pow(y,6))))/
	    (2.*M_PI*pow(1 - pow(rho,2),6.5)*pow(sigma11,7)*
	     pow(sigma22,13));
        }
}

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

void dmvnorm_2d_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 sum1;
    double *x, *y, *dens;

    N = n[0];
    x = malloc(sizeof(double)*N);
    y = malloc(sizeof(double)*N);
    dens = 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];
        }
        dmvnorm_2d(x, y, mu, visigma, detsigma, n, dens);
        for (k = i; k <= N; k++)
            sum1 = sum1 + dens[k - 1];
    }

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





/*************************************************************************
* Double sum of normal density values for Abramson's selector - bivariate
*
* Parameters
* x1 - x values
* x2 - y values
* visigma - vec of inverse Sigma
* n - number of values
* sum - contains sum of density values
*************************************************************************/

void dmvnorm_2d_sum_ab(double *x1, double *x2, double *sigma1, double *sigma2,
		       int *n, int *which, double *sum)
{
  int i, j, N, d, w;
  double visigma[] = {0.0, 0.0, 0.0, 0.0};
  double sum1, norm, detsigma;
  double *x;
  
  d = 2;
  N = n[0];
  w = which[0];
  sum1 = 0.0;
  x = malloc(sizeof(double)*d);

  for (i = 1; i <= N; i++)
  {
    for (j = 1; j <= N; j++)
    {
      x[0] = x1[i - 1] - x1[j - 1];
      x[1] = x2[i - 1] - x2[j - 1];
      if (w==1)
      {
	visigma[0] = 1/(sigma1[i - 1] + sigma1[j - 1]);
	visigma[3] = 1/(sigma2[i - 1] + sigma2[j - 1]);
      }
      else 
      {
	visigma[0] = 1/sigma1[j - 1];
	visigma[3] = 1/sigma2[j - 1];
      }
      detsigma = 1/(visigma[0]*visigma[3]); 
      norm = 1/sqrt(pow(2*M_PI, d) * detsigma); 
      if ((w==0) & (i==j))
	;
      else 
	sum1 = sum1 + norm * exp(-0.5 * mult(x, visigma, x, d));
    }
  }

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


/*************************************************************************
* Double sum of normal density values for re-clustered selector - bivariate
*
* Parameters
* x1 - x values
* x2 - y values
* visigma - vec of inverse Sigma
* n - number of values
* sum - contains sum of density values
*************************************************************************/

void dmvnorm_2d_sum_clust(double *x1, double *x2, double *y1, double *y2, 
			 double *visigma, double *detsigma, int *nx, int *ny, 
			 double *sum)
{
  int i, j, k, Nx, Ny;
  double mu[] = {0.0, 0.0};
  double sum1;
  double *x, *y, *dens;
    
  Nx = nx[0];
  Ny = ny[0];  
  x = malloc(sizeof(double)*Ny);
  y = malloc(sizeof(double)*Ny);
  dens = malloc(sizeof(double)*Ny);
  sum1 = 0.0;
  
  for (i = 1; i <= Nx; i++)
  {
    for (j = 1; j <= Ny; j++)
    {
      x[j - 1] = x1[i - 1] - y1[j - 1];
      y[j - 1] = x2[i - 1] - y2[j - 1];
    }
    dmvnorm_2d(x, y, mu, visigma, detsigma, ny, dens);
    for (k = 1; k <= Ny; k++)
      sum1 = sum1 + dens[k - 1];
  }

  sum[0] = sum1;
  free(x);