func estep(alpha []float32, depths []float32, lambda float32, aik [][]float32) [][]float32 {
	N := len(depths)
	if len(aik) == 0 {
		aik = make([][]float32, len(alpha))
		for i := range aik {
			aik[i] = make([]float32, N)
		}
	}
	denom := make([]float32, len(depths))
	n := len(alpha)
	for k, d := range depths {
		for i := 0; i < n; i++ {
			denom[k] += alpha[i] * float32(pmf(int(d+0.5), float64(float32(i)/2*lambda)))
		}
	}
	for cn := range alpha {
		pdepth(depths, cn, lambda, aik[cn])
		for k, ad := range aik[cn] {
			aik[cn][k] = alpha[cn] * ad / denom[k]
		}
	}
	return aik
}
returns a square alpha ik.
package mops
import "math"
func mean32(a []float32) float32 {
	var sum float32
	for _, v := range a {
		sum += v
	}
	return sum / float32(len(a))
}
func abs(a float64) float64 {
	if a < 0 {
		return -a
	}
	return a
}
func abs32(a float32) float32 {
	if a < 0 {
		return -a
	}
	return a
}
const maxCN = 8
const delta = 0.001
const maxiter = 10
const eps = 0.001
func pmf(k int, pmean float64) float64 {
	return math.Pow(pmean, float64(k)) * math.Exp(-pmean) / math.Gamma(float64(k+1))
}
func pdepth(depths []float32, cn int, lambda float32, out []float32) {
	beta := float64(cn) / 2 * float64(lambda)
	if cn == 0 {
		beta = 0.001 / 2 * float64(lambda)
	}
	for k, d := range depths {
		out[k] = float32(pmf(int(d+0.5), beta))
	}
}