最小二乘法原理及其代码实现

news2024/11/24 10:02:43

一、最小二乘法原理

假设目前我们有一些数据,x是输入,y是与之对应的输出。现在想利用这些已有的数据,从中发现出规律,来预测没有出现过的输入会产生什么样的输出。

假设系统为单输入单输出系统,我们想在这个系统里找到数据背后的规律。规律需要通过模型来进行表征。为了表征规律可以使用不同的手段,不同的手段所建立的模型各有差异,有的模型精确度高但是使用麻烦,有的模型精确度欠缺但是使用简便。下面要介绍典型的建模方法——最小二乘法。

这里说的建模是指建立数学模型,即通过数学表达式来表征规律。通常最容易想到的表达式就是一次函数,二次函数了。

设模型:y=kx+b ,通过数据集(x_i,y_i)(i可以是很多数)来确定k,b

确定k,b只需要两对数据就够了,两点确定一条直线,数据多了会过拟合的。我们采集了很多数据,下面考虑如何将多采集的数据利用起来。

下面有两种思路:

方法一:更换模型,将一次换成高次,如:y=ax^2+bx+c

方法二:仍用原模型,尽可能让各点靠近直线,或者说找到一组k和b使直线周围分布尽可能多的点。

对于方法一,将一次换成二次,三次还可以接受,但是当数据很多时用该方法计算量太大了,不可取。

最小二乘法采用的就是方法二的思想,目标是让各点靠近直线。下面接着思考用什么评价指标来进行衡量。

下面引出了“误差平方和”的概念。为了完成上面的目标,先想到所有点到直线的距离求和,但是若使用点到直线距离公式会出现根号,计算会比较麻烦,所以对各点距离取平方再求和,这样计算量就小很多了。

L=\sum_{i=1}^{n} (y_i-f(x_i))^2

L称为目标函数,即上面所提到的最小平方和。通过求导就能求出它的极小值。

推导过程就省略了,下面直接给出结论:

\theta=(X^TX)^{-1}X^TY

假设有五个样本,要拟合的函数为:

 y=ax^3+bx^2+cx+d

则各参数含义为:

\theta=[d,c,b,a]^T

X = \begin{bmatrix} 1 & x_1 & x_1^2 & x_1^3 \\ 1 & x_2 & x_2^2 & x_2^3 \\ 1 & x_3 & x_3^2 & x_3^3 \\ 1 & x_4 & x_4^2 & x_4^3 \\ 1 & x_5 & x_5^2 & x_5^3 \\ \end{bmatrix}

Y=[y_1,y_2,y_3,y_4]^T

二、MATLAB代码实现

function main()
    % 主程序
    [X, y] = generate_data(100);
    degree = 3;
    X_poly = polynomial_features(X, degree);
    X_new = linspace(0, 2, 100)';
    X_new_poly = polynomial_features(X_new, degree);
    theta_best = least_squares_fit(X_poly, y);

    disp('Theta values:');
    disp(theta_best');

    y_predict = predict(X_new_poly, theta_best);

    plot_data_and_predictions(X, y, X_new, y_predict);
end

function [X, y] = generate_data(num_points)
    % 生成示例数据
    rng(0); % 设置随机数生成器的种子
    X = 2 * rand(num_points, 1);
    y = 4 + 3 * X + randn(num_points, 1);
end

function X_poly = polynomial_features(X, degree)
    % 为给定的X生成多项式特征
    X_poly = [ones(size(X, 1), 1)]; % 添加偏置项(常数项)
    for d = 2:degree+1
        X_poly = [X_poly, X.^d]; % 添加多项式项
    end
end

function theta_best = least_squares_fit(X_poly, y)
    % 使用最小二乘法拟合模型
    theta_best = X_poly' * X_poly \ X_poly' * y;
end

function y_predict = predict(X_poly, theta)
    % 使用拟合的模型进行预测
    y_predict = X_poly * theta;
end

function plot_data_and_predictions(X, y, X_new, y_predict)
    % 可视化数据和预测结果
    figure;
    plot(X_new, y_predict, 'r-', 'DisplayName', 'Prediction');
    hold on;
    plot(X, y, 'b.', 'DisplayName', 'Data');
    hold off;
    xlabel('X');
    ylabel('y');
    title('Polynomial Regression Fit');
    legend('show');
end

三、python代码实现 

import numpy as np
import matplotlib.pyplot as plt


def generate_data(num_points=100):
    """
    生成示例数据
    """
    np.random.seed(0)
    X = 2 * np.random.rand(num_points, 1)
    y = 4 + 3 * X + np.random.randn(num_points, 1)
    return X, y


def polynomial_features(X, degree):
    """
    为给定的X生成多项式特征
    """
    # 在多项式特征中添加偏置项(常数项)
    X_poly = np.c_[np.ones((X.shape[0], 1))]
    # range的范围是1-degree,不包括degree+1
    for d in range(1, degree + 1):
        # np.c_ 是 NumPy 中用于按列合并数组的对象
        X_poly = np.c_[X_poly, X ** d]
    return X_poly


def least_squares_fit(X_poly, y):
    """
    使用最小二乘法拟合模型
    """
    # 使用公式求theta
    theta_best = np.linalg.inv(X_poly.T.dot(X_poly)).dot(X_poly.T).dot(y)
    return theta_best


def predict(X_poly, theta):
    """
    使用拟合的模型进行预测
    """
    return X_poly.dot(theta)


def plot_data_and_predictions(X, y, X_new, y_predict):
    """
    可视化数据和预测结果
    """
    plt.plot(X_new, y_predict, "r-", label="Prediction")
    plt.plot(X, y, "b.", label="Data")
    plt.xlabel("X")
    plt.ylabel("y")
    plt.title("Polynomial Regression Fit")
    plt.legend()
    plt.show()


# 主程序
if __name__ == "__main__":
    # 生成数据
    X, y = generate_data()

    # 设置多项式的次数
    degree = 3

    # 为训练数据和新数据生成多项式特征
    X_poly = polynomial_features(X, degree)
    X_new = np.linspace(0, 2, 100).reshape(100, 1)
    X_new_poly = polynomial_features(X_new, degree)

    # 使用最小二乘法拟合模型
    theta_best = least_squares_fit(X_poly, y)

    # 打印计算出的theta值,按顺序显示常数项、一次项、二次项...系数
    print("Theta values:", theta_best.flatten())

    # 进行预测
    y_predict = predict(X_new_poly, theta_best)

    # 可视化结果
    plot_data_and_predictions(X, y, X_new, y_predict)

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/1813178.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

【TB作品】MSP430F5529 单片机,数字时钟设计与实现,整点时通过蜂鸣器播放音乐进行报时

基于单片机的数字时钟设计与实现 作品名称 基于MSP430单片机的OLED显示数字时钟 作品功能 本作品实现了一个具有时间显示和整点报时功能的数字时钟。通过OLED屏幕显示当前时间,用户可以通过按键设置时间,并在整点时通过蜂鸣器播放音乐进行报时。 作…

Serif Affinity 2.5 (macOS, Windows) - 专业创意软件

Serif Affinity 2.5 (macOS, Windows) - 专业创意软件 Affinity Designer 2, Affinity Photo 2, Affinity Publisher 2 请访问原文链接:Serif Affinity 2.5 (macOS, Windows) - 专业创意软件,查看最新版。原创作品,转载请保留出处。 作者主…

如何优化仓库布局与ERP库存管理

一、引言 随着企业规模的不断扩大,仓库管理和库存控制成为企业运营中不可或缺的一环。优化仓库布局和提高ERP库存管理效率,对于降低企业成本、提高物流效率、增强企业竞争力具有重要意义。 二、优化仓库布局 1. 分析仓库需求 在优化仓库布局之前&…

【学习笔记】Linux

Linux 1、 介绍 1.1、概述 1.2、特点 1.3、Linux的发行版2、 基础篇 —— 文件系统 2.1、文件系统 2.2、目录结构3、 基础篇 —— VI/VIM 编辑器 3.1、概述 3.2、编辑器模式及常用命令4、 基础篇 —— 网络配置 4.1、VMware NetWork …

【Linux】shell——条件判断test,各种运算符,expr

条件判断——test 真——0 假——1 test expression or [ expression ] 整数运算符 字符串运算符 -z 长度是否为0 -n 长度是否不为0 str1 str2 str1 ! str2 补 &&-->逻辑与,前面为真后面才会执行 || -->逻辑或,前面为假后面才…

【算法实战】每日一题:18.1并查集知识点讲解以及算法实战

1.题目 给定一个序列,通过n-1次相邻元素的合并操作,恢复原始序列。 2.涉及知识点 - 并查集 (Union-Find) 并查集 (Union-Find) 详解 概述 并查集(Union-Find),也称为不相交集数据结构,用于处理一些不相…

MFC案例:利用SetTimer函数编写一个“计时器”程序

一、希望达成效果 利用基于对话框的MFC项目,做一个一方面能够显示当前时间;另一方面在点击开始按钮时进行读秒计时,计时结果动态显示,当点击结束时读秒结束并保持最后结果。 二、编程步骤及相关代码、注释 1、启动VS…

OA协同办公系统 iWebPDF插件安装

1、下载压缩文件 iweboffice,并进行解压 链接:https://pan.baidu.com/s/1GQd7000PTZ771ifL5KEflg 提取码:hb56 2、安装iWenpdf2018.exe 3、安装金格中间件外部应用 4、测试了谷歌、360安全,发现安装插件后,只有360极…

BP8519C非隔离降压型恒压芯片

BP8519封装和丝印 BP8519封装和丝印 注意: 该芯片为非隔离ACDC电源芯片,非专业人员请勿使用。专业人员在使用时必须注意防护,避免触电。 非隔离ACDC电源芯片,国内有多家半导体厂商生产,在部分追求低价格的低端仪表、灯…

vivado HW_SIO_GTGROUP、HW_SIO_IBERT

HW_SIO.GTGROUP 描述 GT组与硬件设备上的GT IO组相关,具有可用的数量 GT引脚和组由目标Xilinx FPGA确定。在Kintex-7 xc7k325部件上,用于 例如,有四个GT组,每个组包含四个差分GT引脚对。每个GT pin有自己的接收器hw_sio_rx和发射器…

人工智能GPT-4o?

对比分析 在讨论GPT-4o时,我们首先需要了解其前身,即GPT-4,以及其之前的版本。GPT系列从GPT-1到GPT-4经历了多次迭代,每一次都带来了显著的进步。 GPT-4 vs GPT-4o: 1. **参数规模:** GPT-4o在参数规模上…

PyTorch 张量数据类型

【数据类型】Python 与 PyTorch 常见数据类型对应: 用 a.type() 获取数据类型,用 isinstance(a, 目标类型) 进行类型合法化检测 >>> import torch >>> a torch.randn(2,3) >>> a tensor([[-1.7818, -0.2472, -2.0684],[ 0.…

iOS ------ 对象的本质

一,OC对象本质,用clang编译main.m OC对象结构都是通过基础的C/C结构体实现的,我们通过创建OC文件及对象,将OC对象转化为C文件来探寻OC对象的本质。 代码: interface HTPerson : NSObject property(nonatomic,strong)…

什么是SOLIDWORKS科研版

随着科技的不断进步,工程设计和科学研究变得越来越复杂,需要更强大的工具来满足需求。SOLIDWORKS科研版就是在这样的背景下诞生的,它为科研人员和工程师提供了一套全方面、快捷的解决方案,以应对各种科研和工程挑战。 SOLIDWORKS科…

Surface安装Windows和Ubuntu双系统方法(包括Ubuntu适配触控屏的方法)

这是一个目录0.0 前言让我们从一块砖头开始现在你有了能进入windows系统的surface并且想安装Ubuntu现在Ubuntu也有了再见 前言 之前我的Surface装上Ubuntu了好好的,能用,但是Ubuntu原本的内核是不支持很多Surface的功能的,比如触控屏&#xf…

串口调试助手软件(ATK-XCOM) 版本:v2.0

串口设置 软件启动后,会自动搜索可用的串口,可以显示详细的串口信息,由于兼容性原因某些电脑可能不会显示。 超高波特率接收,在硬件设别支持的情况下,可自定义波特率,点“自定义”即可输入您想要的波特率&…

macOS 15 beta (24A5264n) Boot ISO 原版可引导镜像下载

macOS 15 beta (24A5264n) Boot ISO 原版可引导镜像下载 iPhone 镜像、Safari 浏览器重大更新、备受瞩目的游戏和 Apple Intelligence 等众多全新功能令 Mac 使用体验再升级 请访问原文链接:https://sysin.org/blog/macOS-Sequoia-boot-iso/,查看最新版…

PR基础常识

Pr主要就是用来做视频后期剪辑的。是一款非线性视频剪辑软件。它可以将原有视频作为素材,导入到软件的时间线轨道面板中,对视频进行重新剪辑编排,并可以添加文字、图片、音频等素材文件,也能预设各种效果,让剪辑的视频…

外卖抢单神器

在现代快节奏的生活中,外卖服务已成为许多人日常生活的一部分,给外卖行业带来前所未有的机遇和挑战。随着市场竞争的加剧,许多外卖员开始寻求方法以提升接单效率。但在此过程中,道德和合规性是业务持续性的关键。 正直的经营不仅…

Apple Intelligence 带来的十大影响:人工智能的iPhone时刻到来

引言 在最近的WWDC大会上,Apple发布了全新的Apple Intelligence,引起了全球的广泛关注。这次发布被誉为“人工智能的iPhone时刻”,标志着我们每个人都将拥有第一个AI助理,并将引领AI Agent进入红海时代。本文将详细分析Apple Int…