R/ppi_W11.R

In scorematchingad: Score Matching Estimation by Automatic Differentiation

#' @noRd
#' @param p The number of dimensions
#' @param z The square root of the data values
#' @param h A single column matrix of weights for each row of data
#' @param ind A matrix of index combinations. This can be generated by `indexcombinations(p - 1)`
#' @param qind The number of different index combinations (i.e. columns of `ind`).
#'  (i.e. the weight function evaluated at each measurement)
calcW11 <- function(p, z, h, ind, qind, w = rep(1, nrow(z))){
	sp=p-1

	A1=matrix(0,1,p-1)
	for (j in 1:sp)
	{
		A1[j]=weighted.mean(16*(h^2)*z[,j]^6, w = w)
	}

	x=c(1:sp)


	C1=matrix(0,1,p-1)
	for (j in 1:sp)
	{
		C1[j]=weighted.mean(8*(h^2)*z[,j]^4, w = w)
	}

	AA=diag(as.vector(A1))


	AB=matrix(0,sp,qind)
	for (i in 1:sp)
	{
		for (j in 1:qind)
		{
			if (i==ind[1,j]){AB[i,j]=weighted.mean(16*z[,ind[1,j]]^4*z[,ind[2,j]]^2*h^2, w=w)}
			if (i==ind[2,j]){AB[i,j]=weighted.mean(16*z[,ind[1,j]]^2*z[,ind[2,j]]^4*h^2, w=w)}
		}
	}

	BB=matrix(0,qind,qind)
	for (i in 1:qind)
	{
		for (j in 1:qind)
		{
			if (ind[1,i]==ind[1,j] & ind[2,i]==ind[2,j]){BB[i,j]=weighted.mean(16*z[,ind[1,j]]^4*z[,ind[2,j]]^2*h^2, w=w)+weighted.mean(16*z[,ind[1,j]]^2*z[,ind[2,j]]^4*h^2, w=w)}
			else if (ind[1,i]==ind[1,j]){BB[i,j]=weighted.mean((h^2)*16*z[,ind[2,j]]^2*z[,ind[2,i]]^2*z[,ind[1,j]]^2, w=w)}
			else if (ind[2,i]==ind[2,j]){BB[i,j]=weighted.mean((h^2)*16*z[,ind[1,j]]^2*z[,ind[1,i]]^2*z[,ind[2,j]]^2, w=w)}
			else if (ind[1,i]==ind[2,j]){BB[i,j]=weighted.mean((h^2)*16*z[,ind[2,i]]^2*z[,ind[1,j]]^2*z[,ind[2,j]]^2, w=w)}
			else if (ind[2,i]==ind[1,j]){BB[i,j]=weighted.mean((h^2)*16*z[,ind[1,i]]^2*z[,ind[2,j]]^2*z[,ind[1,j]]^2, w=w)}

		}
	}

	AC=diag(as.vector(C1))



	BC=matrix(0,qind,sp)
	for (i in 1:qind)
	{
		for (j in 1:sp)
		{
			if (j==ind[1,i]){BC[i,j]=weighted.mean(h^2*z[,ind[2,i]]^2*z[,ind[1,i]]^2*8, w=w)}
			else if (j==ind[2,i]){BC[i,j]=weighted.mean(h^2*z[,ind[1,i]]^2*z[,ind[2,i]]^2*8, w=w)}
		}
	}

	C2=matrix(0,1,p-1)
	for (j in 1:sp)
	{
		C2[j]=weighted.mean(4*(h^2)*z[,j]^2, w = w)
	}

	CC=diag(as.vector(C2))



	AA2=AA
	for (i in 1:sp)
	{
		for (j in 1:sp)
		{
			AA2[i,j]=16*weighted.mean((h^2)*z[,i]^4*z[,j]^4, w=w)
		}
	}



	AB2=AB
	for (i in 1:sp)
	{
		for (j in 1:qind)
		{
			AB2[i,j]=32*weighted.mean((h^2)*z[,i]^4*z[,ind[1,j]]^2*z[,ind[2,j]]^2, w=w)
		}
	}



	BB2=BB
	for (i in 1:qind)
	{
		for (j in 1:qind)
		{
			BB2[i,j]=64*weighted.mean((h^2)*(z[,ind[1,i]]^2*z[,ind[2,i]]^2)*(z[,ind[1,j]]^2*z[,ind[2,j]]^2), w=w)
		}
	}


	BC2=BC
	for (i in 1:qind)
	{
		for (j in 1:sp)
		{
			BC2[i,j]=16*weighted.mean((h^2)*z[,j]^2*(z[,ind[1,i]]^2*z[,ind[2,i]]^2), w=w)
		}
	}



	AC2=AC
	for (i in 1:sp)
	{
		for (j in 1:sp)
		{
			AC2[i,j]=8*weighted.mean((h^2)*z[,i]^4*z[,j]^2, w=w)
		}
	}

	CC2=AC
	for (i in 1:sp)
	{
		for (j in 1:sp)
		{
			CC2[i,j]=4*weighted.mean((h^2)*z[,i]^2*z[,j]^2, w=w)
		}
	}



	W1=cbind(AA,AB,AC)
	W2=cbind(t(AB),BB,BC)
	W3=cbind(t(AC),t(BC),CC)
	W0=rbind(W1,W2,W3)


	W12=cbind(AA2,AB2,AC2)
	W22=cbind(t(AB2),BB2,BC2)
	W32=cbind(t(AC2),t(BC2),CC2)
	W2=rbind(W12,W22,W32)

	W=W0-W2
	return(W)
}