--- /dev/null
+using System.Threading;\r
+using System;\r
+using BigNum;\r
+\r
+namespace Psychlops\r
+{\r
+\r
+ namespace Solver\r
+ {\r
+ internal static class CONST\r
+ {\r
+ public static readonly uint MAX_ARG = 5;\r
+ }\r
+\r
+ \r
+ static public class Constants\r
+ {\r
+ public static readonly int LIMIT = 5;\r
+ public delegate double Func0(double x);\r
+ public delegate double Func1(double x, double a);\r
+ public delegate double Func2(double x, double a, double b);\r
+ public delegate double Func3(double x, double a, double b, double c);\r
+ public delegate double Func4(double x, double a, double b, double c, double d);\r
+ public delegate double Func5(double x, double a, double b, double c, double d, double e);\r
+ }\r
+\r
+\r
+ static public class Binomial\r
+ {\r
+ public static double factorial(int fact)\r
+ {\r
+ double x = 1;\r
+ for (int i = 1; i <= fact; i++)\r
+ {\r
+ x *= i;\r
+ }\r
+ return x;\r
+ }\r
+ public static double combination(int nn, int kk)\r
+ {\r
+ return (double)(Int64)(factorial(nn) / (factorial(kk) * factorial(nn - kk)));\r
+ /*\r
+ BigInt n = new BigInt((Int64)nn, new PrecisionSpec());\r
+ BigInt k = new BigInt((Int64)kk, new PrecisionSpec());\r
+ BigInt n_k = n - k;\r
+ n.Factorial();\r
+ k.Factorial();\r
+ n_k.Factorial();\r
+ return (double)(Int64)(n / (k * n_k) );\r
+ * */\r
+ }\r
+ public static double likelihood(double pp, int yes, int no)\r
+ {\r
+ double p;\r
+ if (pp < 0) p = 0; else if (pp > 1) p = 1; else p = pp;\r
+ return combination(yes + no, yes) * System.Math.Pow(p, yes) * System.Math.Pow(1 - p, no);\r
+ }\r
+\r
+ public static double likelihood(double pp, int yes, int no, double comb)\r
+ {\r
+ double p;\r
+ if (pp < 0) p = 0; else if (pp > 1) p = 1; else p = pp;\r
+ return comb * System.Math.Pow(p, yes) * System.Math.Pow(1 - p, no);\r
+ }\r
+\r
+ }\r
+\r
+\r
+ public class BernoulliProcess\r
+ {\r
+ public struct Data\r
+ {\r
+ public double x;\r
+ public int pos, neg, comb;\r
+ public void set(double cond, int posit, int negat)\r
+ {\r
+ x = cond;\r
+ pos = posit;\r
+ neg = negat;\r
+ }\r
+\r
+ };\r
+ public Data[] elems;\r
+ public int length;\r
+ public void set(double[][] num)\r
+ {\r
+ elems = new Data[num.GetLength(0)];\r
+ for (int i = 0; i < elems.GetLength(0); i++)\r
+ {\r
+ elems[i].set(num[i][0], (int)num[i][1], (int)num[i][2]);\r
+ }\r
+ }\r
+ };\r
+\r
+\r
+\r
+ class BinomialLikelihoodThread\r
+ {\r
+ public bool finished;\r
+ public static bool started;\r
+ internal Thread th;\r
+ public static BinomialLikelihoodThread current;\r
+\r
+ public Interval[] itvl;\r
+ public double[] step;\r
+ public double[] champ;\r
+ public double champ_like;\r
+ public BernoulliProcess data;\r
+\r
+\r
+ public Constants.Func0 func0;\r
+ public Constants.Func1 func1;\r
+ public Constants.Func2 func2;\r
+ public Constants.Func3 func3;\r
+ public Constants.Func4 func4;\r
+ public Constants.Func5 func5;\r
+\r
+\r
+ public BinomialLikelihoodThread()\r
+ {\r
+ itvl = new Interval[CONST.MAX_ARG];\r
+ step = new double[CONST.MAX_ARG];\r
+ champ = new double[CONST.MAX_ARG];\r
+ data = new BernoulliProcess();\r
+ }\r
+ public void waitLoop()\r
+ {\r
+ finished = false;\r
+ for (int i = 0; i < CONST.MAX_ARG; i++)\r
+ {\r
+ champ[i] = 0;\r
+ }\r
+ current = this;\r
+ started = false;\r
+ }\r
+\r
+ public void loop1() { waitLoop(); th = new Thread(new ThreadStart(loop1_)); th.Start(); }\r
+ public void loop2() { waitLoop(); th = new Thread(new ThreadStart(loop2_)); th.Start(); }\r
+ public void loop3() { waitLoop(); th = new Thread(new ThreadStart(loop3_)); th.Start(); }\r
+ public void loop4() { waitLoop(); th = new Thread(new ThreadStart(loop4_)); th.Start(); }\r
+ public void loop5() { waitLoop(); th = new Thread(new ThreadStart(loop5_)); th.Start(); }\r
+ void loop1_()\r
+ {\r
+ started = true;\r
+ double p;\r
+ double like = 1.0;\r
+ champ_like=0.0;\r
+ int L = data.length;\r
+ for (double a = itvl[0].begin.value; a < itvl[0].end.value; a += step[0])\r
+ {\r
+ like = 1.0;\r
+ for(int i=0; i<L; i++) {\r
+ p = func1(data.elems[i].x, a);\r
+ like *= Binomial.likelihood(p, data.elems[i].pos, data.elems[i].neg, data.elems[i].comb);\r
+ }\r
+ if(like > champ_like) {\r
+ champ_like = like;\r
+ champ[0] = a;\r
+ }\r
+ }\r
+ finished = true;\r
+ }\r
+ void loop2_()\r
+ {\r
+ started = true;\r
+ double p = 0.0;\r
+ double like = 1.0;\r
+ champ_like = 0.0;\r
+ int L = data.length;\r
+ for (double a = itvl[0].begin.value; a < itvl[0].end.value; a += step[0])\r
+ {\r
+ for (double b = itvl[1].begin.value; b < itvl[1].end.value; b += step[1])\r
+ {\r
+ like = 1.0;\r
+ for (int i = 0; i < L; i++)\r
+ {\r
+ p = func2(data.elems[i].x, a, b);\r
+ like *= Binomial.likelihood(p, data.elems[i].pos, data.elems[i].neg, data.elems[i].comb);\r
+ }\r
+ if (like > champ_like)\r
+ {\r
+ champ_like = like;\r
+ champ[0] = a;\r
+ champ[1] = b;\r
+ }\r
+ }\r
+ }\r
+ finished = true;\r
+ }\r
+ void loop3_()\r
+ {\r
+ started = true;\r
+ double p;\r
+ double like = 1.0;\r
+ champ_like = 0.0;\r
+ int L = data.length;\r
+ for (double a = itvl[0].begin.value; a < itvl[0].end.value; a += step[0])\r
+ {\r
+ for (double b = itvl[1].begin.value; b < itvl[1].end.value; b += step[1])\r
+ {\r
+ for (double c = itvl[2].begin.value; c < itvl[2].end.value; c += step[2])\r
+ {\r
+ like = 1.0;\r
+ for (int i = 0; i < L; i++)\r
+ {\r
+ p = func3(data.elems[i].x, a, b, c);\r
+ like *= Binomial.likelihood(p, data.elems[i].pos, data.elems[i].neg, data.elems[i].comb);\r
+ }\r
+ if (like > champ_like)\r
+ {\r
+ champ_like = like;\r
+ champ[0] = a;\r
+ champ[1] = b;\r
+ champ[1] = c;\r
+ }\r
+ }\r
+ }\r
+ }\r
+ finished = true;\r
+ }\r
+ void loop4_()\r
+ {\r
+ started = true;\r
+ finished = true;\r
+ }\r
+ void loop5_()\r
+ {\r
+ started = true;\r
+ finished = true;\r
+ }\r
+ };\r
+\r
+\r
+\r
+ public class BinomialLikelihood\r
+ {\r
+ /*\r
+ public static void showWindow(Constants.Func1 f)\r
+ {\r
+ Main.canvas.api_canvas.Dispatcher.BeginInvoke(new Action<Constants.Func1>(showWindow_), f);\r
+ }\r
+ internal static void showWindow_(Constants.Func1 f)\r
+ {\r
+ System.Windows.Controls.ChildWindow page = new PsychlopsSilverlight4.Pages.BinomialSolver(f);\r
+ page.Show();\r
+ }\r
+ public static void showWindow(Constants.Func2 f)\r
+ {\r
+ Main.canvas.api_canvas.Dispatcher.BeginInvoke(new Action<Constants.Func2>(showWindow_), f);\r
+ }\r
+ internal static void showWindow_(Constants.Func2 f)\r
+ {\r
+ System.Windows.Controls.ChildWindow page = new PsychlopsSilverlight4.Pages.BinomialSolver(f);\r
+ page.Show();\r
+ }\r
+ public static void showWindow(Constants.Func3 f)\r
+ {\r
+ Main.canvas.api_canvas.Dispatcher.BeginInvoke(new Action<Constants.Func3>(showWindow_), f);\r
+ }\r
+ internal static void showWindow_(Constants.Func3 f)\r
+ {\r
+ System.Windows.Controls.ChildWindow page = new PsychlopsSilverlight4.Pages.BinomialSolver(f);\r
+ page.Show();\r
+ }\r
+ */\r
+\r
+ public int iteration;\r
+\r
+ public Interval[] itvl;\r
+ public double[] step;\r
+ public double[] champ;\r
+ public double champ_like;\r
+ public BernoulliProcess data;\r
+\r
+ public Constants.Func0 func0;\r
+ public Constants.Func1 func1;\r
+ public Constants.Func2 func2;\r
+ public Constants.Func3 func3;\r
+ public Constants.Func4 func4;\r
+ public Constants.Func5 func5;\r
+\r
+ public BinomialLikelihood()\r
+ {\r
+ itvl = new Interval[CONST.MAX_ARG];\r
+ step = new double[CONST.MAX_ARG];\r
+ champ = new double[CONST.MAX_ARG];\r
+ iteration = 2;\r
+ }\r
+\r
+ public void begin(Constants.Func0 d_func)\r
+ { }\r
+ public void begin(Constants.Func1 d_func)\r
+ {\r
+ func1 = d_func;\r
+\r
+ BinomialLikelihoodThread[] l = new BinomialLikelihoodThread[4];\r
+ for (int i = 0; i < 4; i++) { l[i] = new BinomialLikelihoodThread(); }\r
+\r
+ for (int k = 0; k < iteration; k++)\r
+ {\r
+ begin_base(l);\r
+ for (int i = 0; i < 4; i++)\r
+ {\r
+ l[i].data = data;\r
+ l[i].func1 = d_func;\r
+ l[i].loop1();\r
+ }\r
+ end_base(l);\r
+ }\r
+ }\r
+ public void begin(Constants.Func2 d_func)\r
+ {\r
+ func2 = d_func;\r
+\r
+ BinomialLikelihoodThread[] l = new BinomialLikelihoodThread[4];\r
+ for (int i = 0; i < 4; i++) { l[i] = new BinomialLikelihoodThread(); }\r
+\r
+ for (int k = 0; k < iteration; k++)\r
+ {\r
+ begin_base(l);\r
+ for (int i = 0; i < 4; i++)\r
+ {\r
+ l[i].data = data;\r
+ l[i].func2 = d_func;\r
+ l[i].loop2();\r
+ }\r
+ end_base(l);\r
+ }\r
+ }\r
+ public void begin(Constants.Func3 d_func)\r
+ {\r
+ func3 = d_func;\r
+\r
+ BinomialLikelihoodThread[] l = new BinomialLikelihoodThread[4];\r
+ for (int i = 0; i < 4; i++) { l[i] = new BinomialLikelihoodThread(); }\r
+\r
+ for (int k = 0; k < iteration; k++)\r
+ {\r
+ begin_base(l);\r
+ for(int i=0; i<4; i++) {\r
+ l[i].data = data;\r
+ l[i].func3 = d_func;\r
+ l[i].loop3();\r
+ }\r
+ end_base(l);\r
+ }\r
+ }\r
+\r
+ public void begin(Constants.Func4 d_func)\r
+ {\r
+ }\r
+ public void begin(Constants.Func5 d_func)\r
+ {\r
+ }\r
+\r
+ void begin_base(BinomialLikelihoodThread[] l)\r
+ {\r
+ champ_like = 0;\r
+\r
+ data.length = data.elems.GetLength(0);\r
+ for (int i = 0; i < data.elems.GetLength(0); i++)\r
+ {\r
+ data.elems[i].comb = (int)Binomial.combination(data.elems[i].pos + data.elems[i].neg, data.elems[i].pos);\r
+ }\r
+\r
+ for (int j = 0; j < Constants.LIMIT; j++) { step[j] = (itvl[j].end.value - itvl[j].begin.value) / 256.0; }\r
+ for (int i = 0; i < 4; i++)\r
+ {\r
+ l[i].itvl[0] = new Interval((itvl[0].begin.value) + (i * step[0] * 64), (itvl[0].begin.value) + ((i + 1) * step[0] * 64));\r
+ l[i].step[0] = step[0];\r
+ for (int j = 1; j < Constants.LIMIT; j++)\r
+ {\r
+ l[i].itvl[j] = itvl[j];\r
+ l[i].step[j] = step[j];\r
+ }\r
+ l[i].data = data;\r
+ }\r
+ }\r
+\r
+ void end_base(BinomialLikelihoodThread[] l)\r
+ {\r
+ //for (int i = 0; i < 4; i++) { l[i].th.Join(); }\r
+ while (!l[0].finished || !l[1].finished || !l[2].finished || !l[3].finished) { Thread.Sleep(100); } \r
+\r
+ for(int j=0; j<Constants.LIMIT; j++) { champ[j] = 0; }\r
+ champ_like = 0.0;\r
+ for(int i=0; i<4; i++) {\r
+ if(champ_like < l[i].champ_like) {\r
+ champ_like = l[i].champ_like;\r
+ for(int j=0; j<Constants.LIMIT; j++) { champ[j] = l[i].champ[j]; }\r
+ }\r
+ }\r
+\r
+ double r, low, high;\r
+ for (int j = 0; j < Constants.LIMIT; j++)\r
+ {\r
+ r = itvl[j].end.value - itvl[j].begin.value;\r
+ low = champ[j] - r / 8.0 < itvl[j].begin.value ? itvl[j].begin.value : champ[j] - r / 8.0;\r
+ high = champ[j] + r / 8.0 > itvl[j].end.value ? itvl[j].end.value : champ[j] + r / 8.0;\r
+ itvl[j] = new Interval(low, high);\r
+ }\r
+\r
+ }\r
+\r
+\r
+ }\r
+\r
+\r
+ }\r
+\r
+}
\ No newline at end of file
OSDN Git Service
