gnu-octave / statistics

The Statistics package for GNU Octave
GNU General Public License v3.0
24 stars 22 forks source link

Updated `glmfit` #141

Closed ruchikasonagote closed 5 months ago

ruchikasonagote commented 5 months ago

Updates

Issue

(https://github.com/gnu-octave/statistics/issues/126) Fixing glmfit function.

TestCase

rng('default') % For reproducibility
rndvars = randn(100,2);
X = [2 + rndvars(:,1),rndvars(:,2)];
mu = exp(1 + X*[1;2]);
y = poissrnd(mu);
[b,dev] = glmfit(X,y,'poisson')

ScreenShot of output

image

pr0m1th3as commented 5 months ago

Please, stop opening new PRs when you receive a review suggesting modifications are required. This latest PR does not even work. Always test your function before making a PR. The file is named glmfit and your function inside this file is named GLM !! A lot of code is broken. Help docstring is messed up as well compared to earlier versions.

ruchikasonagote commented 5 months ago

Sorry for the mistakes. I have made some changes in the code. I actually have a doubt regarding the deviance of the binomial distribution. For the below example,

x = [2100 2300 2500 2700 2900 3100 ...
     3300 3500 3700 3900 4100 4300]';
n = [48 42 31 34 31 21 23 23 21 16 17 21]';
y = [1 2 0 3 8 8 14 17 19 15 17 21]';
[b, dev] = glmfit (x,[y n],'binomial','Link','probit')

In matlab, deviance = 7.5693 But I am getting deviance = 222.57 could you review the below code and give suggestion.

case "binomial"
          linear_predictor = X * b;
          p = exp (linear_predictor) ./ (1 + exp (linear_predictor));  
          epsilon = 1e-6; 
          if size (y, 2) == 1
            trials = ones(size (y));
          elseif size (y, 2) == 2
            trials = y(:, 2);
          endif
          dev = ...
            -2 * sum (y(:,1) .* log (max (p, epsilon)) + (trials - y(:,1)) .* log (max (1 - p, epsilon)));
ruchikasonagote commented 5 months ago

I got the issue of the deviance and resolved it. I have committed the changes. Please review the changes and let me know your feedback.

pr0m1th3as commented 5 months ago

Please, pay attention in testing the function before making PR. Once again I had to make a large amount of corrections so that it doesn't break the input validation. Pay attention to code style for the statistics package as well. Thanks