r - ユーザー指定のファミリーを使用した glm 回帰のモデル式

Question

バックグラウンド

製品ラインの売上を予測しようとしています (最後のサンプルの y_test)。ある期間の売上は、別の製品 (x_test) の以前のすべての売上と、それらの以前の売上のうちまだ使用されている数に基づいています。ただし、以前に販売された製品がまだ使用されている数を直接測定することはできないため、生存曲線を推測する必要があります。

たとえば、特定のスマートフォン モデル用のアクセサリを作成する場合、アクセサリの売上は少なくとも部分的には、まだ使用されているスマートフォンの数に基づいています。(これは宿題ではありません。)

詳細

glm私はいくつかの時系列データを持っており、または類似のものを使用して回帰モデルを適合させたいと考えています。従属変数と独立変数の関係は次のとおりです。

ここで、p は期間、y _pは従属変数、x _pは独立変数、c ₀と c ₁は回帰係数、F _tは累積分布関数 ( などpgamma)、e _pは残差です。

最初の 3 つの期間を通じて、関数は次のように拡張されます。

#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 ₀、c ₁、c ₂、および c ₃の値を取得したいと考えています。

解決策は、カスタム ファミリを使用して作成することだと思いますが、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