PostgreSQL 聚合函数讲解 - 5 线性回归

1801年，意大利天文学家朱赛普·皮亚齐发现了第一颗小行星谷神星。经过40天的跟踪观测后，由于谷神星运行至太阳背后，使得皮亚齐失去了谷神星的位置。随后全世界的科学家利用皮亚齐的观测数据开始寻找谷神星，但是根据大多数人计算的结果来寻找谷神星都没有结果。时年24岁的高斯也计算了谷神星的轨道。奥地利天文学家海因里希·奥尔伯斯根据高斯计算出来的轨道重新发现了谷神星。

高斯使用的最小二乘法的方法发表于1809年他的著作《天体运动论》中。

法国科学家勒让德于1806年独立发现“最小二乘法”，但因不为世人所知而默默无闻。

勒让德曾与高斯为谁最早创立最小二乘法原理发生争执。

1829年，高斯提供了最小二乘法的优化效果强于其他方法的证明，因此被称为高斯-马尔可夫定理.

研究一个或多个随机变量Y1 ，Y2 ，…，Yi与另一些变量X1、X2，…，Xk之间的关系的统计方法，又称多重回归分析。通常称Y1，Y2，…，Yi为因变量，X1、X2，…，Xk为自变量。回归分析是一类数学模型，特别当因变量和自变量为线性关系时，它是一种特殊的线性模型。最简单的情形是一个自变量和一个因变量，且它们大体上有线性关系，这叫一元线性回归，即模型为Y=a+bX+ε，这里X是自变量，Y是因变量，ε是随机误差，通常假定随机误差的均值为0，方差为σ^2（σ^2大于0）σ^2与X的值无关。若进一步假定随机误差遵从正态分布，就叫做正态线性模型。一般的情形，它有k个自变量和一个因变量，因变量的值可以分解为两部分：一部分是由于自变量的影响，即表示为自变量的函数，其中函数形式已知，但含一些未知参数；另一部分是由于其他未被考虑的因素和随机性的影响，即随机误差。当函数形式为未知参数的线性函数时，称线性回归分析模型；当函数形式为未知参数的非线性函数时，称为非线性回归分析模型。当自变量的个数大于1时称为多元回归，当因变量个数大于1时称为多重回归。

回归分析的主要内容为：

①从一组数据出发，确定某些变量之间的定量关系式，即建立数学模型并估计其中的未知参数。估计参数的常用方法是最小二乘法。

②对这些关系式的可信程度进行检验。

③在许多自变量共同影响着一个因变量的关系中，判断哪个（或哪些）自变量的影响是显著的，哪些自变量的影响是不显著的，将影响显著的自变量入模型中，而剔除影响不显著的变量，通常用逐步回归、向前回归和向后回归等方法。

④利用所求的关系式对某一生产过程进行预测或控制。回归分析的应用是非常广泛的，统计软件包使各种回归方法计算十分方便。

在回归分析中，把变量分为两类。一类是因变量，它们通常是实际问题中所关心的一类指标，通常用Y表示；而影响因变量取值的的另一类变量称为自变量，用X来表示。

回归分析研究的主要问题是：

（1）确定Y与X间的定量关系表达式，这种表达式称为回归方程；

（2）对求得的回归方程的可信度进行检验；

（3）判断自变量X对因变量Y有无影响；

（4）利用所求得的回归方程进行预测和控制。

例如，如果要研究质量和用户满意度之间的因果关系，从实践意义上讲，产品质量会影响用户的满意情况，因此设用户满意度为因变量，记为Y；质量为自变量，记为X。根据图8－3的散点图，可以建立下面的线性关系： Y=A+BX+§

式中：A和B为待定参数，A为回归直线的截距；B为回归直线的斜率，表示X变化一个单位时，Y的平均变化情况；§为依赖于用户满意度的随机误差项。

对于经验回归方程： y=0.857+0.836x

回归直线在y轴上的截距为0.857、斜率0.836，即质量每提高一分，用户满意度平均上升0.836分；或者说质量每提高1分对用户满意度的贡献是0.836分。

因变量 | 自变量 

数据大概是这样的

digoal=> select count(*),regr_r2(t.cnt,test.cnt) from 

(select rn-1 as rn,cnt from test) as t,test where t.rn=test.rn and test.rn>2;

 count |      regr_r2      

-------+-------------------

   362 | 0.835282212765017

(1 row)

digoal=> select count(*),regr_r2(t.cnt,test.cnt) from 

(select rn-1 as rn,cnt from test) as t,test where t.rn=test.rn and test.rn>1;

 count |     regr_r2      

-------+------------------

   363 | 0.32915622582628

(1 row)

斜率

digoal=> select 44248*0.906861594725424+6274.25023499543;

       ?column?        

-----------------------

 46401.062078405991152

(1 row)

digoal=> select 40145/46401.062078405991152;

        ?column?        

------------------------

 0.86517416200873079820

(1 row)

截距

digoal=> select count(*),regr_intercept(t.cnt,test.cnt) from 

(select rn-1 as rn,cnt from test) as t,test where t.rn=test.rn and test.rn>1;

 count |  regr_intercept  

-------+------------------

   363 | 49279.0342891155

(1 row)

斜率

digoal=> select count(*),regr_slope(t.cnt,test.cnt) from     

(select rn-1 as rn,cnt from test) as t,test where t.rn=test.rn and test.rn>1;

 count |    regr_slope     

-------+-------------------

   363 | 0.292250474909646

(1 row)

预测结果

digoal=> select 44248*0.292250474909646+49279.0342891155;

       ?column?        

-----------------------

 62210.533302917516208

(1 row)

准确度

digoal=> select 40145/62210.533302917516208;

        ?column?        

------------------------

 0.64530872616900233730

(1 row)

regr_sxx(y, x)  :  sum(X^2) - sum(X)^2/N 

regr_sxy(y, x)  :  sum(X*Y) - sum(X) * sum(Y)/N

regr_syy(y, x)  :  sum(Y^2) - sum(Y)^2/N 

postgres=# select regr_sxx(y,x), REGR_COUNT(y,x)*VAR_POP(x) from (values(2,400),(6,401),(7,400),(3,400),(1000,488)) as t(x,y);

 regr_sxx |      ?column?       

----------+---------------------

 792833.2 | 792833.200000000000

(1 row)

postgres=# select regr_sxy(y,x), REGR_COUNT(y,x)*COVAR_POP(x,y) from (values(2,400),(6,401),(7,400),(3,400),(1000,488)) as t(x,y);

 regr_sxy | ?column? 

----------+----------

  69885.6 |  69885.6

(1 row)

postgres=# select regr_syy(y,x), REGR_COUNT(y,x)*VAR_POP(y) from (values(2,400),(6,401),(7,400),(3,400),(1000,488)) as t(x,y);

 regr_syy |       ?column?        

----------+-----------------------

   6160.8 | 6160.8000000000000000

(1 row)

=> select * from test order by to_char desc limit 10;
to_char | count
------------+--------
2015030123 | 149496
2015030122 | 165320
2015030121 | 167663
2015030120 | 161071
2015030119 | 145570
2015030118 | 133155
2015030117 | 133962
2015030116 | 130484
2015030115 | 126182
2015030114 | 122998
(10 rows)

=> do language plpgsql $$
declare
r2_1 numeric := 0;
r2_2 numeric := 0;
var int;
inter numeric;
slope numeric;
inter_2 numeric;
slope_2 numeric;
realv numeric;
predicv numeric;
offset_var int := 0; -- 最后一个值的预测值
begin
for i in 1..450 loop
with t1 as (select row_number() over(order by to_char) as rn,count from test order by to_char desc offset offset_var),
t2 as (select row_number() over(order by to_char)-1 as rn,count from test order by to_char desc offset offset_var)
select regr_intercept(t2.count,t1.count),regr_slope(t2.count,t1.count),regr_r2(t1.count,t2.count) into inter,slope,r2_1 from t1,t2 where t1.rn=t2.rn and t1.rn>i;
if r2_1>r2_2 then
inter_2:=inter;
slope_2:=slope;
r2_2:=r2_1;
var:=i;
end if;
end loop;

raise notice '%, %, %, %', var, inter_2,slope_2,r2_2;
select slope_2*count+inter_2 into predicv from test order by to_char desc offset offset_var+1 limit 1;
select count into realv from test order by to_char desc offset offset_var limit 1;
raise notice '%, %', realv, predicv;
end;
$$;
NOTICE: 436, 16599.0041292694, 0.896184690654355, 0.925125327496365
NOTICE: 149496, 164756.257188247368600
DO

=> select 149496/164756.2;
?column?
------------------------
0.90737708201572990880
(1 row)

=> do language plpgsql $$
declare
r2_1 numeric := 0; -- 相关性
r2_2 numeric := 0; -- 最大相关性
var int; -- 样本数量
inter_1 numeric; -- 截距
slope_1 numeric; -- 斜率
inter_2 numeric; -- 最大相关性截距
slope_2 numeric; -- 最大相关性斜率
realv numeric; -- 真实数据
predicv numeric; -- 预测数据
offset_var int := 1; -- 倒数第二个值的预测值, 不停迭代, 最后计算所有的实际值和预测值的corr, 看看相似度如何?
begin
for i in 1..450 loop
with t1 as (select row_number() over(order by to_char) as rn,count from test order by to_char desc offset offset_var),
t2 as (select row_number() over(order by to_char)-1 as rn,count from test order by to_char desc offset offset_var)
select regr_intercept(t2.count,t1.count),regr_slope(t2.count,t1.count),regr_r2(t1.count,t2.count) into inter_1,slope_1,r2_1 from t1,t2 where t1.rn=t2.rn and t1.rn>i;
if r2_1>r2_2 then
inter_2 := inter_1;
slope_2 := slope_1;
r2_2 := r2_1;
var := i;
end if;
end loop;

raise notice '样本数量%, 截距%, 斜率%, 相关性%', var, round(inter_2,4), round(slope_2,4), round(r2_2,4);
select slope_2*count+inter_2 into predicv from test order by to_char desc offset offset_var+1 limit 1;
select count into realv from test order by to_char desc offset offset_var limit 1;
raise notice '真实数据%, 预测数据%, 本次预测偏差,%%%', realv, round(predicv), abs(1-round(predicv/realv,4))*100;
end;
$$;
NOTICE: 样本数量436, 截距10109.8500, 斜率0.9573, 相关性0.9476
NOTICE: 真实数据165320, 预测数据170611, 本次预测偏差,%3.2000
DO

=> create or replace function check_predict(IN ov int, OUT rv numeric, OUT pv numeric, OUT dev numeric) returns record as $$

declare 

  r2_1 numeric := 0; -- 相关性

  r2_2 numeric := 0; -- 最大相关性

  var int;  --  样本数量

  inter_1 numeric;  --  截距

  slope_1 numeric;  --  斜率

  inter_2 numeric;  --  最大相关性截距

  slope_2 numeric;  --  最大相关性斜率

  realv numeric;    --  真实数据

  predicv numeric;  --  预测数据

  offset_var int := ov;   -- 倒数第二个值的预测值, 不停迭代, 最后计算所有的实际值和预测值的corr, 看看相似度如何?

  lps int := 0;

begin

  select count(*)-offset_var-4 into lps from test;  --  循环不要超过总样本数, 同时至少给2个样本.

  for i in 1..lps loop 

    with t1 as (select row_number() over(order by to_char) as rn,to_char,count from test order by to_char desc offset offset_var), 

         t2 as (select row_number() over(order by to_char)-1 as rn,to_char,count from test order by to_char desc offset offset_var) 

      select regr_intercept(t2.count,t1.count),regr_slope(t2.count,t1.count),regr_r2(t1.count,t2.count) into inter_1,slope_1,r2_1 from t1,t2 where t1.rn=t2.rn and t1.rn>i;

    if r2_1>r2_2 then 

      inter_2 := inter_1;

      slope_2 := slope_1;

      r2_2 := r2_1;

      var := i;

    end if;

  end loop;

  raise notice '样本数量%, 截距%, 斜率%, 相关性%', var, round(inter_2,4), round(slope_2,4), round(r2_2,4);

  select slope_2*count+inter_2 into predicv from test order by to_char desc offset offset_var+1 limit 1;

  select count into realv from test order by to_char desc offset offset_var limit 1;

  raise notice '真实数据%, 预测数据%, 本次预测偏差%%%', realv, round(predicv), abs(1-round(predicv/realv,4))*100;

  rv := realv;

  pv := round(predicv);

  dev := abs(1-round(predicv/realv,4));

  return;

end;

$$ language plpgsql;

=> select check_predict(i) from generate_series(1,100) t(i);

预测误差 :