График истинных и расчетных функций уровня опасности - PullRequest
Новая
       93

График истинных и расчетных функций уровня опасности

0 голосов
/ 29 января 2019

Я хочу построить функции истинной и предполагаемой степени опасности, используя следующий код.Но его перестали беспокоить.Я хочу, чтобы такой сюжет с двумя строками показывал оба.

Код 

require(VGAM)
n<-200
k=400
y<-rexp(n,1)

mu<-mean(y)
st<-sd(y)

 h<-0.79 * IQR(y) * length(y) ^ (-1/5)
 x <- seq(min(y) + 0.05, max(y), length =k)

 KLN <- matrix(rep(0, k * n), ncol = k)
 PhiLN<- matrix(rep(0, k * n), ncol = k)

 fhat <- rep(0, k)
 Fhat <- rep(0, k)
 Hazardhat <- rep(0, k)
######################
for(j in 1:k) {
for(i in 1:n) {

KLN[i, j] <- (1/(y[i]*(sqrt(8*pi*log(1+h)))))*exp(-(((log(y[i])-log(x[j]))^2)/(8*log(1+h))))
PhiLN[i, j] <--(sqrt(pi)*exp((sqrt(log(h+1)))/(sqrt(2)))*erfc((-log(x[j])+log(y[i])+sqrt(2+log(h+1))))/(2^(7/4)*(log(h+1))^(1/4)))/(2^(7/4)*(log(h+1))^(1/4))

}
fhat[j] <- 1/n * (sum(KLN[, j]))
Fhat[j] <- 1/n * (sum(PhiLN[, j]))
Hazardhat[j] <- fhat[j]/(1 - (-1*Fhat[j]))
}
d<-dlnorm(y, mu,st)
cdf<-plnorm(y, mu,st)

hz<-d/(1-cdf)
#plot(x,fhat, type = "s", ylab = "Density Function", lty = 1, xlab = "Time")
plot(x, Hazardhat, type = "l",ylab = "Hazard Rate Function", lty = 1, xlab = "Time")
lines(y, hz,type="l",col="red")
legend("topright", c("Real Hazard Density", "Hazard by Lognormal Kernel"),
       col=c("red", "black"), lty=c(1,2))

Еще одна вещь, которую я хочу знать, правильна моя техника или нет (то, как сравниватьфункция опасности с истинным)?

...