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) }