제한된 회귀 식 R

Question

F #의 R 유형 공급자를 사용하여 회귀 관련 R 기능에 액세스하고 있습니다. 회귀 계수에 제약 조건이있을 때 회귀 계수를 제한하여 가중 평균이 0이되도록 계산하고 싶습니다. 가중치의 합이 1입니다. 아래 예는 다양한 가중치가있는 수십 개의 계수가 있으므로 단순화되었습니다. 예상대로이

> lm(y~x1+x2) 

Call: 
lm(formula = y ~ x1 + x2) 

Coefficients: 
(Intercept)   x1   x2 
    0.03468  -0.01460   NA 

의 출력을 제공

y1 <- runif(n = 50,min = 0.02,max=0.05) 
y2 <- runif(n=50,min=0.01,max=0.03) 
y <- c(x1,x2) 
x1 <- c(rep(0,50),rep(1,50)) 
x2 <- c(rep(1,50),rep(0,50)) 
lm(y~x1+x2) 

: 만 아래 R 코드를 보여줍니다. 그러나 나는 x1과 x2에 제약 조건을 두어 가중 평균이 (0.5 * x1 + 0.5 * x2) = 0이되도록하고 싶습니다. 이 경우 절편은 mean(y) = 0.02737966이되고 x1 및 x2 계수는이 값 (-0.006 및 +0.007 각각)의 오프셋을 나타냅니다. quadprog 및 mgcv 패키지는 적용 가능하지만 제약 조건을 적용 할 수 없었습니다.

Answer 1

R에서 최적화 작업을 수행해야하므로 질문에 대한 답변이 정확하지 않을 수도 있습니다.하지만 다음 내용이 도움이 될 수 있습니다. 어쨌든 R이 사용하는 NLopt 라이브러리를 사용합니다. MLE를 공식화하는 데 도움이 필요하면 알려주지 만, 가우스 가정과 선형성이없는 선형 모델에 대해서는 내적 복잡성이 없어야합니다.

LN_COBYLA가 사용자 제공 그라디언트를 사용하지 않더라도 cFunc 및 oFunc의 패턴 일치는 무시합니다. 나는 LD_LBFGS를 시도했지만 AddEqualZeroConstraint()를 지원하지 않는다.

[편집]

완전한 예제를 추가하면 템플릿으로 사용할 수 있습니다. 그것의 관용적이고, 아주 못생긴하지만 요점을 설명합니다. 그러나이 예에서 제약 조건은이를 퇴화시킵니다. NLOptNet, MathNet.Numerics, Fsharp Charting이 필요합니다. 어쩌면 F #에서 제한된 최적화를하려는 다른 사람들을 도울 수 있습니다.

open System 
open System.IO 
open FSharp.Core.Operators 
open MathNet.Numerics 
open MathNet.Numerics.LinearAlgebra 
open MathNet.Numerics.LinearAlgebra.Double 
open MathNet.Numerics.Distributions 
open DiffSharp.Numerical.Float64 
open NLoptNet 

let (.*) (m1 : Matrix<float>) (m2 : Matrix<float>) = 
    m1.Multiply(m2) 

let (.+) (m1 : Matrix<float>) (m2 : Matrix<float>) = 
    m1.Add(m2) 

let (.-) (m1 : Matrix<float>) (m2 : Matrix<float>) = 
    m1.Subtract(m2) 


let V = matrix [[1.;  0.5; 0.2] 
       [0.5;  1.; 0.] 
       [0.2;  0.; 1.]] 
let dat = (DenseMatrix.init 200 3 (fun i j -> Normal.Sample(0., 1.))) .* V.Cholesky().Factor 
let y = DenseMatrix.init 200 1 (fun i j -> 0.) 
let x0 = DenseMatrix.init 200 1 (fun i j -> 0.) 
let x1 = DenseMatrix.init 200 1 (fun i j -> 0.) 
for i in 0 .. 199 do 
    y.[i, 0] <- dat.[i, 0] 
    x0.[i, 0] <- dat.[i, 1] 
    x1.[i, 0] <- dat.[i, 2] 

let ll (th : float array) = 
    let t1 = x0.Multiply(th.[0]) .+ x1.Multiply(th.[1]) 
    let res = (y .- t1).PointwisePower(2.) 
    res.ColumnAbsoluteSums().Sum()/200. 

let oFunc (th : float array) (gradvec : float array) = 
    match gradvec with 
    | null ->() 
    | _  -> (grad ll th).CopyTo(gradvec, 0) 

    ll th 

let cFunc (th : float array) (gradvec : float array) = 
    match gradvec with 
    | null ->() 
    | _  -> (grad ll th).CopyTo(gradvec, 0) 

    th.[0] + th.[1] 

let fitFunc() = 
    let solver = new NLoptSolver(NLoptAlgorithm.LN_COBYLA, uint32(2), 1e-7, 100000) 
    solver.SetLowerBounds([|-10.; -10.;|]) 
    solver.SetUpperBounds([|10.; 10.;|]) 
    //solver.AddEqualZeroConstraint(cFunc) 
    solver.SetMinObjective(oFunc) 
    let initialValues = [|1.; 2.;|] 
    let objReached, finalScore = solver.Optimize(initialValues) 
    objReached |> printfn "%A" 
    let fittedParams = initialValues 
    fittedParams |> printfn "%A" 
    fittedParams 

let fittedParams = fitFunc() |> DenseVector 
let yh = DenseMatrix.init 200 1 (fun i j -> 0.) 
for i in 0 .. 199 do 
    yh.[i, 0] <- dat.[i, 1] * fittedParams.[0] + dat.[i, 2] * fittedParams.[1] 


Chart.Combine([Chart.Line(y.Column(0), Name="y") 
       Chart.Line(yh.Column(0), Name="yh") 
       |> Chart.WithLegend(Title="Model", Enabled=true)]) 
       |> Chart.Show