Some problems with complied function

Question

(* Some quantity needed *)

B1997 = {14,10,6,3,3,6,2,3,2,1,1,1,0,3,2,0,3,1,1,1,0,0,1,2,0,0,1};
kB1997 = 50;

(*Gauss Hermite*)
myghxw[n_Integer] := Block[
    {x, 
    xis = Cases[NSolve[HermiteH[n, x] == 0, WorkingPrecision -> 30], _?NumericQ, Infinity]
    },
    {xis, 2^(n - 1)*Gamma[n+1]*Sqrt[Pi]/(n^2*HermiteH[n - 1, xis]^2)}
];
{xsgh, wsgh} = myghxw[30];

Now the compiled function:

Clear[mytest];
mytest["LNB"] = Compile[
    {
        {data, _Real, 1},
        {kdata, _Integer},
        {f0},
        {mu},
        {sigma}
    },
    With[
        {
            fj = Prepend[PadRight[data, kdata], f0], 

    (* Use Gauss Hermite to approximate the integral *)
            pj = wsgh.Table[
                Gamma[kdata + 1]/Gamma[j + 1]/Gamma[kdata - j + 1]/
                Sqrt[Pi]
                *
                (1/(1 + Exp[-xsgh[[i]]*sigma*Sqrt[2] - mu]))^j
                *
                (1 - 1/(1 + Exp[-xsgh[[i]]*sigma*Sqrt[2] - mu]))^(kdata - j), 
            {i, 1, Length[xsgh]}, {j, 0, kdata}]
        },

        LogGamma[Total[fj] + 1] - Total[LogGamma[fj + 1]] + Total[fj*Log[pj]]
    ],
    CompilationOptions -> {"InlineExternalDefinitions" -> True}
];


Needs["CompiledFunctionTools`"]
CompilePrint[mytest["LNB"]]

Still has TWO calls to MainEvaluate.

Is there a way (tricks) to handle this overflow/underflow?

mytest["LNB"][B1997, kB1997, 0.5, -4, 9]

CompiledFunction::cfn: Numerical error encountered at instruction 17; proceeding with uncompiled evaluation. >>

I am not sure why this particular choice of mu=-4 and sigma=9 will have a problem. Even when using NIntegrate, for any other datasets, too. The NIntegrate is a direct way of getting pj. But the complied function version, we use an approximation. But it's a bit wired, the same mu=-4 and sigma=9 will cause problems in both cases. (for any dataset).

For information, the expression for pjs is

Now the running problem:

NMaximize[{mytest["LNB"][B1997, kB1997, f0, mu, sigma], f0 > 0, sigma > 0}, {f0, mu, sigma}, Method -> {"RandomSearch", "SearchPoints" -> 20}]

CompiledFunction::cfsa: Argument f0 at position 3 should be a machine-size real number. >>

(* {-42.3482, {f0 -> 14.8041, mu -> -2.79366, sigma -> 1.40936}} *)

Did the result still come from the complied version or NOT???

A work around is

new2["LNB"][data_, kdata_, f0_?NumericQ, mu_?NumericQ, sigma_?NumericQ] := mytest["LNB"][data, kdata, f0, mu, sigma];
NMaximize[{new2["LNB"][B1997, kB1997, f0, mu, sigma], f0 > 0, sigma > 0}, {f0, mu, sigma}, Method -> {"RandomSearch", "SearchPoints" -> 20}]

This takes really really really long to run!!

Can this be improved?