C#，数值计算——分类与推理Gaumixmod的计算方法与源程序

1 文本格式

using System;
using System.Collections.Generic;

namespace Legalsoft.Truffer
{
    public class Gaumixmod
    {
        private int nn { get; set; }
        private int kk { get; set; }
        private int mm { get; set; }
        private double[,] data { get; set; }
        private double[,] means { get; set; }
        private double[,] resp { get; set; }
        private double[] frac { get; set; }
        private double[] lndets { get; set; }
        private double[,,] sig { get; set; }
        private double loglike { get; set; }

        public Gaumixmod(double[,] ddata, double[,] mmeans)
        {
            int mmstat = ddata.GetLength(1);

            this.nn = ddata.GetLength(0);
            this.kk = mmeans.GetLength(0);
            this.mm = mmstat;
            this.data = Globals.CopyFrom(ddata);
            this.means = Globals.CopyFrom(mmeans);
            this.resp = new double[nn, kk];
            this.frac = new double[kk];
            this.lndets = new double[kk];
            this.sig = new double[kk, mmstat, mmstat];

            for (int k = 0; k < kk; k++)
            {
                frac[k] = 1.0 / kk;
                for (int i = 0; i < mm; i++)
                {
                    for (int j = 0; j < mm; j++)
                    {
                        sig[k, i, j] = 0.0;
                    }
                    sig[k, i, i] = 1.0e-10;
                }
            }

            estep();
            mstep();
        }

        public double estep()
        {
            double[] u = new double[mm];
            double[] v = new double[mm];
            double oldloglike = loglike;
            for (int k = 0; k < kk; k++)
            {
                //Cholesky choltmp = new Cholesky(sig[k]);
                Cholesky choltmp = new Cholesky(Globals.CopyFrom(k, sig));
                lndets[k] = choltmp.logdet();
                for (int n = 0; n < nn; n++)
                {
                    for (int m = 0; m < mm; m++)
                    {
                        u[m] = data[n, m] - means[k, m];
                    }
                    choltmp.elsolve(u, v);
                    double sum = 0.0;
                    for (int m = 0; m < mm; m++)
                    {
                        sum += Globals.SQR(v[m]);
                    }
                    resp[n, k] = -0.5 * (sum + lndets[k]) + Math.Log(frac[k]);
                }
            }
            loglike = 0;
            for (int n = 0; n < nn; n++)
            {
                double max = -99.9e99;
                for (int k = 0; k < kk; k++)
                {
                    if (resp[n, k] > max)
                    {
                        max = resp[n, k];
                    }
                }
                double sum = 0.0;
                for (int k = 0; k < kk; k++)
                {
                    sum += Math.Exp(resp[n, k] - max);
                }
                double tmp = max + Math.Log(sum);
                for (int k = 0; k < kk; k++)
                {
                    resp[n, k] = Math.Exp(resp[n, k] - tmp);
                }
                loglike += tmp;
            }
            return loglike - oldloglike;
        }

        public void mstep()
        {
            for (int k = 0; k < kk; k++)
            {
                double wgt = 0.0;
                for (int n = 0; n < nn; n++)
                {
                    wgt += resp[n, k];
                }
                frac[k] = wgt / nn;
                for (int m = 0; m < mm; m++)
                {
                    double sum = 0.0;
                    for (int n = 0; n < nn; n++)
                    {
                        sum += resp[n, k] * data[n, m];
                    }
                    means[k, m] = sum / wgt;
                    for (int j = 0; j < mm; j++)
                    {
                        sum = 0.0;
                        for (int n = 0; n < nn; n++)
                        {
                            sum += resp[n, k] * (data[n, m] - means[k, m]) * (data[n, j] - means[k, j]);
                        }
                        sig[k, m, j] = sum / wgt;
                    }
                }
            }
        }
    }
}
 

2 代码格式

using System;
using System.Collections.Generic;namespace Legalsoft.Truffer
{public class Gaumixmod{private int nn { get; set; }private int kk { get; set; }private int mm { get; set; }private double[,] data { get; set; }private double[,] means { get; set; }private double[,] resp { get; set; }private double[] frac { get; set; }private double[] lndets { get; set; }private double[,,] sig { get; set; }private double loglike { get; set; }public Gaumixmod(double[,] ddata, double[,] mmeans){int mmstat = ddata.GetLength(1);this.nn = ddata.GetLength(0);this.kk = mmeans.GetLength(0);this.mm = mmstat;this.data = Globals.CopyFrom(ddata);this.means = Globals.CopyFrom(mmeans);this.resp = new double[nn, kk];this.frac = new double[kk];this.lndets = new double[kk];this.sig = new double[kk, mmstat, mmstat];for (int k = 0; k < kk; k++){frac[k] = 1.0 / kk;for (int i = 0; i < mm; i++){for (int j = 0; j < mm; j++){sig[k, i, j] = 0.0;}sig[k, i, i] = 1.0e-10;}}estep();mstep();}public double estep(){double[] u = new double[mm];double[] v = new double[mm];double oldloglike = loglike;for (int k = 0; k < kk; k++){//Cholesky choltmp = new Cholesky(sig[k]);Cholesky choltmp = new Cholesky(Globals.CopyFrom(k, sig));lndets[k] = choltmp.logdet();for (int n = 0; n < nn; n++){for (int m = 0; m < mm; m++){u[m] = data[n, m] - means[k, m];}choltmp.elsolve(u, v);double sum = 0.0;for (int m = 0; m < mm; m++){sum += Globals.SQR(v[m]);}resp[n, k] = -0.5 * (sum + lndets[k]) + Math.Log(frac[k]);}}loglike = 0;for (int n = 0; n < nn; n++){double max = -99.9e99;for (int k = 0; k < kk; k++){if (resp[n, k] > max){max = resp[n, k];}}double sum = 0.0;for (int k = 0; k < kk; k++){sum += Math.Exp(resp[n, k] - max);}double tmp = max + Math.Log(sum);for (int k = 0; k < kk; k++){resp[n, k] = Math.Exp(resp[n, k] - tmp);}loglike += tmp;}return loglike - oldloglike;}public void mstep(){for (int k = 0; k < kk; k++){double wgt = 0.0;for (int n = 0; n < nn; n++){wgt += resp[n, k];}frac[k] = wgt / nn;for (int m = 0; m < mm; m++){double sum = 0.0;for (int n = 0; n < nn; n++){sum += resp[n, k] * data[n, m];}means[k, m] = sum / wgt;for (int j = 0; j < mm; j++){sum = 0.0;for (int n = 0; n < nn; n++){sum += resp[n, k] * (data[n, m] - means[k, m]) * (data[n, j] - means[k, j]);}sig[k, m, j] = sum / wgt;}}}}}
}