Understanding the Jeffrey prior penalty term in the glmfit function in matlab R2024a

Question

Moosejr 2024-9-4

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2149979-understanding-the-jeffrey-prior-penalty-term-in-the-glmfit-function-in-matlab-r2024a

编辑： Moosejr 2024-11-20，15:57

Hi,

Consider the example separable dataset

x_i = 800:5:900;
y_i = [0 0 0 0 0, ...
    0 0 0 0 0, ...
    1 1 1 1 1, ...
    1 1 1 1 1, ...
    1];

In the R2024a version of matlab it is now possible to fit a probit curve to these data using Jeffrey's prior as a penalty function.

%Standard way - singular coefficients
b_standard = glmfit(x_i', y_i', 'binomial','link', 'probit');
Warning: Iteration limit reached.
Warning: The estimated coefficients perfectly separate failures from successes. This means the theoretical best estimates are not finite. For the fitted linear combination XB of the predictors, the sample proportions P of Y=N in the data satisfy:
   XB<-0.0197012: P=0
   XB>-0.0197012: P=1
To get proper estimates, fit with 'LikelihoodPenalty' set to 'jeffreys-prior'.
%Jeffrey prior - nonsingular coefficients
b_jeff = glmfit(x_i', y_i', 'binomial','link', 'probit',LikelihoodPenalty="jeffreys-prior");

%Plot the standard and Jeffrey-prior coefficients

figure(1)

hold on

plot(b_standard(1),b_standard(2),'r*')

plot(b_jeff(1),b_jeff(2),'k*')

This is good because instead of obtaining divergent coefficients we obtain finite coefficients which are shifted towards the origin.

My questions are related to what is going on under the hood of the glmfit function.

Is there a simple way to plot the Jeffrey prior that is being used (as a function of possible coefficient values), and if so how?
Is it also possible to plot the resulting likelihood function (as a function of possible coefficient values) when we used the Jeffrey prior?

In the case where the Jeffrey prior is not being used, I can compute the likelihood function as

%Plot the likelihood

%I here transformed from b_0 and b_1 to the mean and slope parameters

%because the plots often look nicer in the latter coordinates

mu = linspace(840,860,201);

sigma = linspace(-1,5,101);

LikelihoodList = [];

for i = 1:length(mu)

for j = 1:length(sigma)

LikelihoodList = [LikelihoodList; mu(i), sigma(j), LikelihoodFunction(-mu(i)./sigma(j),1./sigma(j),x_i,y_i)];

end

X = LikelihoodList(:,1); Y = LikelihoodList(:,2); [x,y]=meshgrid(X,Y); Z = LikelihoodList(:,3);

figure(2)

plot3(X,Y,Z)

function likelihood = LikelihoodFunction(b_0,b_1,x_i,y_i)

%Computes the likelihood

%Take in a single value of b0 and b1 and a list of experimental data x_i and y_i

p_i = normcdf(b_0 + b_1*x_i);

likelihood_product_list = p_i.^y_i.*(1-p_i).^(1-y_i);

likelihood = prod(likelihood_product_list);

end

Related texts can be found here:

brglm: Bias Reduction in Binomial-Response Generalized Linear Models (r-project.org)

brglm2: Bias Reduction in Generalized Linear Models (r-project.org)

[P] Fit Logistic Regression with Jeffreys Prior : r/MachineLearning (reddit.com)

Firth's bias correction as a Bayesian model | Random effect

neither are using matlab.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Arnav 2024-11-5，8:28

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2149979-understanding-the-jeffrey-prior-penalty-term-in-the-glmfit-function-in-matlab-r2024a#answer_1541060

在 MATLAB Online 中打开

Hi Moosejr,

According to the sources mentioned, Jeffrey’s prior is proportional to the square root of the determinant of the fisher information. Fisher Information can be calculated as shown in the code snippet shown below:

function I = FisherInformation(b_0, b_1, x_i)
    I = zeros(2, 2); 
    for i = 1:length(x_i)  
        eta_i = b_0 + b_1 * x_i(i);  
        p_i = normcdf(eta_i);   
        phi_i = normpdf(eta_i);  
        % Clip probabilities to avoid division by zero    
        p_i = max(min(p_i, 1 - 1e-10), 1e-10);    
        w = (phi_i^2) / (p_i * (1 - p_i));   
        I(1, 1) = I(1, 1) + w;    
        I(1, 2) = I(1, 2) + w * x_i(i);    
        I(2, 1) = I(1, 2);    
        I(2, 2) = I(2, 2) + w * x_i(i)^2;  
    end
    % Regularization term to avoid singular matrix
    I = I + 1e-10 * eye(2);
end

Then, the likelihood with Jeffrey’s prior can be found as shown in the code snippet below:

likelihood = LikelihoodFunction(b_0, b_1, x_i, y_i);
fisher_info = FisherInformation(b_0, b_1, x_i);
jeffrey_prior = sqrt(det(fisher_info));
LikelihoodList(i, j) = likelihood * jeffrey_prior;

The plot of the likelihood is shown below:

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Moosejr 2024-11-20，15:04

编辑：Moosejr 2024-11-20，15:57

在 MATLAB Online 中打开

Hi @Arnav,

Thank you for your excellent response! I'm very sorry for replying so late, I have been travelling without access to my computer.

I have tried to implement your suggestion in a minimal example, I hope I got it right. I think we get the same curve (I'm not 100% sure because it is a bit hard to compare your png/jpg to my matlab file since I can't locally rotate the former).

However, I noticed that the maximum likelihood value differs from the ML estimate provided by matlab's glmfit function (the 'ko'-point differs from the maximum value of the likelihood function that I plotted). I therefore assume that something is amiss. Do you know what is wrong and how it can be fixed?

x_i = 800:5:900;

y_i = [0 0 0 0 0, ...

0 0 0 0 0, ...

1 1 1 1 1, ...

1];

%Jeffrey prior - nonsingular coefficients

b_jeff = glmfit(x_i', y_i', 'binomial','link', 'probit', LikelihoodPenalty= "jeffreys-prior");

LikelihoodList = [];

mu = linspace(830,870,50);

sigma = linspace(-1,15,50);

[X,Y] = meshgrid(sigma, mu);

for i = 1:length(mu)

for j = 1:length(sigma)

likelihood = LikelihoodFunction(-mu(i)./sigma(j), 1./sigma(j), x_i, y_i);

fisher_info = FisherInformation(-mu(i)./sigma(j), 1./sigma(j), x_i);

jeffrey_prior = sqrt(det(fisher_info));

LikelihoodList(i, j) = likelihood * jeffrey_prior;

end

figure(3)

surf(X, Y, LikelihoodList)

hold on

plot3(1/b_jeff(2), -b_jeff(1)/b_jeff(2), max(LikelihoodList(:)), 'ko');

%% Functions

function likelihood = LikelihoodFunction(b_0,b_1,x_i,y_i)

%Computes the likelihood

%Take in a single value of b0 and b1 and a list of experimental data x_i and y_i

p_i = normcdf(b_0 + b_1*x_i);

likelihood_product_list = p_i.^y_i.*(1-p_i).^(1-y_i);

likelihood = prod(likelihood_product_list);

end

function I = FisherInformation(b_0, b_1, x_i)

I = zeros(2, 2);

for i = 1:length(x_i)

eta_i = b_0 + b_1 * x_i(i);

p_i = normcdf(eta_i);

phi_i = normpdf(eta_i);

% Clip probabilities to avoid division by zero

p_i = max(min(p_i, 1 - 1e-10), 1e-10);

w = (phi_i^2) / (p_i * (1 - p_i));

I(1, 1) = I(1, 1) + w;

I(1, 2) = I(1, 2) + w * x_i(i);

I(2, 1) = I(1, 2);

I(2, 2) = I(2, 2) + w * x_i(i)^2;

end

% Regularization term to avoid singular matrix

I = I + 1e-10 * eye(2);

end

请先登录，再进行评论。

Understanding the Jeffrey prior penalty term in the glmfit function in matlab R2024a

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Understanding the Jeffrey prior penalty term in the glmfit function in matlab R2024a

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论