# generate data
set.seed(5732)
test <- function(x)
        {(1e6*(x^11*(1-x)^6)+1e4*(x^3*(1-x)^10))+.1}
x <- (0:150)/150
y <- NULL
for (i in 1:151) {
    mu <- test(x[i])
    repeat {
        t.wk <- rweibull(1,2,mu)
        c.wk <- rweibull(1,1,2*mu)
        x.wk <- min(t.wk,c.wk)
        delta.wk <- t.wk<=c.wk
        z.wk <- rweibull(1,1,mu/2)
        if (x.wk>z.wk) break
    }
    y <- rbind(y,c(x.wk,delta.wk,z.wk))
}

# fit the model and evaluate
weibull.fit <- gssanova(y~x,family="weibull")
est <- predict(weibull.fit,data.frame(x=x),se=TRUE)

# plot the fit
plot(x,y[,1],type="n",xlab="u",ylab=expression(beta(u)==e^{mu(u)}))
for (i in 1:151) {
    lines(c(x[i],x[i]),c(y[i,1],y[i,3]),col=7)
}
points(x,y[,1],pch=c("+","o")[y[,2]+1],col=3)
lines(x,test(x),col=6)
lines(x,exp(est$fit),col=4)
lines(x,exp(est$fit+1.96*est$se),col=5)
lines(x,exp(est$fit-1.96*est$se),col=5)