Golang DenseMatrix.Cols方法代码示例

func MatrixNormal(M, Omega, Sigma *mx.DenseMatrix) func() (X *mx.DenseMatrix) {
	checkMatrixNormal(M, Omega, Sigma)

	Mv := mx.Vectorize(M)
	Cov := mx.Kronecker(Omega, Sigma)
	normal := MVNormal(Mv, Cov)
	return func() (X *mx.DenseMatrix) {
		Xv := normal()
		X = mx.Unvectorize(Xv, M.Rows(), M.Cols())
		return
	}
}

func MatrixNormal_LnPDF(M, Omega, Sigma *mx.DenseMatrix) func(A *mx.DenseMatrix) float64 {
	checkMatrixNormal(M, Omega, Sigma)

	pf := float64(M.Rows())
	mf := float64(M.Cols())

	sinv, err := Sigma.Inverse()
	if err != nil {
		panic(err)
	}
	oinv, err := Omega.Inverse()
	if err != nil {
		panic(err)
	}

	norm := (2 * math.Pi) * (-0.5 * mf * pf)
	norm += Omega.Det() * (-0.5 * mf)
	norm += Sigma.Det() * (-0.5 * pf)

	return func(X *mx.DenseMatrix) (lp float64) {
		lp = norm
		diff, err := X.MinusDense(M)
		if err != nil {
			panic(err)
		}
		inner := oinv

		inner, err = inner.TimesDense(diff.Transpose())
		if err != nil {
			panic(err)
		}

		inner, err = inner.TimesDense(sinv)
		if err != nil {
			panic(err)
		}

		inner, err = inner.TimesDense(diff)
		if err != nil {
			panic(err)
		}

		innerTrace := inner.Trace()

		lp += -0.5 * innerTrace

		return
	}
}

/*
 M is the mean, Omega is the row covariance, Sigma is the column covariance.
*/
func MatrixNormal_PDF(M, Omega, Sigma *mx.DenseMatrix) func(A *mx.DenseMatrix) float64 {
	checkMatrixNormal(M, Omega, Sigma)
	pf := float64(M.Rows())
	mf := float64(M.Cols())

	norm := math.Pow(2*math.Pi, -0.5*mf*pf)
	norm *= math.Pow(Omega.Det(), -0.5*mf)
	norm *= math.Pow(Sigma.Det(), -0.5*pf)

	return func(X *mx.DenseMatrix) (p float64) {
		p = norm

		sinv, err := Sigma.Inverse()
		if err != nil {
			panic(err)
		}
		oinv, err := Omega.Inverse()
		if err != nil {
			panic(err)
		}
		diff, err := X.MinusDense(M)
		if err != nil {
			panic(err)
		}
		inner := oinv

		inner, err = inner.TimesDense(diff.Transpose())
		if err != nil {
			panic(err)
		}

		inner, err = inner.TimesDense(sinv)
		if err != nil {
			panic(err)
		}

		inner, err = inner.TimesDense(diff)
		if err != nil {
			panic(err)
		}

		innerTrace := inner.Trace()

		p *= math.Exp(-0.5 * innerTrace)

		return
	}
}

func MatrixT(M, Omega, Sigma *mx.DenseMatrix, n int) func() (T *mx.DenseMatrix) {
	checkMatrixT(M, Omega, Sigma, n)

	fmt.Println("M:", M)
	fmt.Println("Sigma:", Sigma)
	fmt.Println("Omega:", Omega)

	p := M.Rows()
	m := M.Cols()

	OmegaInv, err := Omega.Inverse()
	if err != nil {
		panic(err)
	}

	Sdist := Wishart(n+p-1, OmegaInv)

	Xdist := MatrixNormal(mx.Zeros(p, m), mx.Eye(p), Sigma)

	return func() (T *mx.DenseMatrix) {
		S := Sdist()
		Sinv, err := S.Inverse()
		if err != nil {
			panic(err)
		}
		Sinvc, err := Sinv.Cholesky()
		if err != nil {
			panic(err)
		}
		X := Xdist()
		fmt.Println("Sinvc:", Sinvc)
		fmt.Println("X:", X)
		T, err = Sinvc.Transpose().TimesDense(X)
		if err != nil {
			panic(err)
		}
		err = T.AddDense(M)
		if err != nil {
			panic(err)
		}
		return
	}
}

func MatrixTLnPDF(M, Omega, Sigma *mx.DenseMatrix, n int) func(T *mx.DenseMatrix) (ll float64) {
	checkMatrixT(M, Omega, Sigma, n)

	nf := float64(n)
	p := M.Rows()
	pf := float64(p)
	m := M.Cols()
	mf := float64(m)

	var norm float64 = 0

	norm += logΓpr(p, 0.5*(nf+mf+pf-1), 0.5*(nf+pf-1))
	norm += pow(π, -0.5*mf*pf)
	norm += pow(Omega.Det(), -0.5*mf)
	norm += pow(Sigma.Det(), -0.5*pf)

	SigmaInv, err := Sigma.Inverse()
	if err != nil {
		panic(err)
	}
	OmegaInv, err := Omega.Inverse()
	if err != nil {
		panic(err)
	}

	return func(T *mx.DenseMatrix) (ll float64) {
		ll = norm

		diff, err := T.MinusDense(M)
		if err != nil {
			panic(err)
		}
		inner := OmegaInv.Copy()
		inner, _ = inner.TimesDense(diff)
		inner, _ = inner.TimesDense(SigmaInv)
		inner, _ = inner.TimesDense(diff.Transpose())

		ll += log(inner.Det()) * -0.5 * (nf + mf + pf - 1)

		return
	}
}

func MatrixT_PDF(M, Omega, Sigma *mx.DenseMatrix, n int) func(T *mx.DenseMatrix) (l float64) {
	checkMatrixT(M, Omega, Sigma, n)

	nf := float64(n)
	p := M.Rows()
	pf := float64(p)
	m := M.Cols()
	mf := float64(m)

	var norm float64 = 1

	norm *= GammaPRatio(p, 0.5*(nf+mf+pf-1), 0.5*(nf+pf-1))
	norm *= math.Pow(math.Pi, -0.5*mf*pf)
	norm *= math.Pow(Omega.Det(), -0.5*mf)
	norm *= math.Pow(Sigma.Det(), -0.5*pf)

	SigmaInv, err := Sigma.Inverse()
	if err != nil {
		panic(err)
	}
	OmegaInv, err := Omega.Inverse()
	if err != nil {
		panic(err)
	}

	return func(T *mx.DenseMatrix) (l float64) {
		l = norm

		diff, err := T.MinusDense(M)
		if err != nil {
			panic(err)
		}
		inner := OmegaInv.Copy()
		inner, _ = inner.TimesDense(diff)
		inner, _ = inner.TimesDense(SigmaInv)
		inner, _ = inner.TimesDense(diff.Transpose())

		l *= math.Pow(inner.Det(), -0.5*(nf+mf+pf-1))

		return
	}
}

func checkMatrixNormal(M, Omega, Sigma *mx.DenseMatrix) {
	p := M.Rows()
	m := M.Cols()
	if Omega.Rows() != p {
		panic(fmt.Sprintf("Omega.Rows != M.Rows, %d != %d", Omega.Rows(), p))
	}
	if Omega.Cols() != Omega.Rows() {
		panic("Omega is not square")
	}
	if Sigma.Rows() != m {
		panic(fmt.Sprintf("Sigma.Cols != M.Cols, %d != %d", Sigma.Cols(), m))
	}
	if Sigma.Cols() != Sigma.Rows() {
		panic("Sigma is not square")
	}
}

/*
 M is r x c, o x i
 Sigma is r x r, o x o
 Phi is c x c, i x i

 Sigma matches Y o x 1 output dimension
 Phi matches X i x 1 input dimension
*/
func NewKnownVarianceLRPosterior(M, Sigma, Phi *mx.DenseMatrix) (this *KnownVarianceLRPosterior) {
	if M.Rows() != Sigma.Rows() {
		panic("M.Rows != Sigma.Rows")
	}
	if M.Cols() != Phi.Cols() {
		panic("M.Cols != Phi.Cols")
	}
	if Sigma.Rows() != Sigma.Cols() {
		panic("Sigma is not square")
	}
	if Phi.Rows() != Phi.Cols() {
		panic("Phi is not square")
	}
	this = &KnownVarianceLRPosterior{
		M:     M,
		Sigma: Sigma,
		Phi:   Phi,
		XXt:   mx.Zeros(Phi.Cols(), Phi.Cols()),
		YXt:   mx.Zeros(Sigma.Cols(), Phi.Cols()),
	}

	return
}

/*
	If Y ~ N(AX, Sigma, I)
	and A ~ N(M, Sigma, Phi)
	this returns a sampler for P(A|X,Y,Sigma,M,Phi)
*/
func KnownVariancePosterior(Y, X, Sigma, M, Phi *mx.DenseMatrix) func() (A *mx.DenseMatrix) {
	o := Y.Rows()
	i := X.Rows()
	n := Y.Cols()
	if n != X.Cols() {
		panic("X and Y don't have the same number of columns")
	}
	if o != M.Rows() {
		panic("Y.Rows != M.Rows")
	}
	if i != M.Cols() {
		panic("Y.Rows != M.Cols")
	}
	if o != Sigma.Rows() {
		panic("Y.Rows != Sigma.Rows")
	}
	if Sigma.Cols() != Sigma.Rows() {
		panic("Sigma is not square")
	}
	if i != Phi.Rows() {
		panic("X.Rows != Phi.Rows")
	}
	if Phi.Cols() != Phi.Rows() {
		panic("Phi is not square")
	}

	Xt := X.Transpose()

	PhiInv, err := Phi.Inverse()
	if err != nil {
		panic(err)
	}

	XXt, err := X.TimesDense(Xt)
	if err != nil {
		panic(err)
	}

	XXtpPhiInv, err := XXt.PlusDense(PhiInv)
	if err != nil {
		panic(err)
	}

	Omega, err := XXtpPhiInv.Inverse()
	if err != nil {
		panic(err)
	}

	YXtpMPhiInv, err := Y.TimesDense(Xt)
	if err != nil {
		panic(err)
	}

	MPhiInv, err := M.TimesDense(PhiInv)
	if err != nil {
		panic(err)
	}

	err = YXtpMPhiInv.AddDense(MPhiInv)
	if err != nil {
		panic(err)
	}

	Mxy, err := YXtpMPhiInv.TimesDense(Omega)
	if err != nil {
		panic(err)
	}

	if false {
		fmt.Printf("Mxy:\n%v\n", Mxy)
		fmt.Printf("Sigma:\n%v\n", Sigma)
		fmt.Printf("Omega:\n%v\n", Omega)
	}

	return dst.MatrixNormal(Mxy, Sigma, Omega)
}