Understanding the Jeffrey prior penalty term in the glmfit function in matlab R2024a
72 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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');
%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
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:
neither are using matlab.
0 comentarios
Respuestas (0)
Ver también
Categorías
Más información sobre Linear and Nonlinear Regression en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!