eppminvmgf <- function(lambda,intlim=100) {
#  Computes EPPM distribution by inverting moment generating function
#  Gordon Smyth, U of Queensland, gks@maths.uq.edu.au
#  3 Nov 1999

# Count lambda
n <- length(lambda)-1
if(n == 0) return(-lambda)
lambdan <- lambda[1:n]
if(any(lambdan <= 0)) return(-Inf)
if(any(is.na(lambda))) return(NA)
lambda1 <- lambda[n+1]
if(lambda1 < 0) lambda1 <- 0
l0 <- min(lambda)

#  Gauss quadrature nodes
z <- c(
3.922075094838088e-003,
1.176598395614269e-002,
1.960916884718173e-002,
2.745114717060762e-002,
3.529143640331391e-002,
4.312955412612517e-002,
5.096501805348071e-002,
5.879734606310993e-002,
6.662605622569780e-002,
7.445066683453822e-002,
8.227069643517379e-002,
9.008566385502030e-002,
9.789508823297344e-002,
1.056984890489966e-001,
1.134953861536878e-001,
1.212852997978233e-001,
1.290677506618775e-001,
1.368422598855154e-001,
1.446083490970573e-001,
1.523655404429137e-001,
1.601133566169875e-001,
1.678513208900435e-001,
1.755789571390419e-001,
1.832957898764339e-001,
1.910013442794196e-001,
1.986951462191640e-001,
2.063767222899702e-001,
2.140455998384087e-001,
2.217013069924000e-001,
2.293433726902496e-001,
2.369713267096325e-001,
2.445846996965261e-001,
2.521830231940906e-001,
2.597658296714929e-001,
2.673326525526744e-001,
2.748830262450599e-001,
2.824164861682056e-001,
2.899325687823854e-001,
2.974308116171127e-001,
3.049107532995963e-001,
3.123719335831298e-001,
3.198138933754101e-001,
3.272361747667861e-001,
3.346383210584345e-001,
3.420198767904603e-001,
3.493803877699218e-001,
3.567194010987776e-001,
3.640364652017540e-001,
3.713311298541299e-001,
3.786029462094405e-001,
3.858514668270945e-001,
3.930762456999058e-001,
4.002768382815364e-001,
4.074528015138502e-001,
4.146036938541740e-001,
4.217290753024665e-001,
4.288285074283917e-001,
4.359015533982959e-001,
4.429477780020864e-001,
4.499667476800100e-001,
4.569580305493311e-001,
4.639211964309046e-001,
4.708558168756464e-001,
4.777614651908954e-001,
4.846377164666684e-001,
4.914841476018050e-001,
4.983003373300022e-001,
5.050858662457342e-001,
5.118403168300588e-001,
5.185632734763087e-001,
5.252543225156631e-001,
5.319130522426016e-001,
5.385390529402357e-001,
5.451319169055213e-001,
5.516912384743425e-001,
5.582166140464739e-001,
5.647076421104140e-001,
5.711639232680903e-001,
5.775850602594351e-001,
5.839706579868289e-001,
5.903203235394111e-001,
5.966336662172561e-001,
6.029102975554135e-001,
6.091498313478103e-001,
6.153518836710152e-001,
6.215160729078608e-001,
6.276420197709256e-001,
6.337293473258711e-001,
6.397776810146355e-001,
6.457866486784805e-001,
6.517558805808900e-001,
6.576850094303208e-001,
6.635736704028021e-001,
6.694215011643836e-001,
6.752281418934298e-001,
6.809932353027607e-001,
6.867164266616352e-001,
6.923973638175789e-001,
6.980356972180514e-001,
7.036310799319550e-001,
7.091831676709819e-001,
7.146916188107980e-001,
7.201560944120635e-001,
7.255762582412884e-001,
7.309517767915205e-001,
7.362823193028674e-001,
7.415675577828472e-001,
7.468071670265715e-001,
7.520008246367540e-001,
7.571482110435492e-001,
7.622490095242150e-001,
7.673029062226012e-001,
7.723095901684608e-001,
7.772687532965855e-001,
7.821800904657598e-001,
7.870432994775373e-001,
7.918580810948349e-001,
7.966241390603457e-001,
8.013411801147671e-001,
8.060089140148458e-001,
8.106270535512367e-001,
8.151953145661747e-001,
8.197134159709602e-001,
8.241810797632528e-001,
8.285980310441792e-001,
8.329639980352463e-001,
8.372787120950650e-001,
8.415419077358792e-001,
8.457533226399019e-001,
8.499126976754554e-001,
8.540197769129165e-001,
8.580743076404638e-001,
8.620760403796267e-001,
8.660247289006368e-001,
8.699201302375783e-001,
8.737620047033380e-001,
8.775501159043530e-001,
8.812842307551568e-001,
8.849641194927209e-001,
8.885895556905923e-001,
8.921603162728260e-001,
8.956761815277105e-001,
8.991369351212875e-001,
9.025423641106625e-001,
9.058922589571079e-001,
9.091864135389557e-001,
9.124246251642806e-001,
9.156066945833716e-001,
9.187324260009921e-001,
9.218016270884277e-001,
9.248141089953197e-001,
9.277696863612855e-001,
9.306681773273248e-001,
9.335094035470081e-001,
9.362931901974517e-001,
9.390193659900746e-001,
9.416877631811375e-001,
9.442982175820649e-001,
9.468505685695475e-001,
9.493446590954261e-001,
9.517803356963541e-001,
9.541574485032415e-001,
9.564758512504756e-001,
9.587354012849222e-001,
9.609359595747021e-001,
9.630773907177476e-001,
9.651595629501329e-001,
9.671823481541837e-001,
9.691456218663593e-001,
9.710492632849130e-001,
9.728931552773253e-001,
9.746771843875128e-001,
9.764012408428104e-001,
9.780652185607274e-001,
9.796690151554774e-001,
9.812125319442815e-001,
9.826956739534435e-001,
9.841183499241993e-001,
9.854804723183394e-001,
9.867819573236031e-001,
9.880227248588493e-001,
9.892026985790009e-001,
9.903218058797662e-001,
9.913799779021416e-001,
9.923771495366982e-001,
9.933132594276648e-001,
9.941882499768158e-001,
9.950020673471962e-001,
9.957546614667240e-001,
9.964459860317572e-001,
9.970759985107923e-001,
9.976446601486404e-001,
9.981519359718442e-001,
9.985977947971757e-001,
9.989822092480798e-001,
9.993051557937031e-001,
9.995666148626408e-001,
9.997665712667109e-001,
9.999050164693313e-001,
9.999819727039625e-001)

# Gauss quadrature weights
w <- c(
7.844109968033452e-003,
7.843627313765159e-003,
7.842662034926711e-003,
7.841214190912478e-003,
7.839283870809467e-003,
7.836871193391699e-003,
7.833976307113131e-003,
7.830599390098595e-003,
7.826740650132155e-003,
7.822400324645335e-003,
7.817578680701681e-003,
7.812276014980861e-003,
7.806492653759896e-003,
7.800228952893708e-003,
7.793485297792679e-003,
7.786262103399253e-003,
7.778559814162230e-003,
7.770378904009650e-003,
7.761719876319248e-003,
7.752583263887805e-003,
7.742969628898372e-003,
7.732879562885471e-003,
7.722313686698756e-003,
7.711272650464982e-003,
7.699757133547775e-003,
7.687767844505989e-003,
7.675305521050034e-003,
7.662370929996410e-003,
7.648964867220783e-003,
7.635088157608744e-003,
7.620741655005201e-003,
7.605926242161670e-003,
7.590642830682278e-003,
7.574892360967354e-003,
7.558675802155705e-003,
7.541994152064951e-003,
7.524848437130195e-003,
7.507239712340739e-003,
7.489169061175419e-003,
7.470637595535511e-003,
7.451646455676788e-003,
7.432196810139011e-003,
7.412289855674160e-003,
7.391926817172973e-003,
7.371108947589088e-003,
7.349837527862561e-003,
7.328113866840492e-003,
7.305939301197013e-003,
7.283315195350443e-003,
7.260242941379995e-003,
7.236723958939596e-003,
7.212759695170740e-003,
7.188351624613478e-003,
7.163501249115665e-003,
7.138210097740447e-003,
7.112479726672305e-003,
7.086311719121217e-003,
7.059707685225427e-003,
7.032669261952043e-003,
7.005198112996485e-003,
6.977295928680243e-003,
6.948964425846757e-003,
6.920205347755724e-003,
6.891020463975865e-003,
6.861411570276030e-003,
6.831380488514855e-003,
6.800929066528495e-003,
6.770059178016801e-003,
6.738772722428458e-003,
6.707071624843735e-003,
6.674957835856068e-003,
6.642433331452300e-003,
6.609500112890769e-003,
6.576160206578465e-003,
6.542415663945907e-003,
6.508268561321538e-003,
6.473720999803489e-003,
6.438775105130458e-003,
6.403433027550780e-003,
6.367696941690461e-003,
6.331569046418982e-003,
6.295051564714206e-003,
6.258146743525603e-003,
6.220856853635890e-003,
6.183184189521322e-003,
6.145131069210685e-003,
6.106699834142454e-003,
6.067892849020788e-003,
6.028712501670110e-003,
5.989161202888002e-003,
5.949241386297149e-003,
5.908955508195213e-003,
5.868306047404153e-003,
5.827295505117213e-003,
5.785926404745461e-003,
5.744201291762239e-003,
5.702122733546628e-003,
5.659693319225430e-003,
5.616915659514061e-003,
5.573792386555545e-003,
5.530326153758912e-003,
5.486519635635703e-003,
5.442375527635528e-003,
5.397896545980103e-003,
5.353085427496272e-003,
5.307944929447458e-003,
5.262477829364106e-003,
5.216686924872713e-003,
5.170575033523709e-003,
5.124144992618194e-003,
5.077399659033086e-003,
5.030341909045694e-003,
4.982974638156415e-003,
4.935300760910802e-003,
4.887323210720078e-003,
4.839044939680790e-003,
4.790468918392931e-003,
4.741598135777476e-003,
4.692435598892106e-003,
4.642984332746608e-003,
4.593247380116269e-003,
4.543227801355065e-003,
4.492928674207121e-003,
4.442353093617368e-003,
4.391504171541151e-003,
4.340385036752798e-003,
4.288998834652972e-003,
4.237348727075186e-003,
4.185437892091362e-003,
4.133269523816058e-003,
4.080846832210197e-003,
4.028173042883296e-003,
3.975251396895209e-003,
3.922085150556601e-003,
3.868677575228504e-003,
3.815031957121228e-003,
3.761151597091960e-003,
3.707039810441893e-003,
3.652699926711929e-003,
3.598135289477985e-003,
3.543349256145396e-003,
3.488345197742041e-003,
3.433126498711093e-003,
3.377696556702791e-003,
3.322058782365290e-003,
3.266216599134939e-003,
3.210173443025371e-003,
3.153932762416415e-003,
3.097498017841734e-003,
3.040872681775918e-003,
2.984060238420776e-003,
2.927064183491137e-003,
2.869888023999494e-003,
2.812535278040487e-003,
2.755009474574298e-003,
2.697314153209491e-003,
2.639452863985340e-003,
2.581429167153195e-003,
2.523246632957860e-003,
2.464908841417419e-003,
2.406419382103285e-003,
2.347781853919304e-003,
2.288999864880338e-003,
2.230077031890253e-003,
2.171016980519458e-003,
2.111823344781990e-003,
2.052499766911750e-003,
1.993049897138660e-003,
1.933477393464208e-003,
1.873785921436455e-003,
1.813979153924635e-003,
1.754060770893585e-003,
1.694034459177511e-003,
1.633903912253591e-003,
1.573672830015436e-003,
1.513344918545886e-003,
1.452923889890413e-003,
1.392413461829867e-003,
1.331817357653855e-003,
1.271139305934282e-003,
1.210383040300020e-003,
1.149552299212806e-003,
1.088650825745705e-003,
1.027682367365440e-003,
9.666506757214250e-004,
9.055595064455076e-004,
8.444126189708643e-004,
7.832137763839228e-004,
7.219667453370428e-004,
6.606752960773896e-004,
5.993432027083248e-004,
5.379742439462696e-004,
4.765722050100428e-004,
4.151408823342003e-004,
3.536840961310278e-004,
2.922057280352296e-004,
2.307098554425047e-004,
1.692013733654036e-004,
1.076903810218175e-004,
4.626372417712238e-005)

# Find canonical parameter theta to make the tilted mean = 1
theta <- min(l0-1,mean(lambda)-n-1)
dif <- sum(1/(lambda-theta)) - 1
while (dif > 1e-15) {
	ddif <- sum(1/(lambda-theta)^2)
	theta <- theta - dif / ddif
	dif <- sum(1/(lambda-theta)) - 1
}

# Cumulants of tilted distribution
k2 <- sum(1/(lambda-theta)^2)
k4 <- sum(1/(lambda-theta)^4)
k6 <- sum(1/(lambda-theta)^6)

#  Estimate infinity
tol <- 20
coef <- c(tol,-k2/2,k4/24,-k6/720)
r <- polyroot(coef)
if(is.null(intlim)) intlim <- max(Re(r[abs(Im(r))<1e-14]))

# Gauss quadrature
z <- z*intlim*1i
lam <- matrix(1,200,1) %*% (lambda-theta)
x <- 1
f <- Re( exp(log(w) - x*z + log(lam/(lam-z))%*%rep.int(1,n+1) ))
plot(Im(z),f,type="l")
f <- sum(rev(f))/pi*intlim

log(f)-theta+sum(log(lambdan/(lambdan-theta)))-log(lambda1-theta)
}
