ks.c

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

	    6*pow(sigma11,2)*pow(sigma22,2)*pow(x,2)*
	    ((-1 + pow(rho,2))*pow(sigma22,2) + 
	     pow(rho,2)*pow(y,2)) - 4*rho*pow(sigma11,3)*sigma22*x*y*
	    (3*(-1 + pow(rho,2))*pow(sigma22,2) + 
	     pow(rho,2)*pow(y,2)) + pow(sigma11,4)*
	    (3*pow(-1 + pow(rho,2),2)*pow(sigma22,4) + 
	     6*pow(rho,2)*(-1 + pow(rho,2))*pow(sigma22,2)*pow(y,2) + 
	     pow(rho,4)*pow(y,4))))/
	  (2.*M_PI*pow(1 - pow(rho,2),4.5)*pow(sigma11,9)*
	       pow(sigma22,5));
      }

  else if ((r[0] == 3) && (r[1] == 1))
    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)))*
	   (sigma11*pow(sigma22,3)*x*
	    (-3*pow(sigma11,2) + pow(x,2))*y - 
	    3*pow(rho,5)*pow(sigma11,4)*pow(sigma22,2)*
	    (pow(sigma22,2) + pow(y,2)) + 
	    pow(rho,4)*pow(sigma11,3)*sigma22*x*y*
	    (9*pow(sigma22,2) + pow(y,2)) + 
	    3*pow(rho,2)*sigma11*sigma22*x*y*
	    (pow(sigma22,2)*pow(x,2) +pow(sigma11,2)*
	     (-2*pow(sigma22,2) + pow(y,2))) - 
	    rho*pow(sigma22,2)*
	    (pow(sigma22,2)*pow(x,4) + 3*pow(sigma11,4)*
	     (pow(sigma22,2) - pow(y,2)) + 
	     3*pow(sigma11,2)*pow(x,2)*
	     (-2*pow(sigma22,2) + pow(y,2))) + 
	    pow(rho,3)*pow(sigma11,2)*
	    (-3*pow(sigma22,2)*pow(x,2)*
	     (2*pow(sigma22,2) + pow(y,2)) + 
	     pow(sigma11,2)*(6*pow(sigma22,4) - pow(y,4))
	     )))/(2.*M_PI*pow(1 - pow(rho,2),4.5)*pow(sigma11,8)*
		  pow(sigma22,6));
      }
  
  else if ((r[0] == 2) && (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,4)*pow(x,4) - 
		2*rho*(1 + pow(rho,2))*sigma11*pow(sigma22,3)*
		pow(x,3)*y - 2*rho*pow(sigma11,3)*sigma22*x*y*
		(2*(-2 + pow(rho,2) + pow(rho,4))*
		 pow(sigma22,2) + (1 + pow(rho,2))*pow(y,2))
		+ pow(sigma11,2)*pow(sigma22,2)*pow(x,2)*
		((-1 - 4*pow(rho,2) + 5*pow(rho,4))*
		 pow(sigma22,2) + 
		 (1 + 4*pow(rho,2) + pow(rho,4))*pow(y,2)) + 
		pow(sigma11,4)*
		((1 - 3*pow(rho,4) + 2*pow(rho,6))*
		 pow(sigma22,4) + 
		 (-1 - 4*pow(rho,2) + 5*pow(rho,4))*
		 pow(sigma22,2)*pow(y,2) + 
	 pow(rho,2)*pow(y,4))))/
	      (2.*M_PI*pow(1 - pow(rho,2),4.5)*pow(sigma11,7)*
	       pow(sigma22,7));
      }

  else if ((r[0] == 1) && (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)))*
	   (-3*pow(rho,5)*pow(sigma11,2)*pow(sigma22,4)*
	    (pow(sigma11,2) + pow(x,2)) + 
	    pow(rho,4)*sigma11*pow(sigma22,3)*x*
	    (9*pow(sigma11,2) + pow(x,2))*y + 
	    pow(sigma11,3)*sigma22*x*y*
	    (-3*pow(sigma22,2) + pow(y,2)) + 
	    pow(rho,3)*pow(sigma22,2)*
	    (-(pow(sigma22,2)*pow(x,4)) - 
	     3*pow(sigma11,2)*pow(x,2)*pow(y,2) + 
	     6*pow(sigma11,4)*(pow(sigma22,2) - pow(y,2))
	     ) + 3*pow(rho,2)*sigma11*sigma22*x*y*
	    (pow(sigma22,2)*pow(x,2) + 
	     pow(sigma11,2)*
	     (-2*pow(sigma22,2) + pow(y,2))) - 
	    rho*pow(sigma11,2)*
	    (3*pow(sigma22,2)*pow(x,2)*
	     (-pow(sigma22,2) + pow(y,2)) + 
	     pow(sigma11,2)*
	     (3*pow(sigma22,4) - 
	      6*pow(sigma22,2)*pow(y,2) + pow(y,4)))))/
	  (2.*M_PI*pow(1 - pow(rho,2),4.5)*pow(sigma11,6)*
	   pow(sigma22,8));
      }
  
    else if ((r[0] == 0) && (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,4)*pow(x,4) - 
	      4*pow(rho,3)*sigma11*pow(sigma22,3)*pow(x,3)*
	      y + 6*pow(rho,2)*pow(sigma11,2)*
	      pow(sigma22,2)*pow(x,2)*
	      ((-1 + pow(rho,2))*pow(sigma22,2) + pow(y,2))
	      - 4*rho*pow(sigma11,3)*sigma22*x*y*
	      (3*(-1 + pow(rho,2))*pow(sigma22,2) +pow(y,2)) + 
	      pow(sigma11,4)*
		(3*pow(-1 + pow(rho,2),2)*pow(sigma22,4) + 
		 6*(-1 + pow(rho,2))*pow(sigma22,2)*
		 pow(y,2) + pow(y,4))))/
	    (2.*M_PI*pow(1 - pow(rho,2),4.5)*pow(sigma11,5)*pow(sigma22,9));
	}
}


/*************************************************************************
* Bivariate normal density - 5th order derivatives
*
* Parameters
* x1 - x values
* x2 - y values
* v - vec Sigma
* r - (r1, r2) partial derivative
* x - number of values
* derivt - contains density derivative values
*************************************************************************/
void dmvnormd5_2d(double *x1, double *x2, double *vsigma, int *r, int *n, 
		  double *derivt)
{
  /* fifth order derivatives */

  double sigma11, sigma22, rho, x, y;
  int i;
  
  sigma11 = sqrt(vsigma[0]);
  sigma22 = sqrt(vsigma[3]);
  rho = vsigma[1] /(sigma11 * sigma22);
    
  if ((r[0] == 5) && (r[1] == 0))
    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)))*
      (sigma22*x - rho*sigma11*y)*(240*pow(-1 + pow(rho,2),2)*pow(sigma11,4)*
         pow(sigma22,8) + 160*(-1 + pow(rho,2))*pow(sigma11,2)*pow(sigma22,6)*
         pow(sigma22*x - rho*sigma11*y,2) + 
        pow(2*pow(sigma22,2)*x - 2*rho*sigma11*sigma22*y,4)))/
	(32.*M_PI*pow(1 - pow(rho,2),5.5)*pow(sigma11,11)*pow(sigma22,10));
    }
  else if ((r[0] == 4) && (r[1] == 1))
    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)))*
      (-3*pow(-1 + pow(rho,2),2)*pow(sigma11,4)*pow(sigma22,4)*
         (rho*sigma22*x - sigma11*y) - 
        6*(-1 + pow(rho,2))*pow(sigma11,2)*pow(sigma22,2)*(rho*sigma22*x - sigma11*y)*
         pow(sigma22*x - rho*sigma11*y,2) + 
        (-(rho*sigma22*x) + sigma11*y)*pow(sigma22*x - rho*sigma11*y,4) + 
        12*rho*pow(-1 + pow(rho,2),2)*pow(sigma11,4)*pow(sigma22,4)*
         (-(sigma22*x) + rho*sigma11*y) + 
        4*rho*(-1 + pow(rho,2))*pow(sigma11,2)*pow(sigma22,2)*
         pow(-(sigma22*x) + rho*sigma11*y,3)))/
	(2.*M_PI*pow(1 - pow(rho,2),5.5)*pow(sigma11,10)*pow(sigma22,7));
    }
  else if ((r[0] == 3) && (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,5)*pow(x,5) - 
        rho*(2 + 3*pow(rho,2))*sigma11*pow(sigma22,4)*pow(x,4)*y - 
        rho*pow(sigma11,3)*pow(sigma22,2)*pow(x,2)*y*
         (15*(-1 + pow(rho,4))*pow(sigma22,2) + 
           (3 + 6*pow(rho,2) + pow(rho,4))*pow(y,2)) + 
        pow(sigma11,2)*pow(sigma22,3)*pow(x,3)*
         ((-1 - 8*pow(rho,2) + 9*pow(rho,4))*pow(sigma22,2) + 
           (1 + 6*pow(rho,2) + 3*pow(rho,4))*pow(y,2)) - 
        rho*pow(sigma11,5)*y*(3*pow(-1 + pow(rho,2),2)*(3 + 2*pow(rho,2))*
            pow(sigma22,4) + (-3 - 4*pow(rho,2) + 7*pow(rho,4))*pow(sigma22,2)*
            pow(y,2) + pow(rho,2)*pow(y,4)) + 
        pow(sigma11,4)*sigma22*x*(3*pow(-1 + pow(rho,2),2)*(1 + 4*pow(rho,2))*
            pow(sigma22,4) + 3*(-1 - 6*pow(rho,2) + 5*pow(rho,4) + 2*pow(rho,6))*
            pow(sigma22,2)*pow(y,2) + pow(rho,2)*(3 + 2*pow(rho,2))*pow(y,4))))/
	(2.*M_PI*pow(1 - pow(rho,2),5.5)*pow(sigma11,9)*pow(sigma22,8));
    }
  else if ((r[0] == 2) && (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,7)*pow(sigma11,4)*pow(sigma22,5)*x - 
       6*pow(rho,6)*pow(sigma11,3)*pow(sigma22,4)*(2*pow(sigma11,2) + pow(x,2))*
        y - pow(sigma11,3)*pow(sigma22,2)*(pow(sigma11,2) - pow(x,2))*y*
        (3*pow(sigma22,2) - pow(y,2)) + 
       pow(rho,4)*sigma11*pow(sigma22,2)*y*
        (-2*pow(sigma22,2)*pow(x,4) + 
          3*pow(sigma11,4)*(7*pow(sigma22,2) - 3*pow(y,2)) - 
          3*pow(sigma11,2)*pow(x,2)*(5*pow(sigma22,2) + pow(y,2))) + 
       pow(rho,5)*pow(sigma11,2)*pow(sigma22,3)*x*
        (-3*pow(sigma11,2)*(pow(sigma22,2) - 5*pow(y,2)) + 
          pow(x,2)*(7*pow(sigma22,2) + pow(y,2))) + 
       pow(rho,3)*(pow(sigma22,5)*pow(x,5) + 
          2*pow(sigma11,2)*pow(sigma22,3)*pow(x,3)*
           (-2*pow(sigma22,2) + 3*pow(y,2)) + 
          3*pow(sigma11,4)*sigma22*x*(-4*pow(sigma22,4) + pow(y,4))) - 
       pow(rho,2)*sigma11*y*(3*pow(sigma22,4)*pow(x,4) + 
          6*pow(sigma11,2)*pow(sigma22,2)*pow(x,2)*
           (-3*pow(sigma22,2) + pow(y,2)) + 
          pow(sigma11,4)*(6*pow(sigma22,4) - 8*pow(sigma22,2)*pow(y,2) + pow(y,4)))
         + rho*pow(sigma11,2)*sigma22*x*
        (3*pow(sigma22,2)*pow(x,2)*(-pow(sigma22,2) + pow(y,2)) + 
          pow(sigma11,2)*(9*pow(sigma22,4) - 15*pow(sigma22,2)*pow(y,2) + 
             2*pow(y,4)))))/
	(2.*M_PI*pow(1 - pow(rho,2),5.5)*pow(sigma11,8)*pow(sigma22,9));
    }
  else if ((r[0] == 1) && (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(sigma11,2)*(sigma22*x - rho*sigma11*y)*
          (3*pow(-1 + pow(rho,2),2)*pow(sigma22,4) + 
            (6*(-1 + pow(rho,2))*pow(sigma22,2)*pow(rho*sigma22*x - sigma11*y,2))/
             pow(sigma11,2) + pow(rho*sigma22*x - sigma11*y,4)/pow(sigma11,4))) + 
       4*rho*(1 - pow(rho,2))*pow(sigma22,2)*
        (pow(rho,3)*pow(sigma22,3)*x*(3*pow(sigma11,2) + pow(x,2)) - 
          3*pow(rho,2)*sigma11*pow(sigma22,2)*(pow(sigma11,2) + pow(x,2))*y - 
          pow(sigma11,3)*y*(-3*pow(sigma22,2) + pow(y,2)) + 
          3*rho*pow(sigma11,2)*sigma22*x*(-pow(sigma22,2) + pow(y,2)))))/
	(2.*M_PI*pow(1 - pow(rho,2),5.5)*pow(sigma11,5)*pow(sigma22,10));
    }
  else if ((r[0] == 0) && (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)))*
      (-(rho*sigma22*x) + sigma11*y)*(240*pow(-1 + pow(rho,2),2)*pow(sigma11,8)*
         pow(sigma22,4) + 160*(-1 + pow(rho,2))*pow(sigma11,6)*pow(sigma22,2)*
         pow(rho*sigma22*x - sigma11*y,2) + 
        pow(2*rho*sigma11*sigma22*x - 2*pow(sigma11,2)*y,4)))/
	(32.*M_PI*pow(1 - pow(rho,2),5.5)*pow(sigma11,10)*pow(sigma22,11));
    }
}

/*************************************************************************
* Bivariate normal density - 6th order derivatives
*
* Parameters
* x1 - x values
* x2 - y values
* v - vec Sigma
* r - (r1, r2) partial derivative
* x - number of values
* derivt - contains density derivative values
*************************************************************************/
void dmvnormd6_2d(double *x1, double *x2, double *vsigma, int *r, int *n, 
		  double *derivt)
{
  /* sixth order derivatives */

  double sigma11, sigma22, rho, x, y;
  int i;

    sigma11 = sqrt(vsigma[0]);
    sigma22 = sqrt(vsigma[3]);
    rho = vsigma[1] /(sigma11 * sigma22);
    if ((r[0] == 6) && (r[1] == 0))
      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(sigma22,6)*pow(x,6) - 
	      6*rho*sigma11*pow(sigma22,5)*pow(x,5)*y + 
	      15*pow(sigma11,2)*pow(sigma22,4)*pow(x,4)*
	      ((-1 + pow(rho,2))*pow(sigma22,2) + 
	       pow(rho,2)*pow(y,2)) - 
	      20*rho*pow(sigma11,3)*pow(sigma22,3)*pow(x,3)*
	      y*(3*(-1 + pow(rho,2))*pow(sigma22,2) + 
		 pow(rho,2)*pow(y,2)) + 
	      15*pow(sigma11,4)*pow(sigma22,2)*pow(x,2)*
	      (3*pow(-1 + pow(rho,2),2)*pow(sigma22,4) + 
	       6*pow(rho,2)*(-1 + pow(rho,2))*
	       pow(sigma22,2)*pow(y,2) + 
	       pow(rho,4)*pow(y,4)) - 
	      6*rho*pow(sigma11,5)*sigma22*x*y*
	      (15*pow(-1 + pow(rho,2),2)*pow(sigma22,4) + 
	       10*pow(rho,2)*(-1 + pow(rho,2))*
	       pow(sigma22,2)*pow(y,2) + 
	       pow(rho,4)*pow(y,4))+ pow(sigma11,6)*
	      (15*pow(-1 + pow(rho,2),3)*pow(sigma22,6) + 
	       45*pow(rho,2)*pow(-1 + pow(rho,2),2)*
	       pow(sigma22,4)*pow(y,2) + 
	       15*pow(rho,4)*(-1 + pow(rho,2))*
	       pow(sigma22,2)*pow(y,4) + 
	       pow(rho,6)*pow(y,6))))/
	    (2.*M_PI*pow(1 - pow(rho,2),6.5)*pow(sigma11,13)*
	     pow(sigma22,7));
        }
    
    else if ((r[0] == 5) && (r[1] == 1))
      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)))*
	     (sigma11*pow(sigma22,5)*x*
	      (15*pow(sigma11,4) - 
	       10*pow(sigma11,2)*pow(x,2) + pow(x,4))*y - 
	      5*pow(rho,7)*pow(sigma11,6)*pow(sigma22,2)*
	      (3*pow(sigma22,4) + 
	       6*pow(sigma22,2)*pow(y,2) + pow(y,4)) + 
	      pow(rho,6)*pow(sigma11,5)*sigma22*x*y*
	      (75*pow(sigma22,4) + 
	       30*pow(sigma22,2)*pow(y,2) + pow(y,4)) + 
	      5*pow(rho,2)*sigma11*pow(sigma22,3)*x*y*
	      (pow(sigma22,2)*pow(x,4) + 
	       pow(sigma11,4)*
	       (9*pow(sigma22,2) - 6*pow(y,2)) + 
	       2*pow(sigma11,2)*pow(x,2)*
	       (-4*pow(sigma22,2) + pow(y,2))) + 
	      rho*pow(sigma22,4)*
	      (-(pow(sigma22,2)*pow(x,6)) + 
	       15*pow(sigma11,6)*
	       (pow(sigma22,2) - pow(y,2)) + 
	       5*pow(sigma11,2)*pow(x,4)*
	       (3*pow(sigma22,2) - pow(y,2)) + 
	       15*pow(sigma11,4)*pow(x,2)*
	       (-3*pow(sigma22,2) + 2*pow(y,2))) - 
	      5*pow(rho,3)*pow(sigma11,2)*pow(sigma22,2)*
	      (pow(sigma22,2)*pow(x,4)*
	       (3*pow(sigma22,2) + 2*pow(y,2)) + 
	       pow(sigma11,4)*
	       (9*pow(sigma22,4) - 2*pow(y,4)) - 
	       2*pow(sigma11,2)*pow(x,2)*
	       (9*pow(sigma22,4) + 
		3*pow(sigma22,2)*pow(y,2) - pow(y,4))) + 
	      5*pow(rho,4)*pow(sigma11,3)*sigma22*x*y*
	      (2*pow(sigma22,2)*pow(x,2)*
	       (5*pow(sigma22,2) + pow(y,2)) + 
	       pow(sigma11,2)*
	       (-27*pow(sigma22,4) + pow(y,4))) + 
	      pow(rho,5)*pow(sigma11,4)*
	      (-5*pow(sigma22,2)*pow(x,2)*
	       (9*pow(sigma22,4) + 
		12*pow(sigma22,2)*pow(y,2) + pow(y,4)) + 
	       pow(sigma11,2)*
	       (45*pow(sigma22,6) + 
		45*pow(sigma22,4)*pow(y,2) - 
		5*pow(sigma22,2)*pow(y,4) - pow(y,6)))))/
	    (2.*M_PI*pow(1 - pow(rho,2),6.5)*pow(sigma11,12)*
	     pow(sigma22,8));