Initial commit for distribution package. Contains Exponential, Laplace, ...

...Normal, and Uniform
Your idea of using the values of the distribution has merit. Do you know if that works for all distributions in the exponential family (i.e. the sufficient statistics can always be related directly to the parameters of the distribution)? Most statistics packages do not represent the sufficient statistics in this way (see the wikipedia article on conjugate distributions for an example), so it's hard to use other implementations to check. I chose to represent them in this way because I personally think it makes way more sense to say one saw 6 values with mean 2 and standard deviation 12 than it does to say I saw 6 values with total sum 12 and total sum of squares 150. 

If we used your method, how would I abstract the idea of fitting a series of datasets with the same prior?

type FitVarer interface {
   Variance() float64
   FitPrior(x,w,p, weight []float64) (np, nw []float64)
}

func FitVariances(data [][]float64, distribution FitPriorable, priorValues, priorWeights []float64) (variances []float64){
    variances := make([]float64, len(data)) 
    for i, set := range data {
         distribution.FitPrior(set, nil, priorValues, priorWeights)
         variances[i] = distribution.Variance()
     }
}

If the distribution is changed in place, then you need some way to "reset" the prior on each pass of the for loop. The only way I can see to do that is if there's some abstract way to set the parameters of the distribution (like UnmarshalSlice...)