C#，数值计算——Multinormaldev的计算方法与源程序

 

1 文本格式

using System;

namespace Legalsoft.Truffer
{
    public class Multinormaldev : Ran
    {
        public Cholesky chol { get; set; } = null;

        private int mm { get; set; }
        private double[] mean { get; set; }
        private double[,] xvar { get; set; }
        private double[] spt { get; set; }
        private double[] pt { get; set; }

        public Multinormaldev(ulong j, double[] mmean, double[,] vvar) : base(j)
        {
            this.mm = mmean.Length;
            this.mean = Globals.CopyFrom(mmean);
            this.xvar = Globals.CopyFrom(vvar);
            this.chol = new Cholesky(xvar);
            this.spt = new double[mm];
            this.pt = new double[mm];

            if (xvar.GetLength(1) != mm || xvar.GetLength(0) != mm)
            {
                throw new Exception("bad sizes");
            }
        }

        public double[] dev()
        {
            double u;
            double v;
            double q;
            for (int i = 0; i < mm; i++)
            {
                do
                {
                    u = doub();
                    v = 1.7156 * (doub() - 0.5);
                    double x = u - 0.449871;
                    double y = Math.Abs(v) + 0.386595;
                    q = Globals.SQR(x) + y * (0.19600 * y - 0.25472 * x);
                } while (q > 0.27597 && (q > 0.27846 || Globals.SQR(v) > -4.0 * Math.Log(u) * Globals.SQR(u)));
                spt[i] = v / u;
            }
            chol.elmult(spt, pt);
            for (int i = 0; i < mm; i++)
            {
                pt[i] += mean[i];
            }
            return pt;
        }

    }
}
 

2 代码格式

using System;namespace Legalsoft.Truffer
{public class Multinormaldev : Ran{public Cholesky chol { get; set; } = null;private int mm { get; set; }private double[] mean { get; set; }private double[,] xvar { get; set; }private double[] spt { get; set; }private double[] pt { get; set; }public Multinormaldev(ulong j, double[] mmean, double[,] vvar) : base(j){this.mm = mmean.Length;this.mean = Globals.CopyFrom(mmean);this.xvar = Globals.CopyFrom(vvar);this.chol = new Cholesky(xvar);this.spt = new double[mm];this.pt = new double[mm];if (xvar.GetLength(1) != mm || xvar.GetLength(0) != mm){throw new Exception("bad sizes");}}public double[] dev(){double u;double v;double q;for (int i = 0; i < mm; i++){do{u = doub();v = 1.7156 * (doub() - 0.5);double x = u - 0.449871;double y = Math.Abs(v) + 0.386595;q = Globals.SQR(x) + y * (0.19600 * y - 0.25472 * x);} while (q > 0.27597 && (q > 0.27846 || Globals.SQR(v) > -4.0 * Math.Log(u) * Globals.SQR(u)));spt[i] = v / u;}chol.elmult(spt, pt);for (int i = 0; i < mm; i++){pt[i] += mean[i];}return pt;}}
}