Making a few changes in your code, namely

- loglik:=f0*ln(IntegrationTools:-Combine(expand(data.(pj)/f0)));

above:=lnGAMMA(N+1):

- loglik:=unapply(above-below+loglik,f0,mu,sigma);

- pj := seq(binomial(K, j)*(evalf(Int((1-1/(1+exp(-x)))^(K-j)*exp(-(x-mu)^2/(2*sigma^2))/((1+exp(-x))^j*sigma*sqrt(2*Pi)), x = -infinity .. infinity))), j = 0 .. K);,

and using

DirectSearch[Search]('loglik(f0,mu,sigma)',{f0>=0,sigma>=0,mu<=0,mu>=-50,f0<=50,sigma<=50},maximize);

I obtain

test1.mw