# Some problems with complied function

(* 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[
{

(* 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?

• "Still has TWO calls to MainEvaluate." - That would have something to do with PadRight`? Compilable functions – dr.blochwave May 19 '15 at 21:41