EM算法求解高斯混合模型

高斯混合模型用EM求解

public class GausssianEM {
    public static void main(String[] args) {
        int k=4,n=10,i=1;
        Gaussian[] gaussians=new Gaussian[k];
        double[][] R=new double[n][k];
        gaussians[0]=new Gaussian(0.1,2,1);
        gaussians[1]=new Gaussian(0.4,0,0.4);
        gaussians[2]=new Gaussian(0.3,3,2);
        gaussians[3]=new Gaussian(0.2,-1,6);
        Gaussian[] newGs=gaussians; 
        double[] Y={1,-0.3,4.5,2,1.2,8,10,-2.7,3.1,5};
        while(i<=1000){
            R=Estep(Y,gaussians);
            newGs=Mstep(R,Y,gaussians);
            if(distance(newGs, gaussians)<0.0000001){
                gaussians=newGs;break;
            }
            gaussians=newGs;i++;
        }
        System.out.println("迭代次数:"+i);
        for(Gaussian gs:gaussians)
            System.out.println(gs.toString());
    }
    public static double[][] Estep(double[] Y,Gaussian[] gaussians){
        double[][] R=new double[Y.length][gaussians.length];
        double sum;
        for(int j=0;j<Y.length;j++){
            sum=0;
            for(int k=0;k<gaussians.length;k++){
                R[j][k]=gaussians[k].valueAt(Y[j]);
                sum+=R[j][k];
            }
            for(int k=0;k<gaussians.length;k++) R[j][k]=R[j][k]/sum;
        }
        return R;
    }
    public static Gaussian[] Mstep(double[][] R,double[] Y,Gaussian[] gaussians){
        Gaussian[] newGs=new Gaussian[gaussians.length];
        double sumRk,sumRk_Y,sumRk_YU;
        for(int k=0;k<gaussians.length;k++){
            sumRk=0;sumRk_Y=0;sumRk_YU=0;
            for(int j=0;j<Y.length;j++){
                sumRk+=R[j][k];
                sumRk_Y+=R[j][k]*Y[j];
                sumRk_YU+=R[j][k]*Math.pow(Y[j]-gaussians[k].avgValue, 2);
            }
            newGs[k]=new Gaussian(sumRk/Y.length, sumRk_Y/sumRk, sumRk_YU/sumRk);
        }           
        return newGs;
    }
    public static double distance(Gaussian[] news,Gaussian[] olds){
        double distance=0;
        for(int i=0;i<olds.length;i++){
            distance+=Math.pow(news[i].weight-olds[i].weight, 2);
            distance+=Math.pow(news[i].var-olds[i].var, 2);
            distance+=Math.pow(news[i].avgValue-olds[i].avgValue, 2);
        }
        return Math.pow(distance, 0.5);
    }
}
class Gaussian{
    double weight,var,avgValue;
    public Gaussian(double weight, double avgValue,double var) {
        this.weight = weight;this.var = var;this.avgValue = avgValue;
    }
    public double valueAt(double y){
        double result=Math.exp(-(double)Math.pow((y-avgValue),2)/(2*var));
        result=(1.0/(Math.pow(6.28*var,0.5)))*result;
        return result;
    }
    public String toString() {
        return "Gaussian [weight=" + weight + ", var=" + var + ", avgValue="
                + avgValue + "]";
    }
}