myghxw[n_Integer] := Module[
    {xis, wis, getw, x},

    (* get xi *)
    xis = NSolve[HermiteH[n, x] == 0];
    xis = Table[xis[[i]][[1]][[2]], {i, 1, Length[xis]}];


    (* get wi *)
    getw[xi_] := 2^(n - 1)*n!*Sqrt[Pi]/(n^2*HermiteH[n - 1, xi]^2);
    wis = getw[#] & /@ xis;

    {xis, wis}

]

{xs,ws}=myghxw[20];


voles1={29, 15, 15, 16, 27};
kvoles1= 5;

(* m7 = Log Gamma Model *)
mylikLogGamma[data_, kdata_, f0_?NumericQ, alpha_?NumericQ, beta_?NumericQ] := Module[
    {K, pj, fj, N0, loglik, above, below},

    K = kdata;

    (* Using Gauss Hermite *)
    (* Evaluted at each xs *)

    pj = Table[
         Exp[-xs[[i]]^2./beta]^j
         *
         (1.-Exp[-xs[[i]]^2./beta])^(K-j)
         *
         xs[[i]]^(2.*alpha-2.)
         , {i, 1, Length[xs]}
         , {j, 0, K}
       ];

    (* Put weights on *)
    pj = ws.pj;

    pj = Table[Binomial[K, j], {j, 0, K}]*pj;

    fj = Prepend[PadRight[data,K], f0];

    N0 = Sum[fj[[j]], {j, 1, Length[fj]}] ;
    above = LogGamma[N0 + 1] ;
    below = Sum[LogGamma[fj[[j]] + 1], {j, 1, Length[fj]}] ;
    loglik = fj.Log[pj] ;

    loglik = above - below + loglik;

    N[loglik]

]

(* Test for fitting*)

myLogGammaFit[data_, kdata_]:= Module[
    {testfun,ans},

    testfun[f0_, alpha_, beta_]:= Module[
       {ans2},

       ans2 = mylikLogGamma[data,kdata, f0, alpha, beta];

       If[ Im[ans2] !=0, ans2=-9999];

       ans2
    ];

    ans=NMaximize[
       {testfun[f0,alpha,beta], 
         f0 >= 0
         && alpha>0
         && beta>0
       },
       {f0, alpha, beta},
       Method -> {"RandomSearch", "SearchPoints" -> 200}
    ]//Quiet;

    ans

]


myLogGammaFit[voles1,kvoles1]

Spent two days, still can't figure our WHY????

It does not seem to be evaluating???