使用最小二乘法拟合由分段函数构成的广义贝塔分布的参数

本文记录了用MATLAB中的lsqcurvefit函数拟合广义贝塔分布中的α和β参数的过程
1

广义贝塔分布的CDF（累积分布函数）如下：

其中l和r为偏离中心的程度，正规贝塔分布中没有l和r参数

而
即

代码如下

func.m文件
1

function [y]=func(beta0,xdata)
a=beta0(1);
b=beta0(2);

l=-0.06;
r=0.05;

%分段函数①
sub1=xdata<=l;
y(sub1)=0.0;

%分段函数②
sub2=xdata>=r;
y(sub2)=1.0;

%分段函数③
%求积分
syms t;
f=((t-l).^(a-1)).*((r-t).^(b-1))./((r-l).^(a+b-1));

sub3=find((xdata>l)&(xdata<r));
for i=sub3(1):sub3(length(sub3))
y(i)=int(f,t,l,xdata(i))*gamma(a+b)/gamma(a)*gamma(b);
end

y=y';%得到行向量，转置成列向量，因为要与xdata和ydata的维度相同，否则下面的lsqcurvefit会报错
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

main.m文件
1

clear all;clc;close all;
data=xlsread("数据.xlsx");
xdata=data(:,1);%从excel表格中读取第一列数据作为xdata
ydata=data(:,2);%从excel表格中读取第二列数据作为ydata
beta0=[6 4] ;%α和β的初始值，随意给定
disp(xdata');
disp(ydata');

[p]=lsqcurvefit(@(beta0,xdata)func(beta0,xdata),beta0,xdata,ydata);%使用最小二乘法进行拟合
%p中保存的就是拟合得到的α和β的值
disp(p');
yy=func(p,xdata);%把拟合得到的α和β代入算
figure;
plot(xdata,ydata,'ko');%画实际CDF
hold on;
plot(xdata,yy,'-');%用拟合值画CDF
grid on;
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

程序运行后会在底部命令行窗口输出拟合得到的α和β的值

拟合结果如图