# bayesian logistic regression - slicesample - finding Machine learning parameters

8 ビュー (過去 30 日間)
Matthias 2016 年 7 月 27 日
コメント済み: Tom Lane 2016 年 7 月 28 日
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));
% http://stats.stackexchange.com/questions/46429/transform-data-to-desired-mean-and-standard-deviation
featAdd1s = [ones(size(featOr,1),1), scaleFe(featOr(:,1)), scaleFe(featOr(:,2)), scaleFe(featOr(:,3)), scaleFe(featOr(:,4))];
% The number of tests
% 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
%%Slice Sampling
% initial = [1 1 1 1 1];
initial = [1 1 1 1 1];
post(initial)
% 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...

サインインしてコメントする。

### 回答 (1 件)

Tom Lane 2016 年 7 月 27 日
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.
##### 2 件のコメント表示非表示 1 件の古いコメント
Tom Lane 2016 年 7 月 28 日
1. Bayes methods are not limited by the number of samples. If you use the 'logpdf' approach that I illustrated, you should be okay.
2. If you only want to get estimates and use them for prediction, you could take the mean of the trace values, possibly omitting some top rows to avoid the effects of the initial values before the traces settle down. You are right that you would have to transform the new X features using the same scaling that you used during fitting. That is, scale using the mean and std of the X from fitting, not by separately scaling new X values based on their own mean and std.

サインインしてコメントする。

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by