(*******************************************************************
This file was generated automatically by the Mathematica front end.
It contains Initialization cells from a Notebook file, which
typically will have the same name as this file except ending in
".nb" instead of ".m".

This file is intended to be loaded into the Mathematica kernel using
the package loading commands Get or Needs.  Doing so is equivalent
to using the Evaluate Initialization Cells menu command in the front
end.

DO NOT EDIT THIS FILE.  This entire file is regenerated
automatically each time the parent Notebook file is saved in the
Mathematica front end.  Any changes you make to this file will be
overwritten.
***********************************************************************)

<<Statistics`ContinuousDistributions`
<<Statistics`ClusterAnalysis`
<<Graphics`Graphics`

\!\(\*
  RowBox[{
    RowBox[{\(mclustv[data_, ind_, g_, c_, c0_]\), ":=", 
      RowBox[{"Module", "[", 
        RowBox[{\({clust, var, n, npdf, phi, p, cond1, tau, taumat, num, bic, aic, TM, kai, term1, term2, termmu, termsig, ll, EN, iclbic, aic3}\), ",", "\[IndentingNewLine]", 
          RowBox[{
            RowBox[{"clust", "=", 
              StyleBox[\(FindClusters[data, g]\),
                FormatType->StandardForm]}], 
            StyleBox[";",
              FormatType->StandardForm], 
            StyleBox["\[IndentingNewLine]",
              FormatType->StandardForm], \(z[0] = Table[Length[clust[\([i]\)]]\/Length[data], {i, 1, g}]\), ";", "\[IndentingNewLine]", \(mu[0] = Table[Mean[clust[\([i]\)]], {i, 1, g}]\), ";", "\[IndentingNewLine]", \(var[x_] := Variance[x] /; Length[x] > 1\), ";", "\[IndentingNewLine]", \(var[x_] := 10. ^\((\(-2\))\) /; Length[x] \[Equal] 1\), ";", "\[IndentingNewLine]", \(sig2[0] = Table[2*var[clust[\([i]\)]], {i, 1, g}]\), ";", "\[IndentingNewLine]", \(n = Length[data]\), ";", "\[IndentingNewLine]", \(npdf[x_, a1_, b1_] = PDF[NormalDistribution[a1, \@b1], x]\), ";", "\[IndentingNewLine]", \(phi[z_] = PDF[NormalDistribution[0, 1], z]\), ";", "\[IndentingNewLine]", \(p[i_, j_, k_] := npdf[data[\([i]\)], \(mu[k]\)[\([j]\)], \(sig2[k]\)[\([j]\)]]*\((1 - ind[\([i, 1]\)])\)*\((1 - ind[\([i, 2]\)])\) + 0.5*Erf[1\/\@2*\((\(c - \(mu[k]\)[\([j]\)]\)\/\@\(sig2[k]\)[\([j]\)])\), \[Infinity]]*ind[\([i, 1]\)]*\((1 - ind[\([i, 2]\)])\) + 0.5*N[Erf[\(-\[Infinity]\), 1\/\@2*\((\(c0 - \(mu[k]\)[\([j]\)]\)\/\@\(sig2[k]\)[\([j]\)])\)]]*\((1 - ind[\([i, 1]\)])\)*ind[\([i, 2]\)]\), ";", "\[IndentingNewLine]", \(Loglikeli[k_] := \[Sum]\+\(i = 1\)\%n Log[\[Sum]\+\(j = 1\)\%g\( z[k]\)[\([j]\)]*p[i, j, k]]\), ";", "\[IndentingNewLine]", \(cond1[k_] := Abs[\((Loglikeli[k]/Loglikeli[k - 1])\) - 1] /; k \[GreaterEqual] 1\), ";", "\[IndentingNewLine]", \(cond1[0] = 1\), ";", "\[IndentingNewLine]", \(tau[i_, j_, k_] := \(p[i, j, k]*\(z[k]\)[\([j]\)]\)\/\(\[Sum]\+\(s = 1\)\%g p[i, s, k]*\(z[k]\)[\([s]\)]\)\), ";", "\[IndentingNewLine]", \(taumat[k_] := Table[tau[i, j, k], {i, 1, n}, {j, 1, g}]\), ";", "\[IndentingNewLine]", \(term1[k_, j_] := Exp[\(-Log[Erf[1\/\@2*\((\(c - \(mu[k]\)[\([j]\)]\)\/\@\(sig2[k]\)[\([j]\)])\), \[Infinity]]]\) - 0.5*Log[\[Pi]\/2] - 0.5 \((\(c - \(mu[k]\)[\([j]\)]\)\/\@\(sig2[k]\)[\([j]\)])\)\^2]\), ";", "\[IndentingNewLine]", \(term2[k_, j_] := Exp[\(-Log[N[Erf[\(-\[Infinity]\), 1\/\@2*\((\(c0 - \(mu[k]\)[\([j]\)]\)\/\@\(sig2[k]\)[\([j]\)])\)]]]\) - 0.5*Log[\[Pi]\/2] - 0.5 \((\(c0 - \(mu[k]\)[\([j]\)]\)\/\@\(sig2[k]\)[\([j]\)])\)\^2]\), ";", "\[IndentingNewLine]", \(termmu[i_, j_, k_] := data[\([i]\)] /; \((1 - ind[\([i, 1]\)])\)*\((1 - ind[\([i, 2]\)])\) \[Equal] 1\), ";", "\[IndentingNewLine]", \(termmu[i_, j_, k_] := \((\(mu[k]\)[\([j]\)] + \@\(sig2[k]\)[\([j]\)]*term1[k, j])\) /; ind[\([i, 1]\)] \[Equal] 1\), ";", "\[IndentingNewLine]", \(termmu[i_, j_, k_] := \((\(mu[k]\)[\([j]\)] - \@\(sig2[k]\)[\([j]\)]*term2[k, j])\) /; ind[\([i, 2]\)] \[Equal] 1\), ";", "\[IndentingNewLine]", \(termsig[i_, j_, k_] := \((data[\([i]\)] - \(mu[k + 1]\)[\([j]\)])\)\^2 /; \((1 - ind[\([i, 1]\)])\)*\((1 - ind[\([i, 2]\)])\) \[Equal] 1\), ";", "\[IndentingNewLine]", \(termsig[i_, j_, k_] := \((\(mu[k + 1]\)[\([j]\)] - \(mu[k]\)[\([j]\)])\)\^2 - 2*\((\(mu[k + 1]\)[\([j]\)] - \(mu[k]\)[\([j]\)])\)*\@\(sig2[k]\)[\([j]\)]*term1[k, j] + \(sig2[k]\)[\([j]\)]*\((1 + \(c - \(mu[k]\)[\([j]\)]\)\/\@\(sig2[k]\)[\([j]\)]*term1[k, j])\) /; ind[\([i, 1]\)] \[Equal] 1\), ";", "\[IndentingNewLine]", \(termsig[i_, j_, k_] := \((\(mu[k + 1]\)[\([j]\)] - \(mu[k]\)[\([j]\)])\)\^2 + 2*\((\(mu[k + 1]\)[\([j]\)] - \(mu[k]\)[\([j]\)])\)*\@\(sig2[k]\)[\([j]\)]*term2[k, j] + \(sig2[k]\)[\([j]\)]*\((1 - \(c0 - \(mu[k]\)[\([j]\)]\)\/\@\(sig2[k]\)[\([j]\)]*term2[k, j])\) /; ind[\([i, 2]\)] \[Equal] 1\), ";", "\[IndentingNewLine]", \(num = 100\), ";", "\[IndentingNewLine]", \(For[k = 0, k \[LessEqual] num && cond1[k] > 10^\((\(-7\))\) && Min[sig2[k]] > 10^\((\(-6\))\), k = k + 1, \[IndentingNewLine]TM = taumat[k]; \[IndentingNewLine]z[k + 1] = Table[\(1\/n\) \(\[Sum]\+\(i = 1\)\%n TM[\([i, j]\)]\), {j, 1, g}]; \[IndentingNewLine]mu[k + 1] = Table[\(\[Sum]\+\(i = 1\)\%n TM[\([i, j]\)]*termmu[i, j, k]\)\/\(\[Sum]\+\(i = 1\)\%n TM[\([i, j]\)]\), {j, 1, g}]; \[IndentingNewLine]\[IndentingNewLine]sig2[k + 1] = Table[\(\[Sum]\+\(i = 1\)\%n TM[\([i, j]\)]*termsig[i, j, k]\)\/\(\[Sum]\+\(i = 1\)\%n TM[\([i, j]\)]\), {j, 1, g}]; \[IndentingNewLine]rep = k + 1;]\), ";", "\[IndentingNewLine]", \(ll = Loglikeli[rep]\), ";", "\[IndentingNewLine]", \(bic = \(-2\)*ll + \((Log[n])\)*\((Length[mu[0]] + Length[sig2[0]] + Length[z[0]] - 1)\)\), ";", "\[IndentingNewLine]", \(aic = \(-2\)*ll + 2*\((Length[mu[0]] + Length[sig2[0]] + Length[z[0]] - 1)\)\), ";", "\[IndentingNewLine]", \(aic3 = aic + \((Length[mu[0]] + Length[sig2[0]] + Length[z[0]] - 1)\)\), ";", "\[IndentingNewLine]", \(EN = \(-\(\[Sum]\+\(i = 1\)\%n\(\[Sum]\+\(j = 1\)\%g tau[i, j, rep]*Log[tau[i, j, rep]]\)\)\)\), ";", "\[IndentingNewLine]", \(iclbic = bic + 2*EN\), ";", "\[IndentingNewLine]", \(If[Min[sig2[rep]] < 10^\((\(-6\))\), {z[rep], mu[rep], sig2[rep], bic, aic, aic3, iclbic, 0}, {z[rep], mu[rep], sig2[rep], bic, aic, aic3, iclbic, 1}]\)}]}], "]"}]}], ";"}]\)

\!\(\*
  RowBox[{
    RowBox[{\(mcluste[data_, ind_, g_, c_, c0_]\), ":=", 
      RowBox[{"Module", "[", 
        RowBox[{\({clust, var, n, npdf, phi, p, cond1, tau, taumat, num, bic, aic, aic3, TM, kai, term1, term2, termmu, termsig, ll, EN, iclbic}\), ",", "\[IndentingNewLine]", 
          RowBox[{
            RowBox[{"clust", "=", 
              StyleBox[\(FindClusters[data, g]\),
                FormatType->StandardForm]}], 
            StyleBox[";",
              FormatType->StandardForm], 
            StyleBox["\[IndentingNewLine]",
              FormatType->StandardForm], \(z[0] = Table[Length[clust[\([i]\)]]\/Length[data], {i, 1, g}]\), ";", "\[IndentingNewLine]", \(mu[0] = Table[Mean[clust[\([i]\)]], {i, 1, g}]\), ";", "\[IndentingNewLine]", \(var[x_] := Variance[x] /; Length[x] > 1\), ";", "\[IndentingNewLine]", \(var[x_] := 10. ^\((\(-2\))\) /; Length[x] \[Equal] 1\), ";", "\[IndentingNewLine]", \(sig2[0] = \(1\/\(Length[data] - 1\)\) \(\[Sum]\+\(i = 1\)\%g\((Length[clust[\([i]\)]] - 1)\)*var[clust[\([i]\)]]\)\), ";", "\[IndentingNewLine]", \(n = Length[data]\), ";", "\[IndentingNewLine]", \(npdf[x_, a1_, b1_] = PDF[NormalDistribution[a1, \@b1], x]\), ";", "\[IndentingNewLine]", \(phi[z_] = PDF[NormalDistribution[0, 1], z]\), ";", "\[IndentingNewLine]", \(p[i_, j_, k_] := npdf[data[\([i]\)], \(mu[k]\)[\([j]\)], sig2[k]]*\((1 - ind[\([i, 1]\)])\)*\((1 - ind[\([i, 2]\)])\) + 0.5*Erf[1\/\@2*\((\(c - \(mu[k]\)[\([j]\)]\)\/\@sig2[k])\), \[Infinity]]*ind[\([i, 1]\)]*\((1 - ind[\([i, 2]\)])\) + 0.5*N[Erf[\(-\[Infinity]\), 1\/\@2*\((\(c0 - \(mu[k]\)[\([j]\)]\)\/\@sig2[k])\)]]*\((1 - ind[\([i, 1]\)])\)*ind[\([i, 2]\)]\), ";", "\[IndentingNewLine]", \(Loglikeli[k_] := \[Sum]\+\(i = 1\)\%n Log[\[Sum]\+\(j = 1\)\%g\( z[k]\)[\([j]\)]*p[i, j, k]]\), ";", "\[IndentingNewLine]", \(cond1[k_] := Abs[\((Loglikeli[k]/Loglikeli[k - 1])\) - 1] /; k \[GreaterEqual] 1\), ";", "\[IndentingNewLine]", \(cond1[0] = 1\), ";", "\[IndentingNewLine]", \(tau[i_, j_, k_] := \(p[i, j, k]*\(z[k]\)[\([j]\)]\)\/\(\[Sum]\+\(s = 1\)\%g p[i, s, k]*\(z[k]\)[\([s]\)]\)\), ";", "\[IndentingNewLine]", \(taumat[k_] := Table[tau[i, j, k], {i, 1, n}, {j, 1, g}]\), ";", "\[IndentingNewLine]", \(term1[k_, j_] := Exp[\(-Log[Erf[1\/\@2*\((\(c - \(mu[k]\)[\([j]\)]\)\/\@sig2[k])\), \[Infinity]]]\) - 0.5*Log[\[Pi]\/2] - 0.5 \((\(c - \(mu[k]\)[\([j]\)]\)\/\@sig2[k])\)\^2]\), ";", "\[IndentingNewLine]", \(term2[k_, j_] := Exp[\(-Log[N[Erf[\(-\[Infinity]\), 1\/\@2*\((\(c0 - \(mu[k]\)[\([j]\)]\)\/\@sig2[k])\)]]]\) - 0.5*Log[\[Pi]\/2] - 0.5 \((\(c0 - \(mu[k]\)[\([j]\)]\)\/\@sig2[k])\)\^2]\), ";", "\[IndentingNewLine]", \(termmu[i_, j_, k_] := data[\([i]\)] /; \((1 - ind[\([i, 1]\)])\)*\((1 - ind[\([i, 2]\)])\) \[Equal] 1\), ";", "\[IndentingNewLine]", \(termmu[i_, j_, k_] := \((\(mu[k]\)[\([j]\)] + \@sig2[k]*term1[k, j])\) /; ind[\([i, 1]\)] \[Equal] 1\), ";", "\[IndentingNewLine]", \(termmu[i_, j_, k_] := \((\(mu[k]\)[\([j]\)] - \@sig2[k]*term2[k, j])\) /; ind[\([i, 2]\)] \[Equal] 1\), ";", "\[IndentingNewLine]", \(termsig[i_, j_, k_] := \((data[\([i]\)] - \(mu[k + 1]\)[\([j]\)])\)\^2 /; \((1 - ind[\([i, 1]\)])\)*\((1 - ind[\([i, 2]\)])\) \[Equal] 1\), ";", "\[IndentingNewLine]", \(termsig[i_, j_, k_] := \((\(mu[k + 1]\)[\([j]\)] - \(mu[k]\)[\([j]\)])\)\^2 - 2*\((\(mu[k + 1]\)[\([j]\)] - \(mu[k]\)[\([j]\)])\)*\@sig2[k]*term1[k, j] + sig2[k]*\((1 + \(c - \(mu[k]\)[\([j]\)]\)\/\@sig2[k]*term1[k, j])\) /; ind[\([i, 1]\)] \[Equal] 1\), ";", "\[IndentingNewLine]", \(termsig[i_, j_, k_] := \((\(mu[k + 1]\)[\([j]\)] - \(mu[k]\)[\([j]\)])\)\^2 + 2*\((\(mu[k + 1]\)[\([j]\)] - \(mu[k]\)[\([j]\)])\)*\@sig2[k]*term2[k, j] + sig2[k]*\((1 - \(c0 - \(mu[k]\)[\([j]\)]\)\/\@sig2[k]*term2[k, j])\) /; ind[\([i, 2]\)] \[Equal] 1\), ";", "\[IndentingNewLine]", \(num = 100\), ";", "\[IndentingNewLine]", \(For[k = 0, k \[LessEqual] num && cond1[k] > 10^\((\(-7\))\), k = k + 1, \[IndentingNewLine]TM = taumat[k]; \[IndentingNewLine]z[k + 1] = Table[\(1\/n\) \(\[Sum]\+\(i = 1\)\%n TM[\([i, j]\)]\), {j, 1, g}]; \[IndentingNewLine]mu[k + 1] = Table[\(\[Sum]\+\(i = 1\)\%n TM[\([i, j]\)]*termmu[i, j, k]\)\/\(\[Sum]\+\(i = 1\)\%n TM[\([i, j]\)]\), {j, 1, g}]; \[IndentingNewLine]sig2[k + 1] = 1\/n*\((\[Sum]\+\(i = 1\)\%n\(\[Sum]\+\(j = 1\)\%g TM[\([i, j]\)]*termsig[i, j, k]\))\); \[IndentingNewLine]rep = k + 1]\), ";", "\[IndentingNewLine]", \(ll = Loglikeli[rep]\), ";", "\[IndentingNewLine]", \(bic = \(-2\)*ll + \((Log[n])\)*\((Length[mu[0]] + Length[z[0]])\)\), ";", "\[IndentingNewLine]", \(aic = \(-2\)*ll + 2*\((Length[mu[0]] + Length[z[0]])\)\), ";", "\[IndentingNewLine]", \(aic3 = aic + \((Length[mu[0]] + Length[z[0]])\)\), ";", "\[IndentingNewLine]", \(EN = \(-\(\[Sum]\+\(i = 1\)\%n\(\[Sum]\+\(j = 1\)\%g tau[i, j, rep]*Log[tau[i, j, rep]]\)\)\)\), ";", "\[IndentingNewLine]", \(iclbic = bic + 2*EN\), ";", "\[IndentingNewLine]", \({z[rep], mu[rep], sig2[rep], bic, aic, aic3, iclbic}\)}]}], "]"}]}], ";"}]\)

hist[data_,s_:"def"]:=Module[{sturges,sturgeswidth,fd,fdwidth,hc},
    sturges[x_]:=(Max[x]-Min[x])/(Log[2,Length[x]]+1);
    sturgeswidth[x_]:=
      Range[Min[x],Min[x]+Ceiling[(Max[x]-Min[x])/sturges[x]]*sturges[x],
        sturges[x]];
    fd[x_]:=2*(Quantile[x,.75]-Quantile[x,.25])/Length[x];
    fdwidth[x_]:=
      Range[Min[x],Min[x]+Ceiling[(Max[x]-Min[x])/fd[x]]*fd[x],fd[x]];
    Which[s===st,hc=sturgeswidth[data],s===frd,hc=fdwidth[data],s===def,
      hc="Automatic"];
    g1v=Histogram[data,HistogramScale\[Rule]1,HistogramRange\[Rule]{0,1},
        HistogramCategories\[Rule]hc,BarStyle\[Rule]GrayLevel[1]];]

\!\(\(npdf2[x_, a2_, b2_] = PDF[NormalDistribution[a2, \@b2], x];\)\[IndentingNewLine]
  \(plugine[x_, para_] := \[Sum]\+\(j = 1\)\%\(Length[para[\([1]\)]]\)para[\([1, j]\)]*npdf2[x, para[\([2, j]\)], para[\([3]\)]];\)\[IndentingNewLine]
  \(pluginv[x_, para_] := \[Sum]\+\(j = 1\)\%\(Length[para[\([1]\)]]\)para[\([1, j]\)]*npdf2[x, para[\([2, j]\)], para[\([3, j]\)]];\)\[IndentingNewLine]
  \(plugin[x_, para_] := Module[{f1}, f1[y_, param_] := \[Sum]\+\(j = 1\)\%\(Length[param[\([1]\)]]\)param[\([1, j]\)]*npdf2[x, param[\([2, j]\)], param[\([3]\)]] /; Length[param[\([3]\)]] \[Equal] 0; \[IndentingNewLine]f1[y_, param_] := \[Sum]\+\(j = 1\)\%\(Length[param[\([1]\)]]\)param[\([1, j]\)]*npdf2[x, param[\([2, j]\)], param[\([3, j]\)]] /; Length[param[\([3]\)]] \[GreaterEqual] 1; \[IndentingNewLine]f1[x, para]];\)\)

\!\(\(qnewtonv[alpha_, para_, initial_] := Module[{kai, g, ncdf, x}, \[IndentingNewLine]g = Length[para[\([1]\)]]; \[IndentingNewLine]ncdf[x_] = CDF[NormalDistribution[0, 1], x]; \[IndentingNewLine]kai = FindRoot[\[Sum]\+\(j = 1\)\%g para[\([1, j]\)]*ncdf[\((x - para[\([2, j]\)])\)/\@para[\([3, j]\)]] - alpha \[Equal] 0, {x, initial}]; \[IndentingNewLine]x /. kai];\)\[IndentingNewLine]
  \(seisekiv[data_, para_, kugiri_] := Module[{percentage1, hyou}, \[IndentingNewLine]percentage1 = Table[Quantile[data, kugiri[\([i]\)]], {i, 1, Length[kugiri]}]; \[IndentingNewLine]hyou = Table[qnewtonv[kugiri[\([i]\)], para, percentage1[\([i]\)]], {i, 1, Length[kugiri]}]; \[IndentingNewLine]hyou];\)\[IndentingNewLine]
  \(qnewton[alpha_, para_, initial_] := Module[{kai, g, ncdf, x}, \[IndentingNewLine]g = Length[para[\([1]\)]]; \[IndentingNewLine]ncdf[x_] = CDF[NormalDistribution[0, 1], x]; \[IndentingNewLine]kai = FindRoot[\[Sum]\+\(j = 1\)\%g para[\([1, j]\)]*ncdf[\((x - para[\([2, j]\)])\)/\@para[\([3]\)]] - alpha \[Equal] 0, {x, initial}]; \[IndentingNewLine]x /. kai];\)\[IndentingNewLine]
  seisekie[data_, para_, kugiri_] := Module[{percentage1, hyou}, \[IndentingNewLine]percentage1 = Table[Quantile[data, kugiri[\([i]\)]], {i, 1, Length[kugiri]}]; \[IndentingNewLine]hyou = Table[qnewton[kugiri[\([i]\)], para, percentage1[\([i]\)]], {i, 1, Length[kugiri]}]; \[IndentingNewLine]hyou]\)