using LinearAlgebra, StatsBase, Optim
function sfran_loglike(mu,sigma2_mu,sigma2,y)
   d, ni = size(y)
   Sigmai = sigma2*diagm(ones(ni)) .+ sigma2_mu*ones(ni,ni)
   l = -(ni*d)/2*log(2*pi) - d/2*log(det(Sigmai))
   for i=1:d
     yi = y[i,:];
     l = l - .5*(yi .- mu)'*(Sigmai\(yi .- mu))
   end
   return l
end

y = [ 22.6 20.5 20.8; 22.6 21.2 20.5; 17.3 16.2 16.6; 
      21.4 23.7 23.2; 20.9 22.2 22.6; 14.5 10.5 12.3; 
      20.8 19.1 21.3; 17.4 18.6 18.6; 25.1 24.8 24.9; 
      14.9 16.3 16.6];
f =  theta -> -sfran_loglike(theta[1],theta[2],theta[3],y)
ybar = mean(y,dims=2)
theta0 = [mean(ybar) var(ybar) mean(var(y,dims=2))]
res = optimize(f,theta0);
thetahat = res.minimizer