首先申明本人是土木专业的,因为有需要要用到matlab中的拟合用途,今天好好学习了一些关于matlab多变量拟合的东西,从网上下载了一些程序,也运行了一下,就举一些实例,附上源程序吧,主要是两个自变量和三个自变量,一个因变量的拟合。让自己也更清楚,以后用起来也方便。
    原理就是给出一个自变量和因变量的矩阵,然后给出一个自己认为的带有未知数的拟合方程,然后付一组初始值,根据matlab返回的初始值和误差在附一组初始值,知道最后的相关系数较大,也就是误差较小时,就能拟合的比较好,写出拟合后的方程了。
1.广义线性回归拟合和源码(两个自变量,一个因变量,非线性拟合)
【例】这里有这样一组数据,涉及三个变量:p,t 和z,要拟合出 z = f(p,t) 的关系式(非线性的)。
 z    p                 0.8                    1                        1.2
t               
 60              9.73875                  20.75                  36.5987         
120            13.5725                  29.6325                50.93875
180            18.97875                36.59875              80.13875
240             2075125                38.22125                90.925
300            22.055                     44.58                 104.7725
为了使得回归分析的结果更加直观,我调用regstats函数,编写了一个更为实用的函数:reglm,代码如下(代码中有调用方法和例子)。
首先写一个M文件:
function stats = reglm(y,X,model,varnames)
% 多重线性回归分析或广义线性回归分析
%
%  reglm(y,X),产生线性回归分析的方差分析表和参数估计结果,并以表格形式显示在屏幕上. 参
%  数X是自变量观测值矩阵,它是n行p列的矩阵. y是因变量观测值向量,它是n行1列的列向量.
%
%  stats = reglm(y,X),还返回一个包括了回归分析的所有诊断统计量的结构体变量stats.
%
%  stats = reglm(y,X,model),用可选的model参数来控制回归模型的类型. model是一个字符串,
%  其可用的字符串如下
%      'linear'          带有常数项的线性模型(默认情况)
%      'interaction'    带有常数项、线性项和交叉项的模型
%      'quadratic'      带有常数项、线性项、交叉项和平方项的模型
%      'purequadratic'  带有常数项、线性项和平方项的模型
%
%  stats = reglm(y,X,model,varnames),用可选的varnames参数指定变量标签. varnames 
%  可以是字符矩阵或字符串元胞数组,它的每行的字符或每个元胞的字符串是一个变量的标签,它的行
%  数或元胞数应与X的列数相同. 默认情况下,用X1,X2,…作为变量标签.
%
%  例:
%  x = [215 250 180 250 180 215 180 215 250 215 215
%      136.5 136.5 136.5 138.5 139.5 138.5 140.5 140.5 140.5 138.5 138.5]';
%  y = [6.2 7.5 4.8 5.1 4.6 4.6 2.8 3.1 4.3 4.9 4.1]';
%  reglm(y,x,'quadratic')
%
%  ----------------------------------方差分析表----------------------------------
%  方差来源    自由度          平方和          均方            F值          p值
%  回归      5.0000        15.0277        3.0055          7.6122      0.0219
%  残差      5.0000        1.9742        0.3948
%  总计      10.0000        17.0018
%
matlab拟合数据%        均方根误差(Root MSE)    0.6284            判定系数(R-Square)    0.8839
%  因变量均值(Dependent Mean)    4.7273        调整的判定系数(Adj R-Sq)    0.7678
%
%  -----------------------------------参数估计-----------------------------------
%      变量              估计值            标准误            t值          p值
%    常数项            30.9428        2011.1117          0.0154      0.9883
%        X1              0.7040          0.6405          1.0992      0.3218
%        X2            -0.8487          29.1537          -0.0291      0.9779
%    X1*X2            -0.0058          0.0044          -1.3132      0.2461
%    X1*X1              0.0003          0.0003          0.8384      0.4400
%    X2*X2              0.0052          0.1055          0.0492      0.9626
%
%  Copyright 2009 - 2010 xiezhh. 
%  $Revision: 1.0.0.0 $  $Date: 2009/12/22 21:41:00 $
if nargin < 2
  error('至少需要两个输入参数');
end
p = size(X,2);    % X的列数,即变量个数
if nargin < 3 || isempty(model)
  model = 'linear';    % model参数的默认值
end
% 生成变量标签varnames
if nargin < 4 || isempty(varnames)
    varname1 = strcat({'X'},num2str([1:p]'));
    varnames = makevarnames(varname1,model);    % 默认的变量标签
else
    if ischar(varnames)
        varname1 = cellstr(varnames);
    elseif iscell(varnames)
        varname1 = varnames(:);
    else
        error('varnames 必须是字符矩阵或字符串元胞数组');
    end
    if size(varname1,1) ~= p
        error('变量标签数与X的列数不一致');
    else
        varnames = makevarnames(varname1,model);    % 指定的变量标签
    end
end
ST = regstats(y,X,model);    % 调用regstats函数进行线性回归分析,返回结构体变量ST
f = ST.fstat;    % F检验相关结果
t = ST.tstat;    % t检验相关结果
% 显示方差分析表
fprintf('\n');
fprintf('------------------------------方差分析表------------------------------');
fprintf('\n');
fprintf('%s%7s%15s%15s%15s%12s','方差来源','自由度','平方和','均方','F值','p值');
fprintf('\n');
fmt = '%s%13.4f%17.4f%17.4f%16.4f%12.4f';
fprintf(fmt,'回归',f.dfr,f.ssr,f.ssr/f.dfr,f.f,f.pval);
fprintf('\n');
fmt = '%s%13.4f%17.4f%17.4f';
fprintf(fmt,'残差',f.dfe,f.sse,f.sse/f.dfe);
fprintf('\n');
fmt = '%s%13.4f%17.4f';
fprintf(fmt,'总计',f.dfe+f.dfr,f.sse+f.ssr);
fprintf('\n');
fprintf('\n');

版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198,我们将在24小时内删除。