Message-ID: <29988342.1075840441159.JavaMail.evans@thyme>
Date: Thu, 16 Aug 2001 09:51:00 -0700 (PDT)
From: s-news-owner@lists.biostat.wustl.edu
To: dmck@u.washington.edu
Subject: Re: [S] confidence intervals for Weibull shape parameter
Cc: s-news@wubios.wustl.edu
Mime-Version: 1.0
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit
Bcc: s-news@wubios.wustl.edu
X-From: s-news-owner@lists.biostat.wustl.edu
X-To: "D. Mckenzie" <dmck@u.washington.edu>
X-cc: s-news@wubios.wustl.edu
X-bcc: 
X-Folder: \mark fischer 7-12-02\Notes Folders\All documents
X-Origin: FISCHER-M
X-FileName: mark fischer 7-12-02 LNotes.nsf

"D. Mckenzie" <dmck@u.washington.edu> writes:

> Splus 2000 r3.
>
> I have fit a 2-parameter Weibull distribution to a data vector (fire-free
> intervals in forests) and would like to compute the probability that the
> shape parameter is greater than 1.  I can use bootstrap() [I presume] to
> compute the standard error.  Is there a closed expression for CIs? (I
> imagine a normal approximation would be pretty lame here).
>
> Thanks in advance for any suggestions.

Profile likelihood would work well here.  Because the conditional mle
for the scale parameter, given the shape parameter, has a closed-form
expression you can evaluate the profile log-likelihood for the shape
parameter explicitly.

I used this as an example in a presentation at the recent Joint
Statistical Meetings on "Using Open Source Software to Teach
Mathematical Statistics".  (As you might guess from the title, I was
primarily discussing another implementation of the S language but the
examples work - with minor modifications - in S-PLUS.)  The PDF file for
the presentation is available as
http://www.stat.wisc.edu/~bates/JSM2001.pdf
Look especially at slide 26 where the profile log-likelihood function
for the shape parameter is defined as
> profilell <- function(alpha)
+    -sum(dweibull(xmp04.30$lifetime, shape = alpha,
+         scale = mean(xmp04.30$lifetimealpha)(1/alpha), log = TRUE))
> profilell(2.15)   # negative profile log-likelihood at estimate
[1] 77.0951
> mle3 <- nlm(profilell, c(alpha = 1.0), hessian = TRUE)
> unlist(mle3[-3])
minimum   estimate    hessian       code iterations
77.095088   2.152001   3.378437   1.000000   7.000000
(This is R, not S-PLUS.  In S-PLUS I think you must evaluate dweibull
without log = TRUE then take the log of the result and you would use
ms rather than nlm.)

Using the approximate distribution of twice the change in the profile
log-likelihood as a chi-square with 1 df (at least I think that is the
correct formula), you can construct confidence limits.
---------------------------------------------------------------------
This message was distributed by s-news@lists.biostat.wustl.edu.  To
unsubscribe send e-mail to s-news-request@lists.biostat.wustl.edu with
the BODY of the message:  unsubscribe s-news