4

背景

我正在尝试预测产品线的销售额（最后样本中的 y_test）。它在一段时间内的销售额基于另一个产品 (x_test) 的所有先前销售额以及这些先前销售额中有多少仍在使用中。但是，无法直接衡量那些仍在使用的先前销售产品的数量，因此需要推断生存曲线。

例如，如果您为特定的智能手机型号制造配件，那么配件销售至少部分取决于仍在使用的智能手机的数量。（这不是作业，顺便说一句。）

细节

我有一些时间序列数据，并希望使用glm或类似的东西来拟合回归模型。因变量和自变量之间的关系是这样的： 回归公式

其中 p 为时间段，y p为因变量，x p为自变量，c 0和 c 1为回归系数，F t为累积分布函数（如pgamma），e p为残差。

通过前三个时间段，该函数将扩展为如下所示：

#y[1] = c0 + c1*(x[1]*(1-integrate(function(q) {pgamma(q, c2, c2/c3)}, 0, 1)$value))
#y[2] = c0 + c1*(x[1]*(1-integrate(function(q) {pgamma(q, c2, c2/c3)}, 1, 2)$value) + x[2]*(1-integrate(function(q) {pgamma(q, c2, c2/c3)}, 0, 1)$value))
#y[3] = c0 + c1*(x[1]*(1-integrate(function(q) {pgamma(q, c2, c2/c3)}, 2, 3)$value) + x[2]*(1-integrate(function(q) {pgamma(q, c2, c2/c3)}, 1, 2)$value) + x[3]*(1-integrate(function(q) {pgamma(q, c2, c2/c3)}, 0, 1)$value))

所以，我有 x p和 y p的历史数据，我想获得最小化残差的系数/参数 c 0、 c 1、 c 2和 c 3的值。

我认为解决方案是使用glm和创建一个自定义系列，但我不知道该怎么做。 我查看了Gamma家庭的代码，但没有走多远。我已经能够使用“手动”进行优化nlminb，但我更喜欢由或类似功能predict提供的简单性和实用性（即和其他） 。glm

以下是一些示例数据：

# Survival function (the integral part):
fsurv<-function(q, par) {
  l<-length(q)
  out<-vapply(1:l, function(i) {1-integrate(function(x) {pgamma(x, par[1], par[1]/par[2])}, q[i]-1, q[i])$value}, FUN.VALUE=0)
  return(out)}

# Sum up the products:
frevsumprod <- function(x,y) {
  l <- length(y)
  out <- vapply(1:l, function(i) sum(x[1:i]*rev(y[1:i])), FUN.VALUE=0)
  return(out)}

# Sample data:
p<-1:24 # Number of periods
x_test<-c(1188, 2742, 4132) # Sample data
y_test<-c(82520, 308910, 749395, 801905, 852310, 713935, 624170, 603960, 640660, 553600, 497775, 444140) # Sample data
c<-c(-50.161147,128.787437,0.817085,13.845487) # Coefficients and parameters, from another method that fit the data

# Pad the data to the correct length:
pad<-function(p,v,padval=0) {
  l<-length(p)
  padv<-l-length(v)
  if(padv>0) (v<-c(v,rep(padval,padv)))
  return(v)
}
x_test<-pad(p,x_test)
y_test<-pad(p,y_test,NA)

y_fitted<-c[0+1]+c[1+1]*frevsumprod(x_test,fsurv(p,c[(2:3)+1])) # Fitted values from regression

library(ggplot2)
ggplot(data.frame(p,y_test,y_fitted))+geom_point(aes(p,y_test))+geom_line(aes(p,y_fitted)) # Plot actual and fit
4

1 回答 1

1

这不能用glm. familyin指定线性预测器如何与 y的glm平均值连接。请参阅?family和维基。特别是，您需要能够family使用（某些）功能编写一个列表，例如：

> fam <- poisson()
> str(fam)
List of 12
 $ family    : chr "poisson"
 $ link      : chr "log"
 $ linkfun   :function (mu)  
 $ linkinv   :function (eta)  
 $ variance  :function (mu)  
 $ dev.resids:function (y, mu, wt)  
 $ aic       :function (y, n, mu, wt, dev)  
 $ mu.eta    :function (eta)  
 $ initialize:  expression({  if (any(y < 0))  stop("negative values not allowed for the 'Poisson' family")  n <- rep.int(1, nobs| __truncated__
 $ validmu   :function (mu)  
 $ valideta  :function (eta)  
 $ simulate  :function (object, nsim)  
 - attr(*, "class")= chr "family"
> 
> fam <- Gamma()
> str(fam)
List of 12
 $ family    : chr "Gamma"
 $ link      : chr "inverse"
 $ linkfun   :function (mu)  
 $ linkinv   :function (eta)  
 $ variance  :function (mu)  
 $ dev.resids:function (y, mu, wt)  
 $ aic       :function (y, n, mu, wt, dev)  
 $ mu.eta    :function (eta)  
 $ initialize:  expression({  if (any(y <= 0))  stop("non-positive values not allowed for the 'gamma' family")  n <- rep.int(1, n| __truncated__
 $ validmu   :function (mu)  
 $ valideta  :function (eta)  
 $ simulate  :function (object, nsim)  
 - attr(*, "class")= chr "family"

其中eta是指线性预测器。即至少您需要指定一个反向链接函数 ，linkinv它仅依赖于通过参数和协变量之间的点积的协变量。你的不是，因为它以非线性方式依赖于 c_2 和 c_3 。

于 2017-10-25T18:16:33.287 回答