- html - 出于某种原因,IE8 对我的 Sass 文件中继承的 html5 CSS 不友好?
- JMeter 在响应断言中使用 span 标签的问题
- html - 在 :hover and :active? 上具有不同效果的 CSS 动画
- html - 相对于居中的 html 内容固定的 CSS 重复背景?
假设以下情况:
我有一个变量 Y 的计数数据,我假设它是泊松分布的。我还有同一时间段内变量 X 的数据,每个观察值代表一个特定事件。我假设 Y 的值来自两个不同的影响,所以我将每个观察 Y_i 分成两个泊松分布的 Y_i1 和 Y_i2,但我仍然只有关于总 Y_i 的观察。我还假设事件(由 X 表示)对 Y_i1 有长期影响,并且我有参数 lambda_i2 的估计量。
所以我的回归公式是fml=Y_i ~ b_1*X_i+....+b_n*X_(i-n+1) + offset(lambda_i2) -1
且 n>=24。
这意味着 X 的最后 24 个(或更多,因为长期效应)值以累加方式影响 Y_i1 的值,我没有截距(b_0=0)。
我制作了一个矩阵 m,其行代表 Y_i,它的所有 24 个(或更多)回归量用于 Y_i 的每个观察值和 lambda_i2 的相应估计量。
现在我用glm(fml, family=poisson(link="identity"), data=m)
并针对不同的 n 值 (=24,48,36,...) 进行了尝试。
总是,一些系数收到负值,这在解释中没有意义。 (X所代表的事件只能对Y的值产生正向或无影响。)
这引出了我的问题:
如何在我的模型中使用约束 b_i >=0
?
在我之前的研究中,我找到了函数 glmc(),但我不确定如何在此处包含我的约束。
作为替代方案,我还考虑过以贝叶斯方式分析该模型,但我还没有找到用于泊松分布的 glm() 的贝叶斯版本,以便我可以自己指定 b_i 的先验。 (然后我可以在之前包括积极性。)
你有什么想法吗?
这是我的数据和代码的摘录:
y=c(279,623, 1025, 1701, 1862, 2544, 2308, 2231, 2234, 2550, 2698, 2805, 3510, 3032, 2746, 2074, 1062, 513, 226, 116, 87, 79, 116, 335, 594, 1081, 1425, 1775, 2056, 2387, 2337, 2354, 2665, 2406, 2433, 2550, 2820, 3655, 4566, 2330, 1267, 531, 280, 148, 92, 89, 141, 458, 852, 1214)
X=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.88349136, 0.54951680, 0.13306912, 0.15321180, 0.00000000, 1.42569128, 0.55808054, 0.65486418, 0.27530564, 0.24813572, 0, 0, 2.09889028, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.18947898, 0.17347032, 0.94538886, 0.03334654, 0.05593732, 0.00000000, 0.99772264, 0.11121918, 0, 1.41673120, 0.27375384, 0, 0, 0, 0, 5.67487576, 0, 0, 0, 0, 0, 0, 0, 0, 1.55642510, 0.98419866, 0.50992652)
lambda=c(253.5, 562.5, 1053.0, 1645.0, 2064.5, 2215.0, 2503.0, 2443.0, 2514.5, 2701.0, 2972.5, 3035.5, 3422.5, 3295.0, 2882.5, 2094.0, 1211.0, 579.5, 265.5, 155.0, 112.5, 82.5, 117.5, 306.0, 627.0, 1021.0, 1463.5, 1722.5, 2017.5, 2146.5, 2209.0, 2231.5, 2265.0, 2320.0, 2442.0, 2507.0, 2957.0, 3674.0, 3345.5, 2285.0, 1265.0, 555.5, 252.0, 145.5, 86.5, 90.5, 148.0, 362.0, 738.0, 1137.5)
regressors=function(n,x){
m=length(x)-n+1;
r=matrix(0,m,n);
for (i in 0:(n-1)){ r[,(i+1)]=x[(n-i):(length(x)-i)]}
return(r);
}
r=regressors(24,X);
reg=cbind(y,data.frame(r),lambda);
fml=as.formula(paste("y~", paste(colnames(reg)[2:25], collapse = "+"), "+offset(lambda)-1"));
g=glm(fml, poisson(link="identity"), data=reg); %this leads to negative coefficients
obj=function(b){-sum(dpois(y, r%*%b, log=TRUE))}
st=coef(lm(fml, data=reg));
opt=optim(st, obj); % this is where the error occurs
regressors() 是我编写的用于计算回归量的函数(它产生一个有 n 列和 50 行的矩阵 r,每行 i 代表 y_i 的回归量)。
最佳答案
尝试首要原则。
# generate random input data
set.seed(123)
n <- 100
x <- 1:n
X <- cbind(1, x)
b <- c(0.1, 3)
y <- rpois(n, X %*% b)
# log likelihood objective function
obj <- function(b) -sum(dpois(y, X %*% b, log = TRUE))
# as a check try with no constraints - these two should give the same coefs
glm(y ~ X + 0, family = poisson("identity"))
st <- coef(lm(y ~ X + 0)); optim(st, obj)
# now add lower bounds ensuring starting value is in feasible region
optim(pmax(st, 1), obj, lower = c(0, 0), method = "L-BFGS-B")
注意 1:如果您的参数估计值位于上述约束示例中的可行区域边界上,请小心。
注意 2: 这里是对后来添加到问题中的代码的修改。 y
、X
、lambda
和 regressors
与问题相同。请注意,添加约束会强制许多系数为零。
r <- regressors(24,X)
reg <- cbind(y,data.frame(r),lambda)
fml <- y ~ . - lambda + offset(lambda)
# check that g and opt give the same coefs
g <- glm(fml, poisson(link = "identity"), data = reg)
obj <- function(b)-sum(dpois(y, r%*%b + lambda, log = TRUE))
st <- coef(lm(fml, data=reg))
opt <- optim(st, obj, method = "BFGS", control = list(maxit = 500))
optim(pmax(coef(g), 1), obj, method = "L-BFGS-B", lower = 0 * st)
关于r - 如何将泊松回归的系数限制为 R 中的正数?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/34976755/
所以我目前正在代码 war 中解决以下问题: 返回一个数组,其中第一个元素是正数的计数,第二个元素是负数的总和。如果输入数组为空或 null,则返回一个空数组。 我想出了下面的代码,它并不漂亮,但我知
我想禁用文本框中的字母,这就是我使用此表达式的原因: /^[0-9]*$/ 问题是我也想允许负数, 这个表达式不允许我使用减号。(-)..我该怎么办? 最佳答案 尝试使用以下正则表达式: /^-?[0
这是一道面试题。我们有一个正整数数组,我们必须重新排列和连接数组元素,以便生成的数字是可以使用此数组形成的最大数字。 例如: [884 88] -> 88884 [20 19 90] -> 90201
所以我有点无奈。我想我了解如何对 IEEE-754 数字进行加法,但我在减法方面遇到了很大的问题。我感觉我的想法是对的,但我想和大家确认一下。 因此我们有以下两个 IEEE-754 编号: x: 0
我想将一个正数转换成相应的数字列表——数字也应该是整数。 转换时,比如1024,应该返回[1,0,2,4] 最佳答案 在 Swift 4.1 或更高版本中 let number = 1024 let
如果 float 可以用来正确表示我将使用的参数范围,我正在尝试为我的音频应用程序弄清楚。 它需要的“最大”掩码用于频率参数,它是正的,并且允许最多两位数字作为尾数(即从 20.00 hz 到 220
我是一名优秀的程序员,十分优秀!