From: Nadia Theron <NadiaT_at_investec.co.za>

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

# is larger that rv for sample LR

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

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 durationsample

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

Could someone please help?

http://www.investec.com/EmailDisclaimer/emaildisclaimer.htm

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.
*