【灵敏性】具有二进制广义线性模型的顺序设计研究（matlab代码实现）

 👨‍🎓个人主页：研学社的博客 

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

---

💥1 概述

广义线性模型
的顺序实验设计 实验设计是关于选择进行测量的位置。例如，选择不同的药物剂量来检查治疗是否成功。关于统计线性模型的实验设计已经写了很多。但是，这些模型通常不能很好地描述问题。常见示例是响应为二进制（成功或失败）或响应为离散计数数据（拟合泊松模型）。通过广义线性模型（GLM）熟悉对此类数据的分析。
与要求研究人员事先确定观察数据的因子设置的单阶段实验计划不同，顺序实验设计允许根据已经观察到的数据更新和改进实验计划。

📚2 运行结果

 

 部分代码：

%% Input Parameters

% ScaleMin & ScaleMax designate the physically allowed values
% for the experiment (the design region)
ScaleMin=1;
ScaleMax=50;

ScaleX=[ScaleMin ScaleMax]; % Scaling from experiment grid to [-1,1]
RangeX=ScaleX(2)-ScaleX(1);

% Prior Distribution
NN=10000; % Number of parameter vectors to represent the prior using a discretization
MedianMean=17;
MuS=lognrnd(log(MedianMean),0.5,NN,1);
MedianSigma=0.7;
SigmaS=lognrnd(log(MedianSigma),1,NN,1);

family='binomial';
link='logit';
model='linear';

%  -----------------------------------------------------------
%              True Model parameters
 TMu=22;
 TSigma=0.3;
 Tbeta=[-TMu/TSigma 1/TSigma];
 % -----------------------------------------------------------

L=40; % Length of the experiment (total number of observations)

% more initialization
N1=500; % N1 - number of vectors to calculate the weighted sum of determinants of information matrices
RealMin=realmin;
clo=clock;
beta=[-MuS./SigmaS  1./SigmaS];
X0=[]; % This is the matrix with the experimental steps already done
RResult=[];
ss=0;
xrange0=quantile(MuS,[0.1 0.9]);% Basic range for plotting the figures

fprintf('Process began at %g:%g:%1.0f \n',clo(4:6));

%% plot prepartion
%---------------------------------
figure(1)
plot([xrange0(1) xrange0(1) xrange0(2)],[1 0 0],'w'); % empty plot - preparation only
scrsz = get(0,'ScreenSize');
set(gcf,'Position',[2*scrsz(3)/3 2*scrsz(4)/3 scrsz(3)/3 scrsz(4)/3]); % small figure upper right of the screen

figure(2)
plot([xrange0(1) xrange0(1) xrange0(2)],[0 0 1],'w'); % empty plot - preparation only
scrsz = get(0,'ScreenSize');
set(gcf,'Position',[2*scrsz(3)/3 scrsz(4)/6 scrsz(3)/3 scrsz(4)/3]); % small figure upper right of the screen

figure(3)
plot([xrange0(1) xrange0(1) xrange0(2)],[1 0 0],'w'); % empty plot - preparation only
set(gcf,'Position',[0 2*scrsz(4)/3 scrsz(3)/3 scrsz(4)/3]); % small figure upper right of the screen

figure(1);

PMSS=500; % Parameter Space plot Sampe Size - how many vectors of coefficient values to plot
CSS=150; % Curves plot Sample Size - How many curves to plot

Weights=ones(NN,1)/NN;

for Rem=L:-1:1
    fprintf('%g ',ss+1); % the number of the current observation

    % find the Weighted Median (the "best guess" for the parameter vector)
    Mu1=sortrows([MuS Weights]);
    Mu1(:,2)=cumsum(Mu1(:,2)); % Cumulative sum of weightes based on likelihood ratios
    Mu1a=abs( Mu1(:,2) - 0.5 );
    MedMu=Mu1(Mu1a==min(Mu1a),1);
    MedMu=median(MedMu); % in case there was more than one Mu that met the criterion

    Sigma1=sortrows([SigmaS Weights]);
    Sigma1(:,2)=cumsum(Sigma1(:,2)); % Cumulative sum of weightes based on likelihood ratios
    Sigma1a=abs( Sigma1(:,2) - 0.5 );
    MedSigma=Sigma1(Sigma1a==min(Sigma1a),1);
    MedSigma=median(MedSigma); % in case there was more than one Mu that met the criterion

    betaMed=[-MedMu/MedSigma 1/MedSigma];

    %% Now, use this "best guess" to choose the next point
    betaMedX(1,1)=betaMed(1)+betaMed(2)*(ScaleX(1)+0.5*(RangeX)); % Scale
    betaMedX(1,2)=betaMed(2)*ScaleX(2)-betaMed(2)*(ScaleX(1)+0.5*RangeX);
    if betaMedX(2)<25 % if not very sharp curve
        Y0=-1+(X0-ScaleX(1))./RangeX.*2;
        X0t=GLMDAUG(Y0,betaMedX',1,ss+2); % Augment 2 points, D-optimaly
        X0t=ScaleX(1)+(X0t+1)./2.*RangeX; % Scale Back
    else % if very sharp curve - rescale
        ScaleX2=[max([MedMu-4*MedSigma ScaleX(1)]) min([MedMu+4*MedSigma ScaleX(2)])];
        RangeX2=ScaleX2(2)-ScaleX2(1);
        betaMedX2(1,1)=betaMed(1)+betaMed(2)*(ScaleX2(1)+0.5*(RangeX2));
        betaMedX2(1,2)=betaMed(2)*ScaleX2(2)-betaMed(2)*(ScaleX2(1)+0.5*RangeX2);
        Y0=-1+(X0-ScaleX2(1))./RangeX2.*2;
        X0t=GLMDAUG(Y0,betaMedX2',1,ss+2);% Augment 2 points, D-optimaly
        X0t=ScaleX2(1)+(X0t+1)./2.*RangeX2;
    end;

    % X0t has the candidates points for the choice of the next point
    X0t=[median(X0t) ;X0t ]; % Add the median of the 2 points augmentation
    BestSoFar=-inf;
    for i=1:3 % there are 3 candidates
        ThisOne=0;
        X0ti=[X0; X0t(i)];
        if ss==0, X0ti=[X0ti ; X0t(2:3,:)]; end;
        F0ti=x2fx(X0ti);
        for j=1:N1   % Calculate each candidates contribution to the posterior
            ThisOnej=max(det(InfoMtrxGLM(F0ti,beta(j,:)',family,link) ) ,RealMin);
            ThisOne=ThisOne + Weights(j)*log( ThisOnej );
        end;
        if ThisOne>BestSoFar
            BestSoFar=ThisOne;

🎉3 参考文献

部分理论来源于网络，如有侵权请联系删除。

[1]徐尤南,陈洁,刘志强.TBM推进机构静刚度建模与灵敏性分析[J].机械强度,2022,44(04):901-910.DOI:10.16579/j.issn.1001.9669.2022.04.020.

[2]Hovav A. Dror and David M. Steinberg (2006). ?Sequential Experimental Design for Generalized Linear Models,? Technical Report RP-SOR-0607, Tel Aviv University. 

🌈4 Matlab代码实现