%%%%%%%%%%%% EM Yangben
n=60;
Beta0=1.2;
Beta1=1.5;
X=mvnrnd([0],[2],n);
X=[ones(n,1) X];
Beta=[Beta0;Beta1];
k=0.3;       %can shu k de qu zhi
sigma=0.2;
za=-5;zb=15;%抽样样本的取值范围
[email=fmax=@(z)(-(k*(1+exp(z))^(-k-1)*exp(z]fmax=@(z)(-(k*(1+exp(z))^(-k-1)*exp(z[/email])));
Maxz=fminunc(fmax,za,zb);
Fmax=-fmax(Maxz);
[email=fmin=@(z)((k*(1+exp(z))^(-k-1)*exp(z]fmin=@(z)((k*(1+exp(z))^(-k-1)*exp(z[/email])));
Minz=fminbnd(fmin,za,zb);
Fmin=fmin(Minz);
i=1;
YangbenZ=zeros(1,n);
while i<=n
    z=(zb-za)*unifrnd(-0.5,0.5);
    fz=(Fmax-Fmin)*unifrnd(0,1);
    Fz=k*(1+exp(z))^(-k-1)*exp(z);
    if fz<Fz
        YangbenZ(i)=z;
        i=i+1;
    else
    end
end
YangbenZ;
Y=X*Beta+sigma*YangbenZ';
data=[Y X(:,2)];
% % data=[0.8529   -0.3224
% %     2.0939   -0.3297
% %     2.8841   -0.3662
% %     3.4096    1.8542
% %     2.5699    0.5504
% %     0.4880   -0.2221
% %     3.0200    1.2764
% %     1.6828    0.1728
% %     0.1656   -1.4509
% %     4.5815    0.8769
% %     0.4294   -0.9333
% %     4.0996    1.8538
% %     0.0889   -0.7251
% %     2.2531   -0.8967
% %     2.5864    0.2264
% %     3.7090    1.0685
% %     0.3846   -1.4019
% %     4.0931    1.2184
% %    -2.5191   -3.0063
% %     3.2405    0.2036
% %     0.3358   -1.3538
% %    -2.0801   -2.1523
% %     4.3346    1.7515
% %     1.2793    0.2926
% %     3.3010    0.8752
% %     3.1298    0.8364
% %    -1.8348   -2.5560
% %     0.1707   -2.2626
% %     2.5584    0.3076
% %     2.9659    0.0945
% %     0.9508   -0.8365
% %    10.8390    4.1479
% %     3.1038    1.1372
% %    -0.6353   -2.0695
% %     5.6398    1.8221
% %    -0.0174   -0.6673
% %     4.1754    0.8394
% %     1.7335    0.1291
% %     1.7241    0.5321
% %     6.2130    1.9791
% %     1.1161   -0.8185
% %     1.4039   -1.8846
% %    -0.3797   -1.4083
% %    -0.0292   -1.7725
% %     2.6802    0.9474
% %    -0.8584   -2.0100
% %     3.0039    0.3591
% %    -1.6624   -1.7871
% %     3.6846    1.5032
% %     4.0515    1.7746
% %     5.3199    1.8271
% %     6.9899    2.4224
% %     3.5667    0.6344
% %     9.3125    4.4110
% %     4.1555    0.7344
% %     4.9533    1.3298
% %    -0.4664   -1.1732
% %     5.1991    0.9892
% %    -0.0627   -0.0273
% %     0.0598   -0.8029];
indict=data(:,1);
S=6;
data1=data(find(indict<S),:);
data2=data(find(indict>=S),:);
n12=size(data2,1);
data2(:,1)=S*ones(n12,1);
proportion=n12/n;
y1=data1(:,1);
y2=data2(:,1);
x1=data1(:,2);
x2=data2(:,2);
X1=[ones(n-n12,1) x1];
X2=[ones(n12,1) x2];
Kk=10; %%%%%%% GIBBS chou yang  start 表示马尔科夫链长度
k=zeros(1,Kk);
sigma=zeros(1,Kk);
beta0=zeros(1,Kk);
beta1=zeros(1,Kk);
c1=0.01;
d1=0.01;
c2=0.01;
d2=0.01;
miu00=1.2;
sigma00=0.01;
miu01=1.5;
sigma01=0.01;
k(1)=0.3;
sigma(1)=1/gamrnd(c2,d2);
beta0(1)=normrnd(miu00,sigma00);
beta1(1)=normrnd(miu01,sigma01);
r=length(x1);
YYY=zeros(1,100);
for  g=1:100
%%%%%%%%%%%%%%%%%%%%  start sampling
i=1;
sigmaa=0.001;sigmab=10; %抽样参数sigma的范围
[email=fsigma=@(sigma)(sigma^(-c2-r-1)*exp(-d2/sigma)*exp(sum((y1-X1*[beta0(i]fsigma=@(sigma)(sigma^(-c2-r-1)*exp(-d2/sigma)*exp(sum((y1-X1*[beta0(i[/email]) beta1(i)]')/sigma))*...
      exp(1)^sum(log((1+exp((y1-X1*[beta0(i) beta1(i)]')/sigma)).^(-k(i)-1)))*...
      exp(1)^sum(log((1+exp((y2-X2*[beta0(i) beta1(i)]')/sigma)).^(-k(i)))));
sigma0=0.1:0.001:1;
F=zeros(1,length(sigma0));
for l=1:length(sigma0)
    F(l)=fsigma(sigma0(l));
end
Fmax=max(F);
Fmin=min(F);
Yangben=zeros(1,1);
%%%%%%%%%%%%%%%%%%%%%
SUM=100;
while SUM>=0
sigma=(sigmab-sigmaa)*rand;
fsigma=(Fmax-Fmin)*rand;
Fsigma=sigma^(-c2-r-1)*exp(-d2/sigma)*exp(sum((y1-X1*[beta0(i) beta1(i)]')/sigma))*...
       exp(1)^sum(log((1+exp((y1-X1*[beta0(i) beta1(i)]')/sigma)).^(-k(i)-1)))*...
       exp(1)^sum(log((1+exp((y2-X2*[beta0(i) beta1(i)]')/sigma)).^(-k(i))));
   if fsigma<Fsigma
       Yangben=sigma;
   else
   end
   SUM=-abs(Yangben);
end
YYY(g)=Yangben;
end
YYY
plot(YYY)