Hi,
I am trying to fit a GLM on below data. While R does provide direct
estimation, I wanted to go with manual calculation as below
dat = structure(list(PurchasedProb = c(0.37212389963679, 0.572853363351896,
0.908207789994776, 0.201681931037456, 0.898389684967697, 0.944675268605351,
0.660797792486846, 0.62911404389888, 0.0617862704675645, 0.205974574899301,
0.176556752528995, 0.687022846657783, 0.384103718213737, 0.769841419998556,
0.497699242085218, 0.717618508264422, 0.991906094830483, 0.380035179434344,
0.777445221319795, 0.934705231105909, 0.212142521282658, 0.651673766085878,
0.125555095961317, 0.267220668727532, 0.386114092543721, 0.0133903331588954,
0.382387957070023, 0.86969084572047, 0.34034899668768, 0.482080115471035,
0.599565825425088, 0.493541307048872, 0.186217601411045, 0.827373318606988,
0.668466738192365, 0.79423986072652, 0.107943625887856, 0.723710946040228,
0.411274429643527, 0.820946294115856, 0.647060193819925, 0.78293276228942,
0.553036311641335, 0.529719580197707, 0.789356231689453, 0.023331202333793,
0.477230065036565, 0.7323137386702, 0.692731556482613, 0.477619622135535,
0.8612094768323, 0.438097107224166, 0.244797277031466, 0.0706790471449494,
0.0994661601725966, 0.31627170718275, 0.518634263193235, 0.6620050764177,
0.406830187188461, 0.912875924259424, 0.293603372760117, 0.459065726259723,
0.332394674187526, 0.65087046707049, 0.258016780717298, 0.478545248275623,
0.766310670645908, 0.0842469143681228, 0.875321330036968, 0.339072937844321,
0.839440350187942, 0.34668348915875, 0.333774930797517, 0.476351245073602,
0.892198335845023, 0.864339470630512, 0.389989543473348, 0.777320698834956,
0.960617997217923, 0.434659484773874, 0.712514678714797, 0.399994368897751,
0.325352151878178, 0.757087148027495, 0.202692255144939, 0.711121222469956,
0.121691921027377, 0.245488513959572, 0.14330437942408, 0.239629415096715,
0.0589343772735447, 0.642288258532062, 0.876269212691113, 0.778914677444845,
0.79730882588774, 0.455274453619495, 0.410084082046524, 0.810870242770761,
0.604933290276676, 0.654723928077146, 0.353197271935642, 0.270260145887733,
0.99268406117335, 0.633493264438584, 0.213208135217428, 0.129372348077595,
0.478118034312502, 0.924074469832703, 0.59876096714288, 0.976170694921166,
0.731792511884123, 0.356726912083104, 0.431473690550774, 0.148211560677737,
0.0130775754805654, 0.715566066093743, 0.103184235747904, 0.446284348610789,
0.640101045137271, 0.991838620044291, 0.495593577856198, 0.484349524369463,
0.173442334868014, 0.754820944508538, 0.453895489219576, 0.511169783771038,
0.207545113284141, 0.228658142732456, 0.595711996313184, 0.57487219828181,
0.0770643802825361, 0.0355405795853585, 0.642795492196456, 0.928615199634805,
0.598092422354966, 0.560900748008862, 0.526027723914012, 0.985095223877579,
0.507641822332516, 0.682788078673184, 0.601541217649356, 0.238868677755818,
0.258165926672518, 0.729309623362496, 0.452570831403136, 0.175126768415794,
0.746698269620538, 0.104987640399486, 0.864544949028641, 0.614644971676171,
0.557159538846463, 0.328777319053188, 0.453131445450708, 0.500440972624347,
0.180866361130029, 0.529630602803081, 0.0752757457084954, 0.277755932649598,
0.212699519237503, 0.284790480975062, 0.895094102947041, 0.4462353233248,
0.779984889784828, 0.880619034869596, 0.413124209502712, 0.0638084805104882,
0.335487491684034, 0.723725946620107, 0.337615333497524, 0.630414122482762,
0.840614554006606, 0.856131664710119, 0.39135928102769, 0.380493885604665,
0.895445425994694, 0.644315762910992, 0.741078648716211, 0.605303446529433,
0.903081611497328, 0.293730155099183, 0.19126010988839, 0.886450943304226,
0.503339485730976, 0.877057543024421, 0.189193622441962, 0.758103052387014,
0.724498892668635, 0.943724818294868, 0.547646587016061, 0.711743867723271,
0.388905099825934, 0.100873126182705, 0.927302088588476, 0.283232500310987,
0.59057315881364, 0.110360604943708, 0.840507032116875, 0.317963684443384,
0.782851336989552, 0.267508207354695, 0.218645284883678, 0.516796836396679,
0.268950592027977, 0.181168327340856, 0.518576137488708, 0.562782935798168,
0.129156854469329, 0.256367604015395, 0.717935275984928, 0.961409936426207,
0.100140846567228, 0.763222689507529, 0.947966354666278, 0.818634688388556,
0.308292330708355, 0.649579460499808, 0.953355451114476, 0.953732650028542,
0.339979203417897, 0.262474110117182, 0.165453933179379, 0.322168056620285,
0.510125206550583, 0.923968471353874, 0.510959698352963, 0.257621260825545,
0.0464608869515359, 0.41785625834018, 0.854001502273604, 0.347230677725747,
0.131442320533097, 0.374486864544451, 0.631420228397474, 0.390078933676705,
0.689627848798409, 0.689413412474096, 0.554900623159483, 0.429624407785013,
0.452720062807202, 0.306443258887157, 0.578353944001719, 0.910370304249227,
0.142604082124308, 0.415047625312582, 0.210925750667229, 0.428750370861962,
0.132689975202084, 0.460096445865929, 0.942957059247419, 0.761973861604929,
0.932909828843549, 0.470678497571498, 0.603588067693636, 0.484989680582657,
0.10880631650798, 0.247726832982153, 0.498514530714601, 0.372866708086804,
0.934691370232031, 0.523986077867448, 0.317144671687856, 0.277966029476374,
0.787540507735685, 0.702462512534112, 0.165027638664469, 0.0644575387705117,
0.754705621628091, 0.620410033036023, 0.169576766667888, 0.0622140523046255,
0.109029268613085, 0.381716351723298, 0.169310914585367, 0.298652542056516,
0.192209535045549, 0.257170021301135, 0.181231822818518, 0.477313709678128,
0.770737042883411, 0.0277871224097908, 0.527310776989907, 0.880319068906829,
0.373063371982425, 0.047959131654352, 0.138628246728331, 0.321492120390758,
0.154831611318514, 0.132228172151372, 0.221305927727371, 0.226380796171725,
0.13141653384082, 0.981563460314646, 0.327013726811856, 0.506939497077838,
0.68144251476042, 0.0991691031958908, 0.118902558228001, 0.0504396595060825,
0.929253919748589, 0.673712232848629, 0.0948578554671258, 0.492596120806411,
0.461551840649918, 0.375216530868784, 0.991099219536409, 0.176350713707507,
0.813435208518058, 0.0684466371312737, 0.40044974675402, 0.141144325723872,
0.193309862399474, 0.841351716779172, 0.719913988374174, 0.267212083097547,
0.495001644827425, 0.0831138978246599, 0.353884240612388, 0.969208805356175,
0.624714189674705, 0.664618249749765, 0.312489656498656, 0.405689612729475,
0.996077371528372, 0.855082356370986, 0.953548396006227, 0.812305092345923,
0.782182115828618, 0.267878128914163, 0.762151529546827, 0.986311589134857,
0.293605549028143, 0.399351106956601, 0.812131523853168, 0.0771516691893339,
0.363696809858084, 0.442592467181385, 0.156714132521302, 0.582205270184204,
0.970162178855389, 0.98949983343482, 0.176452036481351, 0.542130424408242,
0.384303892031312, 0.676164050819352, 0.26929377974011, 0.469250942347571,
0.171800082316622, 0.36918946239166, 0.725405272562057, 0.486149104312062,
0.0638024667277932, 0.784546229988337, 0.418321635806933, 0.981018084799871,
0.282883955864236, 0.847882149508223, 0.0822392308618873, 0.886458750581369,
0.471930730855092, 0.109100963454694, 0.333277984522283, 0.837416569236666,
0.276849841699004, 0.587035141419619, 0.836732269730419, 0.0711540239863098,
0.702778743347153, 0.69882453721948, 0.46396238100715, 0.436931110452861,
0.562176787760109, 0.928483226336539, 0.230466414242983, 0.221813754411414,
0.420215893303975, 0.333520805696025, 0.864807550329715, 0.177194535732269,
0.49331872863695, 0.429713366786018, 0.564263842999935, 0.656162315513939,
0.97855406277813, 0.232161148451269, 0.240811596158892, 0.796836083987728,
0.831671715015545, 0.113507706439123, 0.963312016334385, 0.147322899661958,
0.143626942764968, 0.925229935208336, 0.507035602582619, 0.15485101705417,
0.348302052123472, 0.659821025794372, 0.311772374436259, 0.351573409279808,
0.147845706902444, 0.658877609064803), Age = c(19L, 35L, 26L,
27L, 19L, 27L, 27L, 32L, 25L, 35L, 26L, 26L, 20L, 32L, 18L, 29L,
47L, 45L, 46L, 48L, 45L, 47L, 48L, 45L, 46L, 47L, 49L, 47L, 29L,
31L, 31L, 27L, 21L, 28L, 27L, 35L, 33L, 30L, 26L, 27L, 27L, 33L,
35L, 30L, 28L, 23L, 25L, 27L, 30L, 31L, 24L, 18L, 29L, 35L, 27L,
24L, 23L, 28L, 22L, 32L, 27L, 25L, 23L, 32L, 59L, 24L, 24L, 23L,
22L, 31L, 25L, 24L, 20L, 33L, 32L, 34L, 18L, 22L, 28L, 26L, 30L,
39L, 20L, 35L, 30L, 31L, 24L, 28L, 26L, 35L, 22L, 30L, 26L, 29L,
29L, 35L, 35L, 28L, 35L, 28L, 27L, 28L, 32L, 33L, 19L, 21L, 26L,
27L, 26L, 38L, 39L, 37L, 38L, 37L, 42L, 40L, 35L, 36L, 40L, 41L,
36L, 37L, 40L, 35L, 41L, 39L, 42L, 26L, 30L, 26L, 31L, 33L, 30L,
21L, 28L, 23L, 20L, 30L, 28L, 19L, 19L, 18L, 35L, 30L, 34L, 24L,
27L, 41L, 29L, 20L, 26L, 41L, 31L, 36L, 40L, 31L, 46L, 29L, 26L,
32L, 32L, 25L, 37L, 35L, 33L, 18L, 22L, 35L, 29L, 29L, 21L, 34L,
26L, 34L, 34L, 23L, 35L, 25L, 24L, 31L, 26L, 31L, 32L, 33L, 33L,
31L, 20L, 33L, 35L, 28L, 24L, 19L, 29L, 19L, 28L, 34L, 30L, 20L,
26L, 35L, 35L, 49L, 39L, 41L, 58L, 47L, 55L, 52L, 40L, 46L, 48L,
52L, 59L, 35L, 47L, 60L, 49L, 40L, 46L, 59L, 41L, 35L, 37L, 60L,
35L, 37L, 36L, 56L, 40L, 42L, 35L, 39L, 40L, 49L, 38L, 46L, 40L,
37L, 46L, 53L, 42L, 38L, 50L, 56L, 41L, 51L, 35L, 57L, 41L, 35L,
44L, 37L, 48L, 37L, 50L, 52L, 41L, 40L, 58L, 45L, 35L, 36L, 55L,
35L, 48L, 42L, 40L, 37L, 47L, 40L, 43L, 59L, 60L, 39L, 57L, 57L,
38L, 49L, 52L, 50L, 59L, 35L, 37L, 52L, 48L, 37L, 37L, 48L, 41L,
37L, 39L, 49L, 55L, 37L, 35L, 36L, 42L, 43L, 45L, 46L, 58L, 48L,
37L, 37L, 40L, 42L, 51L, 47L, 36L, 38L, 42L, 39L, 38L, 49L, 39L,
39L, 54L, 35L, 45L, 36L, 52L, 53L, 41L, 48L, 48L, 41L, 41L, 42L,
36L, 47L, 38L, 48L, 42L, 40L, 57L, 36L, 58L, 35L, 38L, 39L, 53L,
35L, 38L, 47L, 47L, 41L, 53L, 54L, 39L, 38L, 38L, 37L, 42L, 37L,
36L, 60L, 54L, 41L, 40L, 42L, 43L, 53L, 47L, 42L, 42L, 59L, 58L,
46L, 38L, 54L, 60L, 60L, 39L, 59L, 37L, 46L, 46L, 42L, 41L, 58L,
42L, 48L, 44L, 49L, 57L, 56L, 49L, 39L, 47L, 48L, 48L, 47L, 45L,
60L, 39L, 46L, 51L, 50L, 36L, 49L)), class = "data.frame", row.names = c(NA,
-400L))
dat[, 'b0'] = 1
LL = function(b0, b1)
sum(apply(as.matrix(dat[, c('PurchasedProb', 'Age')]), 1, function(iROw)
iROw['PurchasedProb'] * log( 1 / (1 + exp(-1 * (b0 + b1 * iROw['Age'])))) +
(1 - iROw['PurchasedProb']) * log(1 - 1 / (1 + exp(-1 * (b0 + b1 *
iROw['Age']))))))
result1 = optim(par = c(0,1), fn = LL, method = "L-BFGS-B")
coef(result1)
Unfortunately my calculation fails. Could you please help me to
understand what went wrong in my calculation?
Additionally, I want to use Newton-CG algorithm in my optimisation as
to follow certain norm which I am bound to. Is there any way to force
optim() function to use Newton-CG algorithm?
I am expecting to get similar result as
coef(glm("PurchasedProb ~ Age", data = dat, family = binomial())) ###
0.268640013 -0.007738782
Using optim() function to find MLE
2 messages · Christofer Bogaso, Ivan Krylov
? Mon, 29 Jul 2024 09:52:22 +0530 Christofer Bogaso <bogaso.christofer at gmail.com> ?????:
LL = function(b0, b1)
help(optim) documents that the function to be optimised takes a single
argument, a vector containing the parameters. Here's how your LL
function can be adapted to this interface:
LL <- function(par) {
b0 <- par[1]
b1 <- par[2]
sum(apply(as.matrix(dat[, c('PurchasedProb', 'Age')]), 1,
function(iROw) iROw['PurchasedProb'] * log( 1 / (1 + exp(-1 * (b0 + b1
* iROw['Age'])))) + (1 - iROw['PurchasedProb']) * log(1 - 1 / (1 +
exp(-1 * (b0 + b1 * iROw['Age']))))))
}
Furethermore, LL(c(0, 1)) results in -Inf. All the methods supported by
optim() require at least the initial parameters to result in finite
values, and L-BFGS-B requires all evaluations to be finite. You're also
maximising the function, and optim() defaults to minimisation, so you
need an additional parameter to adjust that (or rewrite the LL function
further):
result1 <- optim(
par = c(0, 0), fn = LL, method = "L-BFGS-B",
control = list(fnscale = -1)
)
coef(result1)
help(optim) documents the return value of optim() as not having a class or a $coefficients field. You can use result1$par to access the parameters.
Is there any way to force optim() function to use Newton-CG algorithm?
I'm assuming you mean the method documented in <https://docs.scipy.org/doc/scipy/reference/optimize.minimize-newtoncg.html>. optim() doesn't support the truncated (line search) Newton-CG method. See the 'optimx' and 'nloptr' packages for an implementation of a truncated Newton method (not necessarily exactly the same one).
Best regards, Ivan