p2: add a way to estimate impulse function

p2/stats/helper.go

@@ -0,0 +1,82 @@
+package stats
+
+import (
+	"errors"
+	"log"
+
+	"gonum.org/v1/gonum/mat"
+)
+
+func rotateFloats(f []float64, n int) []float64 {
+	flen := len(f)
+	if n < 0 || flen == 0 {
+		return f
+	}
+
+	r := flen - n%flen
+	f = append(f[r:], f[:r]...)
+
+	return f
+}
+
+func constructMatrix(f []float64) ([][]float64, error) {
+	flen := len(f)
+
+	if flen < 1 {
+		return nil, errors.New("ErrSliceZeroLength")
+	}
+
+	m := make([][]float64, 0, flen)
+
+	// add first row unchanged.
+	m = append(m, f)
+
+	rot := rotateFloats(f, 1)
+
+	for i := 0; i < flen-1; i++ {
+		m = append(m, rot)
+		rot = rotateFloats(rot, 1)
+	}
+
+	return m, nil
+}
+
+func invertMatrix(m [][]float64, size int) (*mat.Dense, error) {
+	flM := []float64{}
+	for i := range m {
+		flM = append(flM, m[i]...)
+	}
+
+	a := mat.NewDense(size, size, flM)
+
+	// compute the inverse of A, if possible.
+	var aInv mat.Dense
+
+	err := aInv.Inverse(a)
+	if err != nil {
+		log.Print("matrix is not invertible")
+		return nil, errors.New("matrix is not invertible")
+	}
+
+	return &aInv, nil
+}
+
+func estimateImpulseFunc(gMat mat.Dense, lruu int) ([]float64, error) {
+	g := gMat.RawMatrix().Data
+	lenG := len(g)
+
+	if lenG != lruu {
+		return nil, errors.New("ErrLenGDiffersFromLRuu")
+	}
+
+	h := make([]float64, lenG)
+
+	// set first elem manually.
+	h[0] = g[0]
+
+	for i := 1; i < lenG; i++ {
+		h[i] = h[i-1] + g[i]
+	}
+
+	return h, nil
+}

p2/stats/stats.go

@@ -4,6 +4,7 @@ import (
 	"errors"
 	"math"
 
+	"gonum.org/v1/gonum/mat"
 	"gonum.org/v1/gonum/stat"
 )
 
@@ -93,3 +94,37 @@ func MutCorrelate(f1, f2 []float64, maxShift float64) ([]float64, error) {
 func Covariance(f1, f2 []float64) float64 {
 	return stat.Covariance(f1, f2, nil)
 }
+
+func ImpulseFunction(ruu []float64, ruy []float64) ([]float64, error) {
+	lruu := len(ruu)
+
+	// construct matrix out of ruu.
+	ruuMatrix, err := constructMatrix(ruu)
+	if err != nil {
+		return nil, err
+	}
+
+	// construct inverse matrix out of ruu matrix.
+	invRuuMatrix, err := invertMatrix(ruuMatrix, lruu)
+	if err != nil {
+		return nil, err
+	}
+
+	// matrixify ruy.
+	// ruyVect := mat.NewDense(lruu, 1, ruy)
+	ruyVect := mat.NewDense(1, lruu, ruy)
+
+	// get g - the dot product of ruy vector and ruu inverted matrix.
+	var g mat.Dense
+
+	// g.Mul(invRuuMatrix, ruyVect)
+	g.Mul(ruyVect, invRuuMatrix)
+
+	// calculate h - the impulse function.
+	h, err := estimateImpulseFunc(g, lruu)
+	if err != nil {
+		return nil, err
+	}
+
+	return h, nil
+}