![enter image description here](https://i.stack.imgur.com/0NAXF.png)

```
pj[mu_, sigma_] := 
  Table[
    NIntegrate[
        (1/(1 + Exp[-x]))^j
        *
        (1 - 1/(1 + Exp[-x]))^(10 - j)
        /sigma/Sqrt[2*Pi]
        *
        Exp[-(x-mu)^2/2/sigma^2], 
        {x, -Infinity, Infinity}
    ],
    {j, 0, 10}
]

pj[-1.9575, 0.3432]
pj[-4, 9]

```

I need this expression for a global optimization. But for some values, `NIntegrate` fails. I have tried a couple of options, such as `MaxPoints` and `MaxRecursion`, with no luck.

During the optimization, I still get all kinds of warnings and, sometimes, it just failed.

Below is a copy of the full code, adjusted to helpful comments and suggestion. But I am finding it hard to modify the `pj` term to be able to evaluate quickly (hopefully), without any warnings.

```
ClearAll["Global`*"]

link1 = {84, 54, 36, 21};
klink1 = 4;

skinks = {56, 19, 28, 18, 24, 14, 9};
kskinks = 7;

taxicabsA:={142, 81, 49, 7, 3, 1};
ktaxicabsA:= 10;

B1999:={11,12,10,4,4,1,4,2,3,3,0,2,4,1,1,0,1,1,1,2,0,1,0,0,0,0,0,0,0,1,3};
kB1999:= 50;

(* Modified Version *)
mylikLNB2[data_, kdata_, f0_, mu_, sigma_] := Module[
    {K, pj, fj, j, N0, loglik, above, below},

    K = kdata;

    pj[muin_, sigmain_] :=

    Table[
        With[
            {integrand=Simplify[
                (1/(1 + Exp[-x]))^j 
                * 
                (1 - 1/(1 + Exp[-x]))^(K - j)/sigmain/Sqrt[2*Pi]
                *
                Exp[-(x - muin)^2/2/sigmain^2]
            ]},
            NIntegrate[integrand, {x, -Infinity, Infinity}]
        ]
        ,{j, 0, K}
    ];

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

    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

]

(* Could try other datasets *)
testfun[f0_, mu_, sigma_] := 
mylikLNB2[taxicabsA, ktaxicabsA, f0, mu, sigma];

(* Desired Result *)
testfun[106.1337,-1.9575,0.3432]

ans=NMaximize[
    {testfun[f0, mu, sigma], 
    f0 >= 0 && sigma>0 },
    {f0, mu, sigma},
    Method -> {"RandomSearch", "SearchPoints" -> 50}
]

(* Gets the Desired Result, but not quite 'fast' enough *)

```

Anyone has a solution?