psins.r

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

#############################################################################
## Estimation of psi_r for normal scales
#############################################################################


psins.1d <- function(r, sigma)
{
  ##return(dnorm.deriv(x=0, r=r, sigma=sqrt(2)*sigma))
  
  if (r %%2 ==0)
    psins <- (-1)^(r/2)*factorial(r)/((2*sigma)^(r+1)*factorial(r/2)*pi^(1/2))
  else
    psins <- 0
    
  return(psins)  
}


psins <- function(r, Sigma, complete=FALSE, Sdr.mat)
{
  d <- ncol(Sigma)
  if (complete)
  {
    return(dmvnorm.deriv(x=rep(0,d), mu=rep(0,d), r=r, Sigma=2*Sigma, Sdr.mat=Sdr.mat))
  }
  else
  {  
    if (d==2)
      return(psins.2d(r=r, Sigma=Sigma))
    if (d==3)
      return(psins.3d(r=r, Sigma=Sigma))
    if (d==4)
      return(psins.4d(r=r, Sigma=Sigma))
    if (d==5)
      return(psins.5d(r=r, Sigma=Sigma))
    if (d==6)
      return(psins.6d(r=r, Sigma=Sigma))
  }
}


psins.2d <- function(r, Sigma)
{
  A <- chol2inv(chol(2*Sigma))
  a11 <- A[1,1]
  a22 <- A[2,2]
  a12 <- A[1,2]

  
  if ((r[1] == 0) & (r[2] == 0))   
    psi <- 1
  
  ### second order psi functionals - normal score
  else if ((r[1]==2) & (r[2]==0))
    psi <- -a11
  else if ((r[1]==1) & (r[2]==1))
    psi <- -a12
  else if ((r[1]==0) & (r[2]==2))
    psi <- -a22
  
  ### fourth order psi functionals - normal score
  else if ((r[1]==4) & (r[2]==0))
    psi <-3*a11^2
  else if ((r[1]==3) & (r[2]==1))
    psi <-3*a11*a12
  else if ((r[1]==2) & (r[2]==2))
    psi <-2*a12^2 + a11*a22
  else if ((r[1]==1) & (r[2]==3))
    psi <-3*a12*a22
  else if ((r[1]==0) & (r[2]==4))
    psi <-3*a22^2
  
  ### sixth order psi functionals - normal score 
  else if ((r[1]==6) & (r[2]==0))
    psi <- -15*a11^3
  else if ((r[1]==5) & (r[2]==1))
    psi <- -15*a11^2*a12
  else if ((r[1]==4) & (r[2]==2))
    psi <- -3*a11*(4*a12^2 + a11*a22)
  else if ((r[1]==3) & (r[2]==3))
    psi <- -6*a12^3 - 9*a11*a12*a22
  else if ((r[1]==2) & (r[2]==4))
    psi <- -3*a22*(4*a12^2 + a11*a22)
  else if ((r[1]==1) & (r[2]==5))
    psi <- -15*a12*a22^2
  else if ((r[1]==0) & (r[2]==6))
    psi <- -15*a22^3

  ### eighth  order psi functionals - normal score 
  else if ((r[1]==8) & (r[2]==0))
    psi <- 105*a11^4
  else if ((r[1]==7) & (r[2]==1))
    psi <- 105*a11^3*a12
  else if ((r[1]==6) & (r[2]==2))
    psi <- 15*a11^2*(6*a12^2 + a11*a22)
  else if ((r[1]==5) & (r[2]==3))
    psi <- 15*a11*a12*(4*a12^2 + 3*a11*a22)
  else if ((r[1]==4) & (r[2]==4))
    psi <- 24*a12^4 + 72*a11*a12^2*a22 + 9*a11^2*a22^2
  else if ((r[1]==3) & (r[2]==5))
    psi <- 15*a12*a22*(4*a12^2 + 3*a11*a22)
  else if ((r[1]==2) & (r[2]==6))
    psi <- 15*a22^2*(6*a12^2 + a11*a22)
  else if ((r[1]==1) & (r[2]==7))
    psi <- 105*a12*a22^3
  else if ((r[1]==0) & (r[2]==8))
    psi <- 105*a22^4
  
  return(1/sqrt(2*pi)^2 * sqrt(det(A))*psi) 
}
  
psins.3d <- function(r, Sigma)
{
  A <- chol2inv(chol(2*Sigma[1:3, 1:3]))
  a11 <- A[1,1]
  a22 <- A[2,2]
  a33 <- A[3,3]
  a44 <- 0

  a12 <- A[1,2]
  a13 <- A[1,3]
  a14 <- 0
  a23 <- A[2,3]
  a24 <- 0
  a34 <- 0
  
  r <- c(r, 0) 
  if (sum(r)==6)
    psi <- psins.4d.part6(r, a11, a22, a33, a44, a12, a13, a14, a23, a24, a34)
  else if (sum(r)==8)
    psi <- psins.4d.part8(r, a11, a22, a33, a44, a12, a13, a14, a23, a24, a34)

  return(1/sqrt(2*pi)^3 * sqrt(det(A))*psi)
}

psins.4d <- function(r, Sigma)
{
  A <- chol2inv(chol(2*Sigma))
  a11 <- A[1,1]
  a22 <- A[2,2]
  a33 <- A[3,3]
  a44 <- A[4,4]

  a12 <- A[1,2]
  a13 <- A[1,3]
  a14 <- A[1,4]
  a23 <- A[2,3]
  a24 <- A[2,4]
  a34 <- A[3,4]
  
  if (sum(r)==6)
    psi <- psins.4d.part6(r, a11, a22, a33, a44, a12, a13, a14, a23, a24, a34)
  else if (sum(r)==8)
    psi <- psins.4d.part8(r, a11, a22, a33, a44, a12, a13, a14, a23, a24, a34)

  return(1/sqrt(2*pi)^4 * sqrt(det(A))*psi)
}


psins.5d <- function(r, Sigma)
{
  A <- chol2inv(chol(2*Sigma[1:5,1:5]))
  a11 <- A[1,1]
  a22 <- A[2,2]
  a33 <- A[3,3]
  a44 <- A[4,4]
  a55 <- A[3,3]
  a66 <- 0
  
  a12 <- A[1,2]
  a13 <- A[1,3]
  a14 <- A[1,4]
  a15 <- A[1,5]
  a16 <- 0
  a23 <- A[2,3]
  a24 <- A[2,4]
  a25 <- A[2,5]
  a26 <- 0
  a34 <- A[3,4]
  a35 <- A[3,5]
  a36 <- 0
  a45 <- A[4,5]
  a46 <- 0
  a56 <- 0

  r <- c(r,0)  
  if (sum(r)==6)
   psi <- psins.6d.part6(r, a11, a22, a33, a44, a55, a66, a12, a13, a14, a15,
                          a16, a23, a24, a25, a26, a34, a35, a36, a45, a46, a56)
  else if (sum(r)==8)
  {
    psi	<- psins.6d.part81(r, a11, a22, a33, a44, a55, a66, a12, a13, a14, a15,
                          a16, a23, a24, a25, a26, a34, a35, a36, a45, a46, a56)
    if (is.na(psi))
      psi <- psins.6d.part82(r, a11, a22, a33, a44, a55, a66, a12, a13, a14, a15,
                          a16, a23, a24, a25, a26, a34, a35, a36, a45, a46, a56)
    if (is.na(psi))
      psi <- psins.6d.part83(r, a11, a22, a33, a44, a55, a66, a12, a13, a14, a15,
                          a16, a23, a24, a25, a26, a34, a35, a36, a45, a46, a56)
    if (is.na(psi))
      psi <- psins.6d.part84(r, a11, a22, a33, a44, a55, a66, a12, a13, a14, a15,
                          a16, a23, a24, a25, a26, a34, a35, a36, a45, a46, a56)
  }
  return(1/sqrt(2*pi)^6 * sqrt(det(A))*psi)
}

 
psins.6d <- function(r, Sigma)
{
  A <- chol2inv(chol(2*Sigma))
  a11 <- A[1,1]
  a22 <- A[2,2]
  a33 <- A[3,3]
  a44 <- A[4,4]
  a55 <- A[3,3]
  a66 <- A[4,4]
  
  a12 <- A[1,2]
  a13 <- A[1,3]
  a14 <- A[1,4]
  a15 <- A[1,5]
  a16 <- A[1,6]
  a23 <- A[2,3]
  a24 <- A[2,4]
  a25 <- A[2,5]
  a26 <- A[2,6]
  a34 <- A[3,4]
  a35 <- A[3,5]
  a36 <- A[3,6]
  a45 <- A[4,5]
  a46 <- A[4,6]
  a56 <- A[5,6]
  
  if (sum(r)==6)
   psi <- psins.6d.part6(r, a11, a22, a33, a44, a55, a66, a12, a13, a14, a15,
                          a16, a23, a24, a25, a26, a34, a35, a36, a45, a46, a56)
  else if (sum(r)==8)
  {
    psi	<- psins.6d.part81(r, a11, a22, a33, a44, a55, a66, a12, a13, a14, a15,
                          a16, a23, a24, a25, a26, a34, a35, a36, a45, a46, a56)
    if (is.na(psi))
      psi <- psins.6d.part82(r, a11, a22, a33, a44, a55, a66, a12, a13, a14, a15,
                          a16, a23, a24, a25, a26, a34, a35, a36, a45, a46, a56)
    if (is.na(psi))
      psi <- psins.6d.part83(r, a11, a22, a33, a44, a55, a66, a12, a13, a14, a15,
                          a16, a23, a24, a25, a26, a34, a35, a36, a45, a46, a56)
    if (is.na(psi))
      psi <- psins.6d.part84(r, a11, a22, a33, a44, a55, a66, a12, a13, a14, a15,
                          a16, a23, a24, a25, a26, a34, a35, a36, a45, a46, a56)
  }
  return(1/sqrt(2*pi)^6 * sqrt(det(A))*psi)
}

 


psins.4d.part6 <- function(r, a11, a22, a33, a44, a12, a13, a14, a23, a24, a34)
{
  if ((r[1]==6) & (r[2]==0) & (r[3]==0) & (r[4]==0))
    psi <- -15*a11^3
  if ((r[1]==5) & (r[2]==1) & (r[3]==0) & (r[4]==0))
    psi <- -15*a11^2*a12
  if ((r[1]==5) & (r[2]==0) & (r[3]==1) & (r[4]==0))
    psi <- -15*a11^2*a13
  if ((r[1]==5) & (r[2]==0) & (r[3]==0) & (r[4]==1))
    psi <- -15*a11^2*a14
  if ((r[1]==4) & (r[2]==2) & (r[3]==0) & (r[4]==0))
    psi <- -3*a11*(4*a12^2 + a11*a22)
  if ((r[1]==4) & (r[2]==1) & (r[3]==1) & (r[4]==0))
    psi <- -3*a11*(4*a12*a13 + a11*a23)
  if ((r[1]==4) & (r[2]==1) & (r[3]==0) & (r[4]==1))
    psi <- -3*a11*(4*a12*a14 + a11*a24)
  if ((r[1]==4) & (r[2]==0) & (r[3]==2) & (r[4]==0))
    psi <- -3*a11*(4*a13^2 + a11*a33)
  if ((r[1]==4) & (r[2]==0) & (r[3]==1) & (r[4]==1))
    psi <- -3*a11*(4*a13*a14 + a11*a34)
  if ((r[1]==4) & (r[2]==0) & (r[3]==0) & (r[4]==2))
    psi <- -3*a11*(4*a14^2 + a11*a44)
  if ((r[1]==3) & (r[2]==3) & (r[3]==0) & (r[4]==0))
    psi <- -6*a12^3 - 9*a11*a12*a22
  if ((r[1]==3) & (r[2]==2) & (r[3]==1) & (r[4]==0))
    psi <- -3*(2*a12^2*a13 + a11*a13*a22 + 2*a11*a12*a23)
  if ((r[1]==3) & (r[2]==2) & (r[3]==0) & (r[4]==1))
    psi <- -3*(2*a12^2*a14 + a11*a14*a22 + 2*a11*a12*a24)
  if ((r[1]==3) & (r[2]==1) & (r[3]==2) & (r[4]==0))
    psi <- -3*(2*a12*a13^2 + 2*a11*a13*a23 + a11*a12*a33)
  if ((r[1]==3) & (r[2]==1) & (r[3]==1) & (r[4]==1))
    psi <- -3*(2*a12*a13*a14 + a11*a14*a23 + a11*a13*a24 + a11*a12*a34)
  if ((r[1]==3) & (r[2]==1) & (r[3]==0) & (r[4]==2))
    psi <- -3*(2*a12*a14^2 + 2*a11*a14*a24 + a11*a12*a44)
  if ((r[1]==3) & (r[2]==0) & (r[3]==3) & (r[4]==0))
    psi <- -6*a13^3 - 9*a11*a13*a33
  if ((r[1]==3) & (r[2]==0) & (r[3]==2) & (r[4]==1))
    psi <- -3*(2*a13^2*a14 + a11*a14*a33 + 2*a11*a13*a34)
  if ((r[1]==3) & (r[2]==0) & (r[3]==1) & (r[4]==2))
    psi <- -3*(2*a13*a14^2 + 2*a11*a14*a34 + a11*a13*a44)
  if ((r[1]==3) & (r[2]==0) & (r[3]==0) & (r[4]==3))
    psi <- -6*a14^3 - 9*a11*a14*a44
  if ((r[1]==2) & (r[2]==4) & (r[3]==0) & (r[4]==0))
    psi <- -3*a22*(4*a12^2 + a11*a22)
  if ((r[1]==2) & (r[2]==3) & (r[3]==1) & (r[4]==0))
    psi <- -3*(2*a12*a13*a22 + 2*a12^2*a23 + a11*a22*a23)
  if ((r[1]==2) & (r[2]==3) & (r[3]==0) & (r[4]==1))
    psi <- -3*(2*a12*a14*a22 + 2*a12^2*a24 + a11*a22*a24)
  if ((r[1]==2) & (r[2]==2) & (r[3]==2) & (r[4]==0))
    psi <- -2*a13^2*a22 - 8*a12*a13*a23 - 2*a11*a23^2 - 2*a12^2*a33 - a11*a22*a33
  if ((r[1]==2) & (r[2]==2) & (r[3]==1) & (r[4]==1))
    psi <- -2*a13*a14*a22 - 4*a12*a14*a23 - 4*a12*a13*a24 - 2*a11*a23*a24 - 2*a12^2*a34 - a11*a22*a34
  if ((r[1]==2) & (r[2]==2) & (r[3]==0) & (r[4]==2))
    psi <- -2*a14^2*a22 - 8*a12*a14*a24 - 2*a11*a24^2 - 2*a12^2*a44 - a11*a22*a44
  if ((r[1]==2) & (r[2]==1) & (r[3]==3) & (r[4]==0))
    psi <- -3*(2*a13^2*a23 + 2*a12*a13*a33 + a11*a23*a33)
  if ((r[1]==2) & (r[2]==1) & (r[3]==2) & (r[4]==1))
    psi <- -4*a13*a14*a23 - 2*a13^2*a24 - 2*a12*a14*a33 - a11*a24*a33 - 4*a12*a13*a34 - 2*a11*a23*a34
  if ((r[1]==2) & (r[2]==1) & (r[3]==1) & (r[4]==2))
    psi <- -2*a14^2*a23 - 4*a13*a14*a24 - 4*a12*a14*a34 - 2*a11*a24*a34 - 2*a12*a13*a44 - a11*a23*a44
  if ((r[1]==2) & (r[2]==1) & (r[3]==0) & (r[4]==3))
    psi <- -3*(2*a14^2*a24 + 2*a12*a14*a44 + a11*a24*a44)
  if ((r[1]==2) & (r[2]==0) & (r[3]==4) & (r[4]==0))
    psi <- -3*a33*(4*a13^2 + a11*a33)
  if ((r[1]==2) & (r[2]==0) & (r[3]==3) & (r[4]==1))
    psi <- -3*(2*a13*a14*a33 + 2*a13^2*a34 + a11*a33*a34)
  if ((r[1]==2) & (r[2]==0) & (r[3]==2) & (r[4]==2))
    psi <- -2*a14^2*a33 - 8*a13*a14*a34 - 2*a11*a34^2 - 2*a13^2*a44 - a11*a33*a44
  if ((r[1]==2) & (r[2]==0) & (r[3]==1) & (r[4]==3))
    psi <- -3*(2*a14^2*a34 + 2*a13*a14*a44 + a11*a34*a44)
  if ((r[1]==2) & (r[2]==0) & (r[3]==0) & (r[4]==4))
    psi <- -3*a44*(4*a14^2 + a11*a44)
  if ((r[1]==1) & (r[2]==5) & (r[3]==0) & (r[4]==0))
    psi <- -15*a12*a22^2
  if ((r[1]==1) & (r[2]==4) & (r[3]==1) & (r[4]==0))
    psi <- -3*a22*(a13*a22 + 4*a12*a23)
  if ((r[1]==1) & (r[2]==4) & (r[3]==0) & (r[4]==1))
    psi <- -3*a22*(a14*a22 + 4*a12*a24)
  if ((r[1]==1) & (r[2]==3) & (r[3]==2) & (r[4]==0))
    psi <- -3*(2*a13*a22*a23 + 2*a12*a23^2 + a12*a22*a33)
  if ((r[1]==1) & (r[2]==3) & (r[3]==1) & (r[4]==1))
    psi <- -3*(a14*a22*a23 + a13*a22*a24 + 2*a12*a23*a24 + a12*a22*a34)
  if ((r[1]==1) & (r[2]==3) & (r[3]==0) & (r[4]==2))
    psi <- -3*(2*a14*a22*a24 + 2*a12*a24^2 + a12*a22*a44)
  if ((r[1]==1) & (r[2]==2) & (r[3]==3) & (r[4]==0))