### originally was ~/R/MM/NUMERICS/dpq-functions/pnt-precision-problem.R ### - - - - - - - - - - - - - - - - - - - - - - - - - - - ##==> see also ./t-nonc-tst.R ## ============ library(DPQ) stopifnot(exprs = { require(graphics) require(sfsmisc) # lseq(), p.m(), mult.fig() }) source(system.file(package="DPQ", "test-tools.R", mustWork=TRUE))# ../inst/test-tools.R ## => showProc.time(), ... list_() , loadList() , readRDS_() , save2RDS() (doExtras <- DPQ:::doExtras()) ## save directory (to read from): (sdir <- system.file("safe", package="DPQ")) if(!dev.interactive(orNone=TRUE)) pdf("pnt-precision-2.pdf") op <- options(nwarnings = 1e5) x <- 10^seq(2,12, by= .25) px <- pt(x, df= 0.9, ncp = .01, lower.tail=FALSE, log=TRUE) ## 16 warnings ## full precision was not achieved in 'pnt' plot(x,px, log="x", type="o") #- show catastrophic behavior ## R-devel 2008-11-04: slope is kept, but we have a jump still ## ## extend x --> "the flattening" (at log(P[ ]) ~= -32) happens still curve(pt(x, df= 0.9, ncp = .01, lower.tail=FALSE, log=TRUE), 100, 1e20, log="x") xs <- lseq(100, 1e6, len = 101) pxs <- pt(xs, df= 0.9, ncp = .01, lower.tail=FALSE, log=TRUE) coef(lm(pxs ~ log(xs))) ## (Intercept) log(x) ## -1.1484345 -0.8999986 ## But (central) pt() has now such problem: x <- 10^seq(2,17, by=1/4) mp <- cbind(x, pnt = pt(x, df= 0.9, ncp = .01, lower.tail=FALSE, log=TRUE), pt = pt(x, df= 0.9, lower.tail=FALSE, log=TRUE)) mp p.m(mp, type="l", log="x", xlab="x", main = "pt(x, df = 0.9, lower = F, log=TRUE)") legend("topright", c("ncp = 0.01", "central"), lty=1:2, col=1:2, bty="n") ## even closer: (ncp ~= 0) : --> even bigger problem : pnt(..) becomes -Inf, ## (even in R-devel 2008-11-04) m2 <- cbind(x, pnt = pt(x, df=.9, ncp=1e-4, lower.tail=FALSE, log=TRUE), pnt = pt(x, df=.9, ncp=1e-10, lower.tail=FALSE, log=TRUE), pt = pt(x, df=.9, lower.tail=FALSE, log=TRUE)) m2 p.m(m2, type="l", log="x", xlab="x", main = "pt(x, df = 0.9, lower = F, log=TRUE)") legend("topright", c("ncp = 1e-4", "ncp = 1e-10", "central"), lty=1:3, col=1:3, bty="n") ### Note: The non-central F ( <==> non-central beta ) is sometimes *WORSE*: p.t.vs.F <- function(df, ncp, x1 = 1, x2= 1e5, nout = 201, col = c("blue","red")) { ## Purpose: Tail - Comparison of non-central t & non-central F ## Abramowitz & Stegun 26.6.19, p.947 ###>> This "non-sense" non-central F <==> "two-tail" non-central t ###>> i.e., there's NO direct representation *unless* ncp = 0 ## ---------------------------------------------------------------------- ## Arguments: ## ---------------------------------------------------------------------- ## Author: Martin Maechler, Date: 1 Nov 2008, 17:12 x <- lseq(x1, x2, length=nout) m <- cbind(pt=pt(x, df=df, ncp= ncp , log=TRUE,lower=FALSE), pf=pf(x^2, df1=1, df2=df, ncp= ncp^2, log=TRUE,lower=FALSE)) matplot(x, m, type="l", log="x", ylab = "log(1 - F*(x))", lwd = 2, col=col, main = "upper tail prob. p*(......, log=TRUE, lower.tail=FALSE)") legend("topright", c(sprintf("pt(x, df = %g, ncp=%g )", df,ncp), sprintf("pf(x^2, df1=1, df2= %g, ncp=%g ^2)",df,ncp)), lty=1:2, col=col, lwd = 2, inset=.01) mtext(R.version.string, side=1, line=-1, adj=0.01, cex = 0.6) invisible(cbind(x=x, m)) } p.t.vs.F(3.14, 2) p.t.vs.F( 2, 20,, 1e9) p.t.vs.F(0.2, 20,, 1e50)# pnt()_R-devel has small jump; 2.8.0 big jumps to horizontal ## but here (df << ncp), pf() {ie pnbeta} is clearly better}: p.t.vs.F(2, 200, 10,1e6) p.t.vs.F(2, 200, 10,1e9)# breaks down eventually as well ## Small ncp, *however* one of the two is WAY wrong ! -- "SYSTEMATIC ERROR" : p.t.vs.F(3, 0.1) ## {and pf() breaks earlier for large x} p.t.vs.F(3, 0.5) ## still {slightly} p.t.vs.F(3, 1) ## no visible problem ## but here, for x -> 0: pf() is systematically larger : p.t.vs.F(3, 1, 1e-2,100) showProc.time() ### --------------------------------------------------------------------------------- ## -> pntR() is R-level pnt() function pt (1e10, df=.9, ncp=1e-10, lower.tail=FALSE, log=TRUE) # >> -Inf pntR(1e10, df=.9, ncp=1e-10, lower.tail=FALSE, log=TRUE) # " ## "finis: precision NOT reached" ## back to less extreme ncp: x <- 2^(0:40) px <- pt(x, df = 1., ncp = 1, lower.tail = FALSE, log = TRUE)## 14 warnings plot(x,px, log="x", type="o") #- show "jump" {just small kink: R-devel 25.10.09} ## also the left tail (!) --- this is now (R-devel 25.Oct.09) worse px <- pt(-x, df = 1.2, ncp = 1, log = TRUE) plot(x, px, log="x", type="o") #- jump and wrong ## larger df: no jump, still wrong px <- pt(x, df = 2, ncp = 1, lower.tail = FALSE, log = TRUE) ## 23 warnings plot(x, px, log="x", type="o") #- show bad behavior ## does my one approximation help here ? px. <- pntJW39(x, df = 2, ncp = 1, lower.tail = FALSE, log = TRUE) lines(x,px., col = 2, lwd=2) ## no! It's worse for large x! px <- pt(x, df = 2, ncp = 2, lower.tail = FALSE, log = TRUE) ## 23 warnings plot(x,px, log="x", type="o") #- show bad behavior ## The Mathematica example ## LogLogPlot[ 1 - CDF[NoncentralStudentTDistribution[3, 5], x], {x, 1, 10^6}] ## takes many minutes of computing !!!! px <- pt(x, df = 3, ncp = 5, lower.tail = FALSE, log = TRUE) plot(x, px, log="x", type="o") #- show bad behavior ## equivalent {different x-labels}: plot(log(x), px, type="o") ## BTW: Very slow convergence here {not in outer tail} -- can clearly be improved pntR(20, df=0.9, ncp=30, lower.tail=FALSE, log=TRUE, errmax=1e-15) pntR(20, df= 2, ncp=30, lower.tail=FALSE, log=TRUE, errmax=1e-15) ### "Solve this problem": Find tail formula empirically ### ---------------------------------------- ### log(1 - P(x)) ~= alpha - beta * log(x) ### ---------------------------------------- ## Empirically: alpha = g(ncp) ; beta = h(df) << see below: beta == df ( = nu) !! summary(lm(px ~ log(x), subset=x > 7 & x < 35000))# seems clear, R^2 = 0.9999 showProc.time() ptRTailAsymp <- function(df,ncp, f.x1 = 1e5, f.x0 = 1, nx = 1000, do.plot=FALSE, verbose = do.plot, ask = do.plot) { ## Purpose: Right tail asymptotic of pt_noncentral() ## ---------------------------------------------------------------------- ## Arguments: ## ---------------------------------------------------------------------- ## Author: Martin Maechler, Date: 23 Oct 2008, 15:15 stopifnot(length(df) == 1, is.numeric(df), df > 0) stopifnot(length(ncp) == 1, is.numeric(ncp), ncp >= 0) n.x0 <- round(.25 * nx) # at this index, we anchor the robust line ## Determine an x-range, where pt(*, lower.tail=FALSE, log=TRUE) ## is nearly linear alp <- .999 x0. <- qt(alp, df=df, ncp=ncp)##<< it's absurd to use this here, ## -- but unfortunately, qtAppr() is too unreliable here it <- 1 while(!is.finite(x0.) && it < 20) { it <- it+1 alp <- 1 - 2*(1-alp) x0. <- qt(alp, df=df, ncp=ncp) } stopifnot(is.finite(x0.)) x0 <- x0.* f.x0 x1 <- x0 * f.x1 it <- 1; converged <- FALSE if(do.plot && ask) { op <- par(ask=TRUE); on.exit(par(op)) } repeat { xx <- lseq(x0, x1, len = nx) lx <- log(xx) px <- pt(xx, df=df, ncp=ncp, lower.tail=FALSE, log=TRUE) rob.slope <- median(diff(px), na.rm=TRUE) /mean(diff(lx)) Out <- !is.finite(px) SSY <- sum((px[!Out] - mean(px[!Out]))^2) if(do.plot) { r <- lm.fit(cbind(1, lx[!Out]), px[!Out]) if(it == 1) { ## first plot plot(px ~ lx, xlab = quote(log(x)), xlim = extendrange(r=log(c(x0,x1)), f=0.2), ylim = extendrange(px, f = 0.25), main = sprintf("pt(exp(.), df=%g, ncp=%g, lower.tail=FALSE, log.p=TRUE)", df, ncp)) mtext(sprintf("Search for [x0, x1] s.t. pt(*, log=TRUE) is *linear* in it; starting @ [%g, %g]", log(x0), log(x1))) } else { # it >= 2 ... subsequent iterations lines(px ~ lx) lines(lx[!Out], r$fitted, col= "green2") abline(a = px[n.x0] - rob.slope*lx[n.x0], b = rob.slope, col = "tomato", lwd=2) } } rob.fitted <- px[n.x0] + rob.slope*(lx - lx[n.x0]) rob.resid <- px - rob.fitted R2.rob <- 1 - sum(rob.resid^2) / SSY cat(sprintf("robust R^2 = %12.8f\n", R2.rob)) if(R2.rob < 0.99) { x1 <- x1/10 if(chg0 <- x0 > x1/2) x0 <- x0/2 message(c(sprintf("'robust' fit not ok\n ==> x1 := x1 / 10 = %g", x1), if(chg0)sprintf(" -> x0 changed too, to x0= %g", x0))) } else { ## reasonable R2.rob if(all(sml.res <- abs(rob.resid) < 4*mad(rob.resid))) { message("converged in ", it, " iterations") converged <- TRUE break } ## else if (!sml.res[nx]) { ## largest x[] does NOT have small resid x1 <- exp(max(lx[sml.res])) if(verbose) cat(sprintf("new x1: log(x1)= %12g\n", log(x1))) } else ## throw away small x *ONLY* after large x are fixed if (!sml.res[1]) { ## smallest x[] does NOT have small resid x0 <- exp(min(lx[sml.res])) if(verbose) cat(sprintf("new x0: log(x0)= %12g\n", log(x0))) } else ## no change at both ends.. hmm if(R2.rob > .99999 || it > 10) { message("** R^2-pseudo-converged in ", it, " iterations") converged <- NA break } } if((it <- it+1) > 1000) { if(verbose) warning("too many iterations!") break } }## end{repeat} r <- lm.fit(cbind(1, lx), px) R2 <- 1 - sum(r$resid^2) / SSY cf <- r$coefficients if(do.plot) { legend("topright", c(sprintf("final int. [%g, %g]", log(x0), log(x1)), sprintf("coef. = (%g, %g)", cf[1], cf[2])), pch=NA, bty="n") abline(v = log(c(x0,x1)), col=adjustcolor(1, 1/2), lty=2) } ## browser() c(cf, df=df, ncp=ncp, converged=converged, x0=x0, x1=x1, alp=alp, x0.999=x0., R2 = R2, R2.rob=R2.rob) } p.tailAsymp <- function(r, add.central=FALSE, F.add = FALSE) { x01 <- r[c("x0","x1")] stopifnot(length(x01) == 2, is.numeric(x01)) curve(pt(x, df=r["df"], ncp=r["ncp"], log=TRUE,lower.tail=FALSE), x01[1] / 100, x01[2]*100, log="x", ylab = "", main = sprintf(paste("pt(x, df=%g, ncp=%g, log", if(prod(par("mfrow")> 2)) "..)" else "=TRUE, lower.tail=FALSE)", sep=''), r["df"], r["ncp"])) if(F.add) { ## This "non-sense" non-central F <==> "two-tail" non-central t ## i.e., there's NO direct representation *unless* ncp = 0 colF <- "forestgreen" curve(pf(x^2, df1=1, df2=r["df"], ncp= r["ncp"]^2, log=TRUE,lower.tail=FALSE), add = TRUE, col = colF) x. <- 10^par("usr")[2] text(x., pf(x.^2,df1=1, df2=r["df"],ncp= r["ncp"]^2, log=TRUE,lower.tail=FALSE), "pf(x^2, (1,df), ncp^2, ...)", adj=c(1,-.2), col = colF) } if(add.central) { curve(pt(x, df=r["df"], log=TRUE,lower.tail=FALSE), add=TRUE, col = "midnightblue", lwd=2) } abline(v = x01, col ="gray") cL <- "tomato" mtext(sprintf("slope = %-10.6g", r[2]), line=-1 , adj=1, col=cL, cex=par("cex")) mtext(sprintf("intCpt= %-10.6g", r[1]), line=-2.2, adj=1, col=cL, cex=par("cex")) abline(a=r[1], b = r[2] * log(10), col = cL) } if(FALSE) debug(ptRTailAsymp) r35 <- ptRTailAsymp(df=3, ncp=5) if(interactive()) x11(type = "Xlib") r45 <- ptRTailAsymp(df=4, ncp=5, do.plot=TRUE) r46 <- ptRTailAsymp(df=4, ncp=6, do.plot=TRUE) r410 <- ptRTailAsymp(df=4, ncp=10, do.plot=TRUE) r420 <- ptRTailAsymp(df=4, ncp=20) p.tailAsymp(r420) p.tailAsymp(r410) p.tailAsymp(r46) p.tailAsymp(r45) r440 <- ptRTailAsymp(df=4, ncp=40, do.plot=TRUE)## "wrong" p.tailAsymp(r440) r.810 <- ptRTailAsymp(df=.8, ncp=10, do.plot=TRUE) p.tailAsymp(r.810)# tip top r.518 <- ptRTailAsymp(df=.5, ncp=18, do.plot=TRUE) p.tailAsymp(r.518) ## better (much faster convergence; better range r.518 <- ptRTailAsymp(df=.5, ncp=18, do.plot=TRUE, f.x0=1/10,f.x1 = 100) p.tailAsymp(r.518) r.110 <- ptRTailAsymp(df=.1, ncp=10, do.plot=TRUE, f.x0=1e-6,f.x1 = 1e10) p.tailAsymp(r.110) # easy simple r71 <- ptRTailAsymp(df=7, ncp=1, do.plot=TRUE) p.tailAsymp(r71)# looks good, even though slope = -6.956 ## --converging to *central* t-distrib: r7.01 <- ptRTailAsymp(df=7, ncp=.01, do.plot=TRUE) p.tailAsymp(r7.01)# looks ok r7.001 <- ptRTailAsymp(df=7, ncp=.001, do.plot=TRUE) p.tailAsymp(r7.001, add.central=TRUE)# --- but the slope is not quite -7 ... ## hmm: tt <- lseq(10,100, length=1000) px <- pt(tt, df=7, log=TRUE,lower.tail=FALSE) plot(px ~ log(tt), pch=".", main = "pt(tt, df=7, log=TRUE, lower.t=F)") ll <- lm(px ~ log(tt)) summary(ll) ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 4.595546 0.004271 1076 <2e-16 *** ## log(tt) -6.930061 0.001214 -5707 <2e-16 *** ## ^^^^^^^^^ ....... ##--- significantly different from -7 tt <- lseq(100,1e7, length=1000) # much larger range: t -> oo px <- pt(tt, df=7, log=TRUE,lower.tail=FALSE) plot(px ~ log(tt), pch=".", main = "pt(tt, df=7, log=TRUE, lower.t=F) -- larger range(tt)") ll <- lm(px ~ log(tt)) summary(ll)##--> clearly: slope == -7 -- ok ## Left tail tt <- rev(- lseq(100,1e7, length=1000)) px <- pt(tt, df=4.5, log=TRUE) plot(px ~ log(-tt), pch=".", main = "pt(tt, df=7, log=TRUE, lower.tail=TRUE)") ll <- lm(px ~ log(-tt)) summary(ll)##--> clearly: slope == -4.5 -- ok ### ====> (empirical) Theorem: slope == -df ### ============ ### and same is true for left tail ---- see below (rmat0 <- t(sapply(ls(patt="^r.*[0-9]$"), get))) rmat <- rmat0[ rmat0[,"R2.rob"] > 0.9999 ,] rmat[,2:3] p.tailAsymp(r2.30 <- ptRTailAsymp(df= 2, ncp=30)) # very nice ## "Large" ncp : (the search and resulting slope seem *wrong* ..) p.tailAsymp(r2.50 <- ptRTailAsymp(df= 2, ncp=50, do.plot=TRUE)) showProc.time() ## Now do a larger set: nd <- length(df. <- c(.1, .8, 1, 1.5, 2, 4, 10)) nc <- length(nc. <- c(.01, .1, 1, 2, 5)) ## for indexing etc: c.df <- formatC(df., width=1) c.nc <- formatC(nc., width=1) c.c <- outer(c.df, c.nc, paste, sep="_") indR <- function(i.df, i.nc) paste(c.df[i.df], c.nc[i.nc], sep="_") sfil1 <- file.path(sdir, "pnt-prec-sim1.rds") if(!doExtras && file.exists(sfil1)) { Res <- readRDS_(sfil1) } else { ## do run the simulation [always if(doExtras)] : --------------- r35 <- ptRTailAsymp(df=3, ncp=5) if(names(r35)[1] == "") names(r35)[1] <- "intercpt" Res <- matrix(NA_real_, nrow = length(r35), ## <-- length of output ncol = nd*nc, dimnames=list(names(r35), c.c)) for(i.df in seq_along(df.)) { df <- df.[i.df] cat("\ndf = ", formatC(df)," : \n----------\n") for(i.nc in seq_along(nc.)) { cat(i.nc,"") ncp <- nc.[i.nc] r <- try(ptRTailAsymp(df=df, ncp=ncp)) Res[, indR(i.df,i.nc) ] <- if(inherits(r, "try-error")) NA else r } } attr(Res, "version") <- list(R.version) save2RDS(Res, file=sfil1) }##-- end{ run simulation } -------------------------------------------- dR <- as.data.frame(t(Res)) op <- mult.fig(mfrow=c(nd,nc),marP=-1)$old.par op2 <- par(cex = par("cex") * 3/4) # or even more for(i.df in seq_along(df.)) for(i.nc in seq_along(nc.)) { r <- Res[, indR(i.df,i.nc) ] if(is.na(r["df"])) { frame() } else p.tailAsymp(r) # ------------ } par(op); par(op2) showProc.time() ### ====> Theorem: slope == -df ### ============ [now have prove for ncp=0: central t] ### and same is true for left tail ## but I see extra discontinuous behavior at x = 0 : ## df = 10 and various ncp : curve(pt(x, df=10, ncp=10, log=TRUE), -250,50, ylim=c(-50,0), n=10001) title("pt(x, df=10, ncp= *, log=TRUE), -250,50, ..) for various ncp") ## curve(pt(x, df=10, ncp=0.01, log=TRUE), add=TRUE, n=10001, col=2) curve(pt(x, df=10, ncp=0.0, log=TRUE), add=TRUE, n=10001, col=2) ## curve(pt(x, df=10, ncp=0.1, log=TRUE), add=TRUE, n=10001, col="red2") curve(pt(x, df=10, ncp=0.2, log=TRUE), add=TRUE, n=10001, col="red2") curve(pt(x, df=10, ncp=0.5, log=TRUE), add=TRUE, n=10001, col="red2") ## curve(pt(x, df=10, ncp=1, log=TRUE), add=TRUE, n=10001, col="blue2") curve(pt(x, df=10, ncp=2, log=TRUE), add=TRUE, n=10001, col="blue2") curve(pt(x, df=10, ncp=5, log=TRUE), add=TRUE, n=10001, col="blue2") ## same, df=5: here we see increasing problem at x=0 for ncp -> large curve(pt(x, df=5, ncp=100, log=TRUE), -5000,300, ylim=c(-50,0), n=10001) curve(pt(x, df=5, ncp=0, log=TRUE), add=TRUE, n=10001, col="blue",lty=3, lwd=3) for(n. in c(1:8) / 10) curve(pt(x, df=5, ncp=n., log=TRUE), add=TRUE, n=2001, col="blue2") for(n. in c(1:8)) curve(pt(x, df=5, ncp=n., log=TRUE), add=TRUE, n=2001, col="red2") for(n. in 10*c(1:8)) curve(pt(x, df=5, ncp=n., log=TRUE), add=TRUE, n=2001, col="green2") showProc.time() 1 ## From: Michael E Meredith ## To: Martin Maechler ## Subject: Re: [Rd] Noncentral dt() with tiny 'x' values (PR#8874) ## Date: Thu, 18 May 2006 19:54:09 +0800 ## Hello Martin, ## Thanks for the response. I see what you mean about the plots! ## One thing I should also mention is that I'm getting huge numbers of warnings: ## " full precision was not achieved in 'pnt' " ## I've been hunting for a pbm in my code, but have now found I get it ## with 'power.t.test(*, sd = NULL) # eg. (ptt <- power.t.test(20, 1, power=0.8, sd=NULL, tol = 1e-8)) stopifnot(inherits(ptt, "power.htest"), is.list(ptt), all.equal(1.0999525, ptt$sd, tol = 1e-7)) ## MM: This has been fixed for R 2.4.0, with NEWS entry ## ## o pt() with a very small (or zero) non-centrality parameter could ## give an unduly stringent warning about 'full precision was not ## achieved'. (PR#9171) ## I presume this is another related old buglet which is now appearing ## due to the improved reporting of warnings in R 2.3.0. ## Regards, Mike. ## MM: The note below is ## ## Fixed for R 2.3.1 --- with NEWS entry ## ----- -------- ## o dt(x, df, ncp= not.0) no longer gives erratic values for ## |x| < ~1e-12. (PR#8874) ## At 17:57 18/05/2006, you wrote: ## > >>>>> "MikeMer" == meredith ## > >>>>> on Thu, 18 May 2006 03:52:51 +0200 (CEST) writes: ## > ## > MikeMer> Full_Name: Mike Meredith ## > MikeMer> Version: 2.3.0 ## > MikeMer> OS: WinXP SP2 ## > MikeMer> Submission from: (NULL) (210.195.228.29) ## > ## > ## > ## > MikeMer> Using dt() with a non-centrality parameter and ## > near-zero values for 'x' results ## > MikeMer> in erratic output. Try this: tst <- c(1e-12, 1e-13, 1e-14, 1e-15, 1e-16, 1e-17, 0) dt(tst,16,1) ## > MikeMer> I get: 0.2381019 0.2385462 0.2296557 0.1851817 ## > 0.6288373 3.8163916 (!!) ## > MikeMer> 0.2382217 ## > ## >I get quite different values (several '0', BTW), which just ## >confirms the erratic nature. ## > ## >As often, plots give even a clearer picture: x <- lseq(1e-3, 1e-33, length= 301) plot(x, dt(x, df=16, ncp=1), type = "o", cex=.5, log = "x") plot(x, dt(x, df=16, ncp=0.1), type = "o", cex=.5, log = "x") plot(x, dt(x, df= 3, ncp=0.1), type = "o", cex=.5, log = "x") ## > ## > ## > MikeMer> The 0.238 values are okay, the others nonsense, and ## > MikeMer> they cause confusing spikes on plots of dt() vs 'x' ## > MikeMer> if 'x' happens to include tiny values. (Other ## > MikeMer> values of df and ncp also malfunction, but not all ## > MikeMer> give results out by an order of magnitude!) ## > ## >I think almost all do, once you start looking at plots like the above. ## > ## > MikeMer> I'm using the work-around dt(round(x,10),...), but ## > MikeMer> dt() should really take care of this itself. ## > ## >or actually rather do something more smart; the cutoff at 1e-10 ## >is quite crude. ## > ## >Note that this is not a new bug at all; but rather as old as ## >we have dt(*, ncp= .) in R. ## > ## >Thanks for reporting it! ## >Martin Maechler, ETH Zurich ## =============================================== ## Michael E Meredith ## Wildlife Conservation Society (WCS) Malaysia Program ## 7 Jalan Ridgeway, 93250 Kuching, Malaysia ## Fax: +60-82-252 799 Mobile: +60-19 888 1533 ## email: meredith@easynet.co.uk http://www.mered.org.uk ## Program website: http://www.wcsmalaysia.org