// g02hd Example Program Text // C# version, NAG Copyright 2008 using System; using NagLibrary; using System.Globalization; using System.IO; namespace NagDotNetExamples { public class G02HDE { static string datafile = "ExampleData/g02hde.d"; static void Main(String[] args) { if (args.Length == 1) { datafile = args[0]; } StartExample(); } public static void StartExample() { try { G02.G02HD_CHI chiG02HD = new G02.G02HD_CHI(chi); G02.G02HD_PSI psiG02HD = new G02.G02HD_PSI(psi); DataReader sr = new DataReader(datafile); double beta, eps, psip0, sigma, tol; int i, indw, isigma, j, k, m, maxit, n, nit, nitmon; int ifail; Console.WriteLine("g02hd Example Program Results"); // Skip heading in data file beta = 0.0; sr.Reset(); // Read in the dimensions of X sr.Reset(); n = int.Parse(sr.Next()); m = int.Parse(sr.Next()); double[] rs = new double[n]; double[] theta = new double[m]; double[] wgt = new double[n]; double[,] x = new double[n, m]; double[] y = new double[n]; if (n > 0 && m > 0) { // Read in the x matrix, the y values and set x[i-1,0] to 1 for the // constant term for (i = 1; i <= n; i++) { sr.Reset(); for (j = 2; j <= m; j++) { x[i - 1, j - 1] = double.Parse(sr.Next(), CultureInfo.InvariantCulture); } y[i - 1] = double.Parse(sr.Next(), CultureInfo.InvariantCulture); x[i - 1, 0] = 1.00e0; } // Read in weights for (i = 1; i <= n; i++) { sr.Reset(); wgt[i - 1] = double.Parse(sr.Next(), CultureInfo.InvariantCulture); } betcal(ref n, wgt, ref beta); // Set other parameter values maxit = 50; tol = 0.50e-4; eps = 0.50e-5; psip0 = 1.00e0; // Set value of isigma and initial value of sigma isigma = 1; sigma = 1.00e0; // Set initial value of theta for (j = 1; j <= m; j++) { theta[j - 1] = 0.00e0; } // Change nitmon to a positive value if monitoring information // is required nitmon = 0; // Schweppe type regression indw = 1; // G02.g02hd(chiG02HD, psiG02HD, psip0, beta, indw, isigma, n, m, x, y, wgt, theta, out k, ref sigma, rs, tol, eps, maxit, nitmon, out nit, out ifail); // Console.WriteLine(" "); if ((ifail != 0) && (ifail != 7)) { Console.WriteLine(" "); Console.WriteLine("** g02hd failed with ifail = {0,5}", ifail); } else { if (ifail == 7) { Console.Write(" {0}{1,5}", "g02hd returned ifail = ", ifail); Console.Write(" {0}", "Some of the following results may be unreliable"); } Console.WriteLine(" {0}{1,4}{2}", "g02hd required ", nit, " iterations to converge"); Console.WriteLine(" {0}{1,4}", " K = ", k); Console.WriteLine(" {0}{1,9:f4}", " Sigma = ", sigma); Console.WriteLine(" {0}", " THETA"); for (j = 1; j <= m; j++) { Console.WriteLine(" {0,9:f4}", theta[j - 1]); } Console.WriteLine(" "); Console.WriteLine(" {0}", " Weights Residuals"); for (i = 1; i <= n; i++) { Console.WriteLine(" {0,9:f4}{1,9:f4}", wgt[i - 1], rs[i - 1]); } } } // } catch (Exception e) { Console.WriteLine(e.Message); Console.Write("Exception Raised"); } } // public static double psi(double t) { double psiValue = 0.0; double c = 1.5e0; if (t <= -c) { psiValue = -c; } else if (Math.Abs(t) < c) { psiValue = t; } else { psiValue = c; } return psiValue; } public static double chi(double t) { double chiValue = 0.0; double ps = 0.0; double dchi = 1.5e0; ps = dchi; if (Math.Abs(t) < dchi) { ps = t; } chiValue = ps * ps / 2.00e0; return chiValue; } // public static void betcal(ref int n, double[] wgt, ref double beta) { // Calculate beta for Schweppe type regression const double dchi = 1.50e0; double amaxex, anormc, b, d2, dc, dw, dw2, pc, w2; int i; amaxex = -Math.Log(X02.x02ak()); anormc = Math.Sqrt(X01.x01aa() * 2.00e0); d2 = dchi * dchi; beta = 0.00e0; for (i = 1; i <= n; i++) { w2 = wgt[i - 1] * wgt[i - 1]; dw = wgt[i - 1] * dchi; pc = S.s15ab(dw); dw2 = dw * dw; dc = 0.00e0; if (dw2 < amaxex) { dc = Math.Exp(-dw2 / 2.00e0) / anormc; } b = (-dw * dc + pc - 0.50e0) / w2 + (1.00e0 - pc) * d2; beta = b * w2 / (double)n + beta; } } } }