具有實際概率的邏輯回歸∈(𝑎,𝑏)∈(a,b)in(a,b)在哪裡0<𝑎<𝑏<10<a<b<10<a<b`<1
使用邏輯回歸對概率建模時 1 ，擬合概率的範圍是 (0,1) . 邏輯函數 2 漸近線 0 和 1 ，所以這是一個很好的匹配。然而，在某些應用中，實際範圍 3 的概率可以是 (a,b) 和 0<a<b<1 ，導致尾部出現（潛在的大量）不匹配。
例如，考慮一個所有人都工作並獲得報酬的人群。付款取決於技能、努力和其他我們衡量的東西。每個人還參與抽獎，其貨幣結果正態分佈，期望值為零。個人的收入是工作收入和抽獎結果的總和。對於任何給定的常數 c , 條件概率 p 個人收入大於 c , 以技巧、努力和其他東西為條件，而不是彩票的結果，將滿足 0<a<p<b<1 . （我相信那裡有更好的例子，但這是首先想到的。）
問題：
- 如何調整邏輯回歸模型 a 和 b 給出？（或者不需要調整？）
- 如何調整邏輯回歸模型 a 和 b 沒有給出，但我們知道 0<a<b<1 ?
1 或者就此而言，一個概率模型
2 或者就此而言，標準的普通 CDF
2 感興趣人群中個體的範圍
什麼時候 a 和 b 沒有給出，只需使用通常的邏輯模型（或任何合適的），因為（如果它使用合適的鏈接函數）保證擬合概率不小於下限 0 並且上限不大於 1. 這些界限給出了區間估計 a 和 b.
有趣的問題是什麼時候 a 和 b 是已知的。 您正在娛樂的那種模型似乎如下。您想到了一個單參數分佈系列 F=Fθ 在哪裡 θ 對應於一些“概率”參數。例如， Fθ 可能是伯努利 (θ) 響應時的分佈 Y 是二進制的。
對於與解釋變量向量相關的觀察 x, 響應模型 Yx 然後採取形式
Yx∼Fθ(x);θ(x)=g(xβ)
對於一些“反向鏈接函數” g 我們必須指定：它是模型的一部分。例如，在邏輯回歸中， g 經常被認為是由定義的邏輯函數
g(x)=11+exp(−x).
無論細節如何，製作時 n 獨立觀察 yi （每個都與一個向量相關聯 xi ) 假設符合這個模型，他們的可能性是
L(β)=n∏i=1Pr(Yxi=yi∣θ(xi)=g(xiβ))
你可以像往常一樣繼續最大化它。（豎筆僅表示它後面的參數值決定使用哪個概率函數：它不是條件概率。）
讓 ˆβ 是相關的參數估計。的預測條件分佈 Yi 因此是
Yi∼Fˆθ(xi);ˆθ(xi)=g(xiˆβ).
當圖像 g 包含在區間內 [a,b], 那麼顯然每個 ˆθ(x) 無論如何，也在那個區間 x 也許。（也就是說，這個結論既適用於 x 在數據集中並用於外推到其他 x. )
一個有吸引力的選擇 g 簡單地重新調整通常的邏輯函數，
g(x;a,b)=g(x)−ab−a.
**將此視為出發點：**像往常一樣，探索性分析和擬合優度測試將幫助您確定這是否適合 g.
供以後使用，**請注意 g 和 g(;a,b) 有一個比看起來更複雜的關係，**因為最終它們被用來確定 ˆβ 通過他們的論點 xβ. 因此，這種關係的特徵在於函數 x→y 取決於
g(x)=g(y;a,b)=g(y)−ab−a,
有溶液（如果 g 是可逆的，通常是這樣）
y=g−1((b−a)g(x)+a).
除非 g 最初是線性的，這通常是非線性的。
**為了解決此線程中其他地方表達的問題，**讓我們比較使用獲得的解決方案 g 和 g(;a,b). 考慮最簡單的情況 n=1 觀察和需要估計參數向量的標量解釋變量 β=(β1). 認為 F 是二項式的家庭 (10,θ) 分佈，讓 x1=(1), 並想像 Yi=9 被觀察到。寫作 θ 為了 θ(x1), 可能性是
L(β)=(109)θ9(1−θ)1;θ=g((1)(β1))=g(β1).
L 最大化時 g(β1)=θ=9/10, 用獨特的解決方案ˆβ=g−1(9/10)=log(9/10/(1/10))=log(9)≈2.20.
現在讓我們假設 a=0 和 b=1/2: 也就是說，我們假設 θ≤1/2 不管什麼價值 x 可能有。隨著縮放版本 g 我們和以前一樣計算，只是代入 g(;a,b) 為了 g:
L(β;0,1/2)=(109)θ9(1−θ)1;θ=g((1)(β1);0,1/2)=g(β1;0,1/2).
這不再最大化 θ=9/10, 因為這是不可能的 g(θ;0,1/2) 超過 1/2, 按設計。 L(β;0,1/2) 被最大化的任何 β 那會讓 θ 盡可能接近 9/10; 這發生在 β 任意變大。那麼，使用受限反向鏈接函數的估計是
ˆβ=∞.
顯然兩者都沒有 ˆθ 或者 ˆβ 是原始（無限制）估計的任何簡單函數；特別是，它們與任何重新縮放無關。
**這個簡單的例子暴露了整個程序的一個危險：**當我們假設什麼 a 和 b （以及關於模型的所有其他內容）與數據不一致，我們可能會得到模型參數的奇怪估計 β. 這就是我們付出的代價。
但是，如果我們的假設是正確的，或者至少是合理的呢？ 讓我們用 b=0.95 代替 b=1/2. 這次， ˆθ=9/10 確實最大化似然性，由此估計 β 滿足
910=g(ˆβ;0,0.95)=g(ˆβ)−00.95−0,
所以
g(ˆβ)=0.95×910=0.855,
蘊涵
ˆβ=log(0.855/(1−0.855))≈1.77.
在這種情況下， ˆθ 不變但_ ˆβ 發生了複雜的變化（ 1.77 不是重新縮放的版本 2.20 ).
在這些例子中， ˆθ 當原始估計不在區間內時必須更改 [a,b]. 在更複雜的例子中，它可能必須改變才能改變其他觀察值的估計值 x. 這是效果之一 [a,b] 限制。另一個影響是，即使限制改變了估計的概率， ˆθ, 原始反向鏈接之間的非線性關係 g 和受限鏈接 g(;a,b) 在參數估計中引起非線性（並且可能是複雜的）變化 ˆβ.
**為了說明，**我根據這個模型創建了數據 β=(4,−7) 和限制 a=1/10 和 b=1/2 為了 n 解釋值的等距值 x 之間 0 和 1 包含，然後使用普通邏輯回歸（無約束）擬合它們一次，並使用縮放反向鏈接方法再次使用已知約束。
以下是結果 n=12 二項式 (8,θ(x)) 觀察結果（實際上反映了 12×8=96 獨立的二進制結果）：
這已經提供了洞察力：模型（左圖）預測接近上限的概率 b=1/2 對於小 x. 隨機變化導致某些觀測值的頻率大於 1/2. 在沒有任何限制的情況下，邏輯回歸（中圖）傾向於預測更高的概率。類似的現象發生在大 x.
受限模型極大地改變了估計的斜率 −3.45 到 −21.7 為了使預測保持在 [a,b]. 發生這種情況的部分原因是它是一個小數據集。
*直觀地說，*更大的數據集應該產生更接近底層（真實）數據生成過程的結果。人們也可能期望不受限制的模型運行良好。可以？為了檢查，我創建了一個大一千倍的數據集： n=1200 二項式的觀察 (80,θ(x)) 回复。
當然，正確的模型（右面板）現在非常適合。然而，觀察到的頻率的隨機變化仍然導致普通邏輯模型超出限制。
顯然，當假定值 a 和 b 是（接近）正確的並且鏈接函數的形狀大致正確，最大似然效果很好 - 但它絕對不會產生與邏輯回歸相同的結果。
為了提供完整的文檔，這裡是
R生成第一個圖的代碼。更改12為1200和8生成80第二個圖形。# # Binomial negative log likelihood. # logistic.ab <- function(x, a=0, b=1) { a + (b - a) / (1 + exp(-x)) } predict.ab <- function(beta, x, invlink=logistic.ab) { invlink(cbind(1, x) %*% beta) } Lambda <- function(beta, n, k, x, invlink=logistic.ab, tol=1e-9) { p <- predict.ab(beta, x, invlink) p <- (1-2*tol) * p + tol # Prevents numerical problems - sum((k * log(p) + (n-k) * log(1-p))) } # # Simulate data. # N <- 12 # Number of binomial observations x <- seq(0, 1, length.out=N) # Explanatory values n <- rep(8, length(x)) # Binomial counts per observation beta <- c(4, -7) # True parameter a <- 1/10 # Lower limit b <- 1/2 # Upper limit set.seed(17) p <- predict.ab(beta, x, function(x) logistic.ab(x, a, b)) X <- data.frame(x = x, p = p, n = n, k = rbinom(length(x), n, p)) # # Create a data frame for plotting predicted and true values. # Y <- with(X, data.frame(x = seq(min(x), max(x), length.out=101))) Y$p <-with(Y, predict.ab(beta, x, function(x) logistic.ab(x, a, b))) # # Plot the data. # par(mfrow=c(1,3)) col <- hsv(0,0,max(0, min(1, 1 - 200/N))) with(X, plot(x, k / n, ylim=0:1, col=col, main="Data with True Curve")) with(Y, lines(x, p)) abline(h = c(a,b), lty=3) # # Reference fit: ordinary logistic regression. # fit <- glm(cbind(k, n-k) ~ x, data=X, family=binomial(link = "logit"), control=list(epsilon=1e-12)) # # Fit two models: ordinary logistic and constrained. # for (ab in list(c(a=0, b=1), c(a=a, b=b))) { # # MLE. # g <- function(x) logistic.ab(x, ab[1], ab[2]) beta.hat <- c(0, 1) fit.logistic <- with(X, nlm(Lambda, beta.hat, n=n, k=k, x=x, invlink=g, iterlim=1e3, steptol=1e-9, gradtol=1e-12)) if (fit.logistic$code > 3) stop("Check the fit.") beta.hat <- fit.logistic$estimate # Check: print(rbind(Reference=coefficients(fit), Constrained=beta.hat)) # Plot: Y$p.hat <- with(Y, predict.ab(beta.hat, x, invlink=g)) with(X, plot(x, k / n, ylim=0:1,, col=col, main=paste0("Fit with a=", signif(ab[1], 2), " and b=", signif(ab[2], 2)))) with(Y, lines(x, p.hat, col = "Red", lwd=2)) with(Y, lines(x, p)) abline(h = c(a, b), lty=3) } par(mfrow=c(1,1))
