1

我正在使用RF# 中的类型提供程序来访问一些与回归相关的 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适用，但是我无法应用这些约束。

4

1 回答 1

0

也许不完全是您问题的答案，因为它要求在 R 中进行优化。但也许以下内容会有所帮助。它无论如何都使用 NLopt 库，我认为这是 R 使用的？如果您在制定 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
于 2016-07-23T05:41:57.793 回答