风险函数将是。$\frac{f(t)}{1-F(t)}$

假设寿命是指数的，则在该链接处给出公式；它是，指数的速率参数；您从数据中估计速率参数的方式与对任何指数分布的估计方式完全相同。$\lambda$$\lambda$

每小时死亡率是一小时内危险率的积分，但由于危险是恒定的并且时间以小时为单位，所以它是微不足道的。

--

关于模型适用性的讨论：

如果寿命真的是指数级的，那么危险将是恒定的。但这在我看来并不是指数级的。

虽然可以通过如此少的数据获得风险函数的非参数估计，但我可能会从 Weibull 模型开始，因为它包含您建议的指数寿命模型作为特例，并且可以同时适应增加和减少的风险函数。

survreg(Surv(beelife)~1)
Call:
survreg(formula = Surv(beelife) ~ 1)

Coefficients:
(Intercept) 
   3.416515 

Scale= 0.7265582 

Loglik(model)= -140.5   Loglik(intercept only)= -140.5
n= 33


 plot(ecdf(beelife))
 lines((0:90),pweibull(0:90,1/.726558,exp(3.4165)),col=2)

合身看起来很合理。

此外，如果我们查看对数（比例）系数的标准误差，它表明指数（如果我理解一切正确，对应于 0 的对数比例）似乎与数据不一致：

             Value Std. Error     z         p
(Intercept)  3.417      0.133 25.66 3.50e-145
Log(scale)  -0.319      0.137 -2.32  2.01e-02

Scale= 0.727 

Weibull distribution
Loglik(model)= -140.5   Loglik(intercept only)= -140.5
Number of Newton-Raphson Iterations: 6 
n= 33 

如果我做得对，拟合的危险函数应该如下所示：$h(t) = f(t)/(1-F(t))$

我可能会考虑的另一个模型是伽玛，它也包括指数作为特例。它也有增加和减少的危险。但是，我怀疑在这种情况下它会做得更好。

该问题询问每小时死亡率与寿命指数分布之间的关系。

宇宙把时间分成很短的滴答声。为每只蜜蜂掷骰子。骰子有很多面，大部分是空白的——但有一个是黑色的。每当黑面出现时，蜜蜂就会死去。

我们将用计算机来模拟这个，但是因为懒惰，我不想像死神那样精细地切割时间。因此，与其每秒掷一次骰子，不如说，也许我会每小时掷一次骰子：但我会让它显示黑色面的可能性增加 3600 倍，以弥补掷骰子的不足。因为在一秒钟甚至一小时内死亡的机会是如此之小，所以这是一个很好的近似值。

为了利用计算机一次对数组进行大量计算的能力，我将一次模拟整个蜂巢一只蜜蜂。让骰子一次又一次地滚动，每小时一次，直到蜜蜂死亡。然后继续一只新生（转世？）蜜蜂，重复直到许多蜜蜂出生和死亡。结果将与所有蜜蜂共存一样。

这是R代码。它的输入是蜜蜂寿命的数组，bees，用于估计指数率。

set.seed(17)                          # (Allows reproducible results)
n <- 10^5                             # Number of time slices
units <- 1                            # Number of hours per time slice
rate.hat <- 1 / mean(bees)            # ML estimate of the death rate
deaths <- runif(n) < rate.hat * units # Roll Death's die
times <- which(deaths != 0)           # Note when bees are killed
lifetimes <- (diff(c(0,times))-1) * units # The time differences are the bee lifetimes

这是deaths数组的第一部分，显示了开始掷骰的结果（“*”表示黑面）：

paste(ifelse(head(deaths, 100),"*", "."), collapse="")

……*…………………………………………………………………………………………………………………… ..................* ..................* ..................

黑脸出现在骰子上的最初几次：

head(times)

12 67 82 120 134 147

因此，第一只蜜蜂在第 12 个滴答声中被砍倒（所以我将它的寿命计算为 11 个滴答声），第二个在第 67 个滴答声（一生 54 个滴答声），依此类推：这些是星星的位置上一个输出。

这些数字绘制了这些死亡（每个都是左侧的垂直斜线，将每个生命显示为斜线之间的水平间隙）和生命的分布（在右侧）。后者具有叠加的指数分布函数。这是一个很好的选择。

绘制直方图的代码是

hist(lifetimes, freq=FALSE, ylim=c(0, rate.hat), breaks=25, xlab="Hours")
curve(dexp(x, rate=rate.hat), col="Red", add=TRUE)

为什么直方图是指数的？ 因为平均而言，在某个年龄之后活着的蜜蜂数量在下一个时间片中减少了相同的分数。这意味着直方图必须呈指数下降。

顺便说一句，最大似然估计量的原因1/mean(bees)是指数概率分布函数的形式为$\kappa \exp(-\kappa t)$一辈子$t$. 要理解这一点，请回想一下 PDF 是一种密度：它为我们提供了每单位时间的概率。自从$t$是时间（以小时为单位），$\kappa$必须是每单位时间的（死亡）概率：它是比率。 它的对数等于$\log(\kappa) - \kappa t$. 因此，生命周期数据集的对数似然$(t_1, t_2, \ldots, t_n)$是

$$\log(L(\kappa)) = n \log(\kappa) - \sum_{i=1}^n \kappa\, t_i.$$

微积分向我们展示了（通过对$\kappa$并将其设置为零）

$$\hat\kappa = \frac{n}{\sum_{i=1}^n t_i}$$ 是唯一的关键点（显然是可能性最大化的地方）。这是平均寿命的倒数。 这就是每小时死亡率如何与最大似然机制相适应。

这是一种方法

获取数据的直方图

bl <- c(7.1, 2.3, 9.6, 25.8, 14.6, 12.8, 20.9, 30.0, 71.1, 36.9, 3.9, 18.2, 24.3, 55.8, 84.7,
          13.6, 10.8, 18.2, 19.0, 34.8, 27.1, 54.2, 34.8, 65.8, 33.3, 26.5, 9.5, 44.2, 35.7, 2.3,
          4.0, 25.9, 41.3)

b1h <- hist(bl, freq = FALSE, xlim = c(0, quantile(bl, 0.999)), plot=TRUE)
df <- data.frame(density=b1h$density, mids=b1h$mids)

将直方图绘制为点

plot(b1h$mids,b1h$density,cex=1,xlab='Life (hr)',ylab='Probability')

将指数分布拟合到直方图

fit <- nls(density ~ lambda*exp(-lambda*mids), start=list(lambda=.1), data=df)
lines(df$mids,predict(fit),col='red')

或者，将指数分布拟合到数据

fit1 <- fitdistr(bl, "exponential") 
curve(dexp(x, rate = fit1$estimate), col = "green", add = TRUE)

平均寿命为1/λ，死亡率为λ

cat('from nls (red) the mean life is',1/coef(fit),'and mortality rate is',coef(fit),'per year')
cat('from fitdistr (green) the mean life is',1/coef(fit1),'and mortality rate is',coef(fit1),'per year')

结果是

from nls (red) the mean life is 35.47136 and mortality rate is 0.02819176 per year
from fitdistr (green) the mean life is 27.84848 and mortality rate is 0.0359086 per year

函数 fitdisr 提供参数的最大似然估计。nls 函数提供最小二乘估计。注意，对于指数分布，最大似然估计由 n/sum(xi) 给出，参考推导