[R] Using optim() function to find MLE

[Christofer Bogaso]

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