llPlot {dbd}

Plot the log likelihood surface for the data.

Description

Plot, as a perspective plot or a contour plot, the log likelihood surface for the data set from which parameters are being estimated.

Usage

llPlot(x, distr=c("db","betabinom"),ntop, zeta, size, alim = NULL, blim = NULL,
       ngrid = c(100, 100), plotType = c("persp", "contour", "none"),
       theta = -30, phi = 40, ...)  

Arguments

x

A vector of numeric data purportedly arising from a db or beta binomial distribution.

distr

Character string specifying which of the two relevant distributions (db, or beta binomial) is to be considered.

ntop

See mleDb() and ddb(). Ignored if distr is "betabinom".

zeta

See mleDb() and ddb(). Ignored if distr is "betabinom".

size

Integer scalar specifying the upper limit of the “support” of the beta binomial distribution under consideration. The support is the set of integers {0, 1, ..., size}. (The values of x may sometimes be considered to be the number of “successes” in size trials. The size argument is ignored if distr is "db".

alim

Numeric vector of length 2; the range of alpha values over which the surface is to be plotted. Defaults to c(0,10) if distr is "db" and to c(0,1) if distr is "betabinom".

blim

Numeric vector of length 2; the range of beta values over which the surface is to be plotted. Defaults to c(0,10) if distr is "db" and to c(0,100) if distr is "betabinom".

ngrid

The dimensions of the grid of paramter values at which the log likelihood is to be evaluated in order to plot the surface. Note that ngrid may be supplied as an integer scalar, in which case it is replicated to a vector of length 2.

plotType

Character string specifying the nature of the plot to be produced. If it is "none" then no plot is produced. The value returned may be plotted at a later occasion.

theta

An argument to be passed to persp(). Ignored unless plotType is "persp".

phi

An argument to be passed to persp(). Ignored unless plotType is "persp".

...

Other arguments that may be passed to persp() or to contour()

Details

This function could conceivably be useful in diagnosing problems with parameter estimation should these arise.

Value

A list with entries

x

The vector of values of the first parameter (alpha for distr="db", m for distr="betabinom") over which the surface is to be plotted. There are ngrid[1] such values, ranging from alim[1] to alim[2].

y

The vector of values of the second parameter (beta for distr="db", s for distr="betabinom") over which the surface is to be plotted. There are ngrid[2] such values, ranging from blim[1] to blim[2].

z

An ngrid[1] x ngrid[2] numeric matrix, specifying the surface. the value of z[i,j] is ll(x[i],y[j]) where ll() is the log likelihood function.

dxy

A data frame with columns named "alpha" and "beta" for distr="db" or "m" and "s" for distr="betabinom", and ngrid[1]*ngrid[2] rows. It is formed by applying expand.grid() to the x and y entries of this list.

fxy

A numeric vector of length ngrid[1]*ngrid[2]. Its ith value is the log likelihood evaluated at the ith row of dxy. Its entries are the same as the entries of z.

There is obviously considerable redundancy in the returned value.

The names x and y that are used for the first two entries of this list conform to the names of the arguments of persp() and contour.

If plotType is "persp" or "contour" the value is returned invisibly.

Author(s)

Rolf Turner [e-mail]

See Also

link{mleDb}() link{mleBb}() link{persp}() link{contour}()

Examples

X <- hmm.discnp::SydColDisc
X$y <- as.numeric(X$y)
X <- split(X,f=with(X,interaction(locn,depth)))
x <- X[[19]]$y
srf <- llPlot(x,ntop=5,zeta=FALSE,alim=c(0.5,0.7),blim=c(0.2,0.4),plotType="c")
## Not run: 
if(require(rgl)) {
   with(srf,plot3d(ab$alpha,ab$beta,fab)
# Allows dynamic rotation of the surface.
}

## End(Not run)
# Negative (!) parameters for the db distribution.
set.seed(42)
xs <- rdb(100,-1,-1,5)
fit <- mleDb(xs,5)
llPlot(xs,ntop=5,zeta=FALSE,alim=c(-4,2),blim=c(-4,2),plotType="c",
       main="log likelihood contours")
points(fit[1],fit[2],pch=20,col="red")
points(-1,-1,pch=20,col="blue")
legend("topright",pch=20,col=c("red","blue"),
       legend=c("estimate","true value"),bty="n")