# [R] Weibull maximum likelihood estimates for censored data

Date: Tue, 08 Apr 2008 12:14:45 +0200

Hello!

The problem is that if I type out the code as below the likelihood ratio is greater than one.

> Interest

D C

1 17 1

2 10 0

3 15 0

4 2 0

5 42 0

6 53 0

7 193 0

8 11 0

9 2 0

10 8 0

11 12 1

library(stats4)

dur_ind_test = function (CDMatrix) # Matrix with durations and censores

{

lLnw <- function(b){

D = CDMatrix

NT = nrow(D)

a =((NT-D[1,2]-D[NT,2])/ sum(D[,1]^b))^(1/b)

f = sum(log((a^b)*b*(D[2:(NT-1),1]^(b-1))*exp(-((a*D[2:(NT-1),1])^b))))

fd1 = (a^b)*b*(D[1,1]^(b-1))*exp(-(a*D[1,1])^b)

fdn = (a^b)*b*(D[NT,1]^(b-1))*exp(-(a*D[NT,1])^b)

S1 = exp(-(a*D[1,1])^b)

SN = exp(-(a*D[NT,1])^b)

-(D[1,2]*log(S1)+(1-D[1,2])*log(fd1)+ f + D[NT,2]*log(SN)+ (1-D[NT,2])*log(fdn))

}

lLne <- function(A){

D = CDMatrix

NT = nrow(D)

b=1

f = sum(log(A*b*(D[2:(NT-1),1]^(b-1))*exp(-(A^(1/b)*D[2:(NT-1),1])^b)))

fd1 = A*b*(D[1,1]^(b-1))*exp(-(A^(1/b)*D[1,1])^b)

fdn = A*b*(D[NT,1]^(b-1))*exp(-(A^(1/b)*D[NT,1])^b)

S1 = exp(-(A^(1/b)*D[1,1])^b)

SN = exp(-(A^(1/b)*D[NT,1])^b)

lLw = D[1,2]*log(S1)+(1-D[1,2])*log(fd1)+ f + D[NT,2]*log(SN)+ (1-D[NT,2])*log(fdn)

-(D[1,2]*log(S1)+(1-D[1,2])*log(fd1)+ f + D[NT,2]*log(SN)+ (1-D[NT,2])*log(fdn))

}

fit1 <- mle(lLnw,start = list(b = 0.5)) # Estimate parameters using ml

fit2 <- mle(lLne,start = list(A = 0.02))

```Lw <- lLnw(coef(fit1))                      # Maximum log likelihood :
Weibull

Le <- lLne(coef(fit2))                      # Maximum log likelihood :
```
Exponential
```LR0 <- (Le/Lw)                        # Likelihood ratio with duration
```
sample

NSimM <- cbind(as.matrix(sort(rchisq(nsim,1,0))),runif(nsim,0,1)) # chi-square df1 simulations, uniform rvs

Uniftest <- runif(1,0,1)

```firstrow <- cbind(LR0,Uniftest)                                       #
```
use sample LR as LR

NSimM <- rbind(firstrow,NSimM)

Test <- matrix(rep(0,2*(nsim+1)),nrow=(nsim+1))

NSimM <- cbind(NSimM,Test)

```for(i in 2:nsim+1) {                                                 #
```
indicates the number of simulation above the sample
```                                       if (NSimM[i,1]< LR0)NSimM[i,3]<-
1                          # likelihood ratio

else if(NSimM[i,1]== LR0)if(NSimM[i,2]>= Uniftest)NSimM[i,4]<-1
```
# if equal, only indicate if rv for simulation
```                   }
```

# is larger that rv for sample LR

Gn <- 1-(sum(NSimM[,3]))/nsim + sum(NSimM[,4])/nsim

pval <- (nsim*Gn+1)/(nsim+1)
#Calculate Monte Carlo p-value

out <- c(pval,confint(fit1))

now <- c(Le,Lw)

LR0

}}

> test_1 <- dur_ind_test(CDMatrix = Interest,nsim=1000)

Profiling...

> test_1

A b

42.32406 41.59035

=> likelihood ratio = 1.017641

The disclaimer also provides our corporate information and names of our directors as required by law.

The disclaimer is deemed to form part of this message in terms of Section 11 of the Electronic Communications and Transactions Act 25 of 2002. If you cannot access the disclaimer, please obtain a copy thereof from us by sending an email to: disclaimer_at_investec.co.za

[[alternative HTML version deleted]]

R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Tue 08 Apr 2008 - 10:58:02 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Wed 09 Apr 2008 - 14:30:28 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.