I get a problem when I run this program in R. anybody help me to solving this problem..?
par_1<-cbind(c(5.038159),c(3.899621))
par_2<-cbind(c(2.435457),c(13.89517))
tau<-365
cdf2 <- function(x, help) {
pgamma(x, shape=par_1[1], scale=par_1[2]) *
pgamma(x, shape=par_2[1], scale=par_2[2])-help
}
nextEventTime <- function(censoring) {
randomNumber <- runif(n=1, min=0, max=1)
pnew <- randomNumber * (1 - cdf2(censoring, 0)) + cdf2(censoring, 0)
uniroot(f=cdf2, interval=c(0, 1000*tau), help=pnew)$root
}
hazardRate1 <- function(t) {
dgamma(t, shape=par_1[1], scale=par_1[2]) /
(1 - pgamma(t, shape=par_1[1], scale=par_1[2]))
}
hazardRate2 <- function(t) {
dgamma(t, shape=par_2[1], scale=par_2[2]) /
(1 - pgamma(t,shape=par_2[1], scale=par_2[2]))
}
nextEventType <- function(t) {
p <- hazardRate1(t)/(hazardRate1(t)+hazardRate2(t))
randomNumber <- runif(n=1, min=0, max=1)
if (randomNumber <= p) {1} else {2}
}
baris<-c(1:20000)
nexteventtime<-rep(0,time=20000)
nexteventype<-rep(0,time=20000)
dfnexteventime<-data.frame(baris,nexteventtime,nexteventype)
for(i in 1:nrow(dfnexteventime)){
dfnexteventime$nexteventtime[i]<-nextEventTime(dfnexteventime$nexteventtime[i])
dfnexteventime$nexteventype[i]<-nextEventType(dfnexteventime$nexteventtime[i])
}
View(dfnexteventime)
When I run this program, this program will error & produce output like this
Error in if (randomNumber <= p) { : missing value where TRUE/FALSE needed
I think this problem because t value in nextEventType(t)
function can't zero (t!=0). But nextEventTime(dfnexteventime$nexteventtime[i])
never produce zero value, when I run this part for 10 times,
baris<-c(1:20000)
nexteventtime<-rep(0,time=20000)
nexteventype<-rep(0,time=20000)
dfnexteventime<-data.frame(baris,nexteventtime,nexteventype)
for(i in 1:nrow(dfnexteventime)){
dfnexteventime$nexteventtime[i]<-nextEventTime(dfnexteventime$nexteventtime[i])
}
without nextEventType function. This part never produce 0 value. So, I confuse, what is a problem?.
I want result nextEventType(t)
produce not zero value. because if using zero value will be Error in if(ramdonNumber <= p) { :...
Aucun commentaire:
Enregistrer un commentaire