制限付き回帰式R

Question

F＃のR型プロバイダを使用して、回帰関係のあるR関数にアクセスしています。 回帰係数に制約があると、その加重平均が0になるように回帰を見積もりたいと思います。加重の合計は1です。以下の例は、さまざまな重みを持つ数十の係数を持つため、のみ以下のRコードを示しています。

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) 

予想通りこれは

> lm(y~x1+x2) 

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

Coefficients: 
(Intercept)   x1   x2 
    0.03468  -0.01460   NA 

の出力を提供します。しかし、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