bayesian logistic regression - slicesample - finding Machine learning parameters
10 views (last 30 days)
since I have problems with separation for logistic regression I would like to use bayesian logistic regression
I follow this script bayesian logistic regression
However it is for 1D and my problem has 4 features, not 1. I run into problems where the slicesample function is used, and I can't solve the error message "INITIAL is not in the domain of the target distribution."
I tried different initial, but somehow my posterior distribution "post" is always zero. Did I define it wrong? Shall I use another distribution than the binomial one? Or maybe I made something wrong for "logitp" and its use in "post". I feed in a matrix "featAdd1s", while in the link mentioned above, a vector "weight" is used.
Any help highly appreciated. I attached the GT.mat needed for the code.
% features = predictors
featOr = GT(:,1:end-1); % original feature
% recenter and rescale input features as suggested by Gelman2008
% "A weakly informative default prior distribution for logistic and other regression models"
% such that each input feature has mean 0 and std = 0.5
scaleFe = @(feat) (feat.*(0.5/std(feat)))-mean(feat)*(0.5/std(feat));
% add ones, necessary?
featAdd1s = [ones(size(featOr,1),1), scaleFe(featOr(:,1)), scaleFe(featOr(:,2)), scaleFe(featOr(:,3)), scaleFe(featOr(:,4))];
% The number of tests
total = ones(size(featAdd1s,1),1);
% The number of success
poor = GT(:,end);
%%Logistic Regression Model
logitp = @(b,X) 1./(1+exp(-X*b));
%%use priors according to Gelman2008 - "A weakly informative default prior distribution for logistic and other regression models"
% for intercept
pdIn = makedist('tLocationScale','mu',0,'sigma',10,'nu',1);
prior1 = @(b0) pdf(pdIn,b0);
pd = makedist('tLocationScale','mu',0,'sigma',2.5,'nu',1); % for predictors
prior2 = @(b1) pdf(pd,b1);
prior3 = @(b2) pdf(pd,b2);
prior4 = @(b3) pdf(pd,b3);
prior5 = @(b4) pdf(pd,b4);
By Bayes' theorem, the joint posterior distribution of the model parameters is proportional to the product of the likelihood and priors.
post = @(b) prod(binopdf(poor,total,logitp(b',featAdd1s)))*prior1(b(1))*prior2(b(2))*prior3(b(3))*prior4(b(4))*prior5(b(5)) ; % priors
% initial = [1 1 1 1 1];
initial = [1 1 1 1 1];
% initial = [0.02 0.02 0.02 0.02 0.02];
nsamples = 1000;
% trace = slicesample(initial,nsamples,'pdf',post);
trace = slicesample(initial,nsamples,'pdf',post,'width',[10, 10, 10, 10,10]); % ?! No idea...10 is default...
Tom Lane on 27 Jul 2016
It looks like you have the right idea. But I suspect the calculations are underflowing. You're multiplying thousands of probability values, many perhaps quite small. I was able to make your problem work by using the log posterior instead of the posterior itself:
logpost = @(b) sum(log(binopdf(poor,total,logitp(b',featAdd1s)))) + ...
log(prior1(b(1))) + log(prior2(b(2))) + log(prior3(b(3))) + log(prior4(b(4))*prior5(b(5)))
trace = slicesample(initial,nsamples,'logpdf',logpost,'width',[10, 10, 10, 10,10]);
I may try to get the published example changed to use this technique.
By the way, it's not necessary in your problem, but sometimes setting the slope coefficients to 0 as an initial value, and the intercept coefficient to some moderate value, can give a starting point that will at least be feasible.