提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

前言

提示：这里可以添加本文要记录的大概内容：
matlab数据分析实验，围绕多元线性回顾参数估计展开，本篇记录使用变分贝叶斯法模拟参数（主要记录代码，有问题再改）

提示：以下是本篇文章正文内容，下面案例可供参考

一、原理

多元线性回归（暂时不考虑随机误差项，论文时再加）
变分贝叶斯：（暂时没理解，直接按照老师给的图片做的，后续理解了再用自己语言解释）
文献参考如下：（讲述变分贝叶斯算法中涉及到的参数，记住（2）式（3）式（4）式和logp(y)即可）
算法原理：
先假设Bq(sigma2)>0
迭代（2）式（3）式（4）式
当p(y)小于某个值时停止迭代（注意：这里的小于某个值是指相邻两次p(y)间的差）（因为是差，某些常数在下面代码中就没写了，要写的自己写，也可以问）
（下述求解中，p(y)暂时做不出来，直接对logp(y)进行比较）
||x||：矩阵的迹（主对角线上元素和）（代码：trace（x))

二、使用步骤

1.产生数据

代码如下（示例）：
要运行记得在交互页面，脚本页面运行不了

%产生xy数据
function [y,x]=gendata(n,p,beta)%做一个函数，名字我令为gendata
x=normrnd(0,1,n,p);%x是n*p矩阵
muy=x*beta';%乘，a的列数=b的行数,pinv矩阵的逆
y=normrnd(muy,1,n,1);%y是n*1矩阵
end
%真实值设置为beta数据根据真实值产生，由数据模拟结果与真实值进行对比（此例子中，beta是4维的，p=4)
%beta=[1,-sqrt(3)/3,sqrt(2),0.8];sigma=1;n=100;p=4;
%[y,x]=gendata(50,4,[1,-sqrt(3)/3,sqrt(2),0.8],1)

代码都只是函数，完整的在后面

2.变分贝叶斯算法

代码如下（示例）：（比较复杂，每一步都是慢慢摸索的解释也比较多，可以不要）

 %变分贝叶斯法
function [mbeta,msigma2]=VBmoni(err,Bq,l0,x,y,n,p)
err=log(0.0001);%精度
A=0.001;
B=0.001;
sigmab2=1000;
Bq=0.01;%迭代的初始值
l0=0;%logpymin0,迭代的初始值 %while(l0<err)%如果l0的值小于err时停止迭代

i=1;
for i=1:100 %进行100次循环以此代替迭代
    i=i+1;
    SIGMA=inv((A+n/2)/Bq*x'*x+1/sigmab2*eye(p));%1式,4*4,p*p
    mu=((A+n/2)/Bq)*SIGMA*x'*y;%4*1,p*1,2式
    Bq=B+1/2*(y-x*mu)'*(y-x*mu)+1/2*trace(x'*x*SIGMA);%1*1,3式
    l1=1/2*log(det(SIGMA))-(mu'*mu+trace(SIGMA))/(2*sigmab2)-(A+n/2)*log(Bq);%logpymin1
    
    l0=l1-l0;%logpymin0=logpymin1-logpymin0，py值的差  
   
    %VBmoni(err,Bq,l0,x,y,n,p);%可以用这个来迭代就不要for 但是我没弄出来，求指教
end
mbeta=mu';%matlab出来结果都是列，转成行好看
msigma2=((A+n/2)/Bq)'; 
end

3.总代码

做一次模拟，没有多次，暂时只做到这里

%变分贝叶斯__模拟多元线性回归参数
%定义的数据类
n=500;
p=4;
beta=[1,-sqrt(3)/3,sqrt(2),0.8];% beta=[1.0000,-0.5774,1.4142,0.8000]
sigma=1;
%result1=zeros(100,5);
%true=[beta,sigma];
%产生数据
[y,x]=gendata(n,p,beta);
%运行变分贝叶斯
err=log(0.0001);
Bq=0.01;%迭代的初始值
l0=0;%logpymin0,迭代的初始值
[mbeta,msigma2]=VBmoni(err,Bq,l0,x,y,n,p);
result1=[mbeta,msigma2]

产生数据比较多，截屏一点点就行

总结

本文仅仅只是学习记录，以及分享，有问题的话欢迎提出，我看看我会不会，为了论文奋斗吧。
温馨提示：对于matlab矩阵的这块不熟悉的话，多问问，区分矩阵的逆和转置，我当初是弄错了，改了好久代码的
为了背后一点点摸索的心酸，最后写出来一点点东西，尽管可能不厉害，但我自豪就行，加油干吧