标签: least-squares

我想在不使用"lm"的情况下计算R中的普通最小二乘(OLS)估计,这有几个原因.首先,考虑到数据大小在我的情况下是一个问题,"lm"还计算了许多我不需要的东西(例如拟合值).其次,我希望能够在R中自己实现OLS,然后再使用其他语言(例如,在C中使用GSL).

您可能知道,模型是:Y = Xb + E; 用E~N(0,sigma ^ 2).如下所述,b是具有2个参数的向量,即平均值(b0)和另一个系数(b1).最后,对于我将要进行的每个线性回归,我想要b1(效应大小)的估计,其标准误差,sigma ^ 2的估计(残差方差)和R ^ 2(确定系数).

这是我的数据.我有N个样本(例如个体,N~ = 100).对于每个样本,我有Y个输出(响应变量,Y~ = 10 ^ 3)和X个点(解释变量,X~ = 10 ^ 6).我想分开处理Y输出,即.我想推出Y线性回归:一个用于输出1,一个用于输出2,等等.此外,我想一个一个地使用解释变量:对于输出1,我想在第1点回归,然后在第2点回归,然后......终于在X点了.(我希望它很清楚......!)

这是我的R代码,通过矩阵代数检查"lm"与计算OLS估算的速度.

首先,我模拟虚拟数据:

nb.samples <-  10  # N
nb.points <- 1000  # X
x <- matrix(data=replicate(nb.samples,sample(x=0:2,size=nb.points, replace=T)),
            nrow=nb.points, ncol=nb.samples, byrow=F,
            dimnames=list(points=paste("p",seq(1,nb.points),sep=""),
              samples=paste("s",seq(1,nb.samples),sep="")))
nb.outputs <- 10  # Y
y <- matrix(data=replicate(nb.outputs,rnorm(nb.samples)),
            nrow=nb.samples, ncol=nb.outputs, byrow=T,
            dimnames=list(samples=paste("s",seq(1,nb.samples),sep=""),
              outputs=paste("out",seq(1,nb.outputs),sep="")))

以下是我自己使用的函数:

GetResFromCustomLinReg <- function(Y, xi){ # both Y and xi are N-dim …

我有4个点,非常接近于一个平面 - 它是1,4-二氢吡啶循环.

我需要计算从C3和N1到平面的距离,该平面由C1-C2-C4-C5组成.计算距离是可以的,但拟合平面对我来说很难.

1,4-DHP循环http://i.stack.imgur.com/dhNDo.png

1,4-DHP循环,另一种观点http://i.stack.imgur.com/6Xs0z.png

from array import *
from numpy import *
from scipy import *

# coordinates (XYZ) of C1, C2, C4 and C5
x = [0.274791784, -1.001679346, -1.851320839, 0.365840754]
y = [-1.155674199, -1.215133985, 0.053119249, 1.162878076]
z = [1.216239624, 0.764265677, 0.956099579, 1.198231236]

# plane equation Ax + By + Cz = D
# non-fitted plane
abcd = [0.506645455682, -0.185724560275, -1.43998120646, 1.37626378129]

# creating distance variable
distance =  zeros(4, float)

# calculating distance from point to …

A*(A\D)下面的结果怎么可能不等于D？

它应该产生D- 这是Octave文档的摘录:

线性方程组在数值分析中普遍存在.要求解线性方程组Ax = b,请使用左除法运算符'\':x = A\b

以下是想要尝试的人的代码:

A = [1,1,1;0,0,0;2,1,2;2,1,2;3,5,6]
D = [1;2;3;4;5]
% A is of rank 3:
rank(A)
% therefore the system Ax=D has a unique solution
x = A\D
% but Octave has not given the good solution:
A*x

有人说我Matlab产生完全相同的结果.

编辑10/10/2012:在阅读了答案之后,让我指出我犯了一个严重错误的地方:声明"A是等级3因此系统Ax = D有一个独特的解决方案"是绝对错误的!顺便说一句,上面显示的文档非常令人不安.

我是 python 的新手，并尝试使用 lmfit 包来检查我自己的计算，但是我不确定 (1) 如何包含数据 (sig) 的错误以用于以下测试（和 2）的错误我使用如下所示的 conf_interval2d 获取）：

    import numpy as np
    from lmfit import Parameters, Minimizer, conf_interval, conf_interval2d, minimize, printfuncs


    x=np.array([ 0.18,  0.26,  1.14,  0.63,  0.3 ,  0.22,  1.16,  0.62,  0.84,0.44,  1.24,  0.89,  1.2 ,  0.62,  0.86,  0.45,  1.17,  0.59, 0.85,  0.44])
    data=np.array([ 68.59,  71.83,  22.52,44.587,67.474 ,  55.765,  20.9,41.33783784,45.79 ,  47.88,   6.935,  34.15957447,44.175,  45.89230769,  57.29230769,  60.8,24.24335594,  34.09121287,  42.21504003,  26.61161674])
    sig=np.array([ 11.70309409,  11.70309409,  11.70309409,  11.70309409,11.70309409,  11.70309409,  11.70309409,  11.70309409,11.70309409,  11.70309409,  11.70309409,  11.70309409,11.70309409,  11.70309409,  11.70309409,  11.70309409,11.70309409,  11.70309409,  11.70309409, …

假设我们有两个数据数组：x = [1,2,3] y = [2,4,6]显然，线性拟合将返回2的斜率和0的截距，当然，这两个例程脾气暴躁linalg.lstsq()，polyfit()并且成功。但是他们认为坡度和截距是要搜索的参数。

是否可以保持斜率固定并仅定义截距？

我需要解决一个线性方程组Lx = b，其中x始终是一个向量（3x1数组），L是一个Nx3数组，b是一个Nx1向量。N通常在4到10之类的范围内。

scipy.linalg.lstsq（L，b）

但是，我需要做很多次（大约200x200 = 40000次），因为x实际上是与图像中每个像素相关联的东西。因此，x实际上存储在PxQx3数组中，其中P和Q类似于200-300，最后一个数字“ 3”表示向量x。现在，我只遍历每一行每一行，并逐一求解方程。如何有效地求解那些40000个不同的线性方程组，可能使用某些矢量化技术或其他特殊方法？

谢谢

我想优化 SIR 模型的拟合。如果我只用 60 个数据点拟合 SIR 模型，我会得到一个“好”的结果。“好”意味着，拟合模型曲线接近数据点，直到 t=40。我的问题是，我怎样才能得到更好的拟合，也许基于所有数据点？

ydata = ['1e-06', '1.49920166169172e-06', '2.24595472686361e-06', '3.36377954575331e-06', '5.03793663882291e-06', '7.54533628058909e-06', '1.13006564683911e-05', '1.69249500601052e-05', '2.53483161761933e-05', '3.79636391699325e-05', '5.68567547875179e-05', '8.51509649182741e-05', '0.000127522555808945', '0.000189928392105942', '0.000283447055673738', '0.000423064043409294', '0.000631295993246634', '0.000941024110897193', '0.00140281896645859', '0.00209085569326554', '0.00311449589149717', '0.00463557784224762', '0.00689146863803467', '0.010227347567051', '0.0151380084180746', '0.0223233100045688', '0.0327384810150231', '0.0476330618585758', '0.0685260046667727', '0.0970432959143974', '0.134525888779423', '0.181363340075877', '0.236189247803334', '0.295374180276257', '0.353377036130714', '0.404138746080267', '0.442876028839178', '0.467273954573897', '0.477529937494976', '0.475582401936257', '0.464137179474659', '0.445930281787152', '0.423331710456602', '0.39821360956389', '0.371967226561944', '0.345577884704341', '0.319716449520481', '0.294819942458255', '0.271156813453547', '0.24887641905719', '0.228045466022105', '0.208674420183194', '0.190736203926912', '0.174179448652951', '0.158937806544529', '0.144936441326754', '0.132096533873646', '0.120338367115739', '0.10958340819268', '0.099755679236243', '0.0907826241267504', '0.0825956203546979', '0.0751302384111894', '0.0683263295744258', '0.0621279977639921', '0.0564834809370572', '0.0513449852139111', '0.0466684871328814', '0.042413516167789', '0.0385429293775096', …

我正在尝试使我的数据适合某个功能。我一直将此示例代码用作指南http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html#example-of-solving-a-fitting-problem。我的代码如下：

from scipy.optimize import least_squares
import numpy as np
import matplotlib.pyplot as plt

def model(x, u):
    return -x[0] * np.sqrt((x[1]/u) - 1)

def fun(x, u, y):
    return y - model(x, u)

def jac(x, u, y):
    J = np.empty((u.size, x.size))
    J[:, 0] = np.sqrt((x[1]/u) - 1)
    J[:, 1] = x[0] / (2 * u * np.sqrt((x[1]/u) - 1))
    return J

u = np.array(T_h2)
y = np.array(lnR2)
x0 = np.array([0.1,0.2])
res = least_squares(fun, x0, jac=jac, bounds=(0, 100), args=(u, y), verbose=1) …

我正在 Python 2.7 中运行拟合，并lmfit使用一些测试数据和以下代码。我需要重量为 的加权拟合1/y（使用 Leven-Marq. 例程）。我已经定义了权重并在这里使用它们：

from __future__ import division
from numpy import array, var
from lmfit import Model
from lmfit.models import GaussianModel, LinearModel

import matplotlib.pyplot as plt
import seaborn as sns

xd = array([1267, 1268, 1269, 1270, 1271, 1272, 1273, 1274, 1275, 1276,
    1277, 1278, 1279, 1280, 1281, 1282, 1283, 1284, 1285, 1286, 1287, 1288,
     1289, 1290, 1291, 1292, 1293, 1294, 1295, 1296, 1297, 1298, 1299, 1300,
     1301, 1302, 1303, 1304, 1305, 1306, 1307, …

我有这些价值观：

T_values = (222, 284, 308.5, 333, 358, 411, 477, 518, 880, 1080, 1259) (x values)

C/(3Nk)_values = (0.1282, 0.2308, 0.2650, 0.3120 , 0.3547, 0.4530, 0.5556, 0.6154, 0.8932, 0.9103, 0.9316) (y values)

我知道他们遵循模型：

C/(3Nk)=(h*w/(k*T))**2*(exp(h*w/(k*T)))/(exp(h*w/(k*T)-1))**2

我也知道k=1.38*10**(-23)和h=6.626*10**(-34)。我必须找到最能描述测量数据的w。我想在python中使用最小二乘法来解决这个问题，但是我真的不太了解它是如何工作的。谁能帮我？