title: Can someone please explain dynamic time warping for determining time series similarity?

text: I am trying to grasp the dynamic time warping measure for comparing time series together. I have three time series datasets like this:               T1 <- structure(c(0.000213652387565, 0.000535045478866, 0, 0, 0.000219346347883,      0.000359669104424, 0.000269469145783, 0.00016051364366, 0.000181950509461,      0.000385579332948, 0.00078170803205, 0.000747244535774, 0, 0.000622858922454,      0.000689084895259, 0.000487983408564, 0.000224744353298, 0.000416449765747,      0.000308388157895, 0.000198906016907, 0.000179549331179, 9.06289650172e-05,      0.000253506844685, 0.000582896161212, 0.000386473429952, 0.000179839942451,      0, 0.000275608635737, 0.000622665006227, 0.00036075036075, 0.00029057097196,      0.000353232073472, 0.000394710874285, 0.000207555002076, 0.000402738622634,      0, 0.000309693403531, 0.000506521463847, 0.000226988991034, 0.000414164423276,      9.6590360282e-05, 0.000476689865573, 0.000377572210685, 0.000378967314069,      9.25240562546e-05, 0.000172309813044, 0.000447627573859, 0, 0.000589333071408,      0.000191699415317, 0.000362943471554, 0.000287549122975, 0.000311688311688,      0.000724112961622, 0.000434656621269, 0.00122292103424, 0.00177549812586,      0.00308008213552, 0.00164338537387, 0.00176056338028, 0.00180072028812,      0.00258939580764, 0.00217548948513, 0.00493015612161, 0.00336344416683,      0.00422716412424, 0.00313360554553, 0.00540144648906, 0.00425728829246,      0.0046828437633, 0.00397219463754, 0.00501656412683, 0.00492700729927,      0.00224424911165, 0.000634696755994, 0.00120550276557, 0.00125313283208,      0.00164551010813, 0.00143575017947, 0.00237006940918, 0.00236686390533,      0.00420336269015, 0.00329840900272, 0.00242005185825, 0.00326554846371,      0.006217237596, 0.0037103784586, 0.0038714672861, 0.00455830066551,      0.00361747518783, 0.00304147465438, 0.00476801760499, 0.00569875504121,      0.00583855136233, 0.0050566695728, 0.0042220072126, 0.00408237321963,      0.00255222610833, 0.00123507616303, 0.00178136133508, 0.00147434637311,      0.00126742712294, 0.00186590371937, 0.00177226406735, 0.00249154653853,      0.00549127279859, 0.00349072202829, 0.00348027842227, 0.00229555236729,      0.00336862367661, 0.00383477593952, 0.00273999412858, 0.00349618180145,      0.00376108175875, 0.00383351588171, 0.00368928059028, 0.00480028982882,      0.00388823582602, 0.00745054380406, 0.0103754506287, 0.00822677278011,      0.00778350981989, 0.0041831792162, 0.00537228238059, 0.00723645609231,      0.0144428396845, 0.00893333333333, 0.0106231171714, 0.0158367059652,      0.01811729548, 0.0207095263821, 0.0211700064641, 0.017604180993,      0.0165804327375, 0.0188679245283, 0.0191859923629, 0.0269251008595,      0.0351239669421, 0.0283510318573, 0.0346557651212, 0.0270022042616,      0.0260845175767, 0.0349758630112, 0.0207069247809, 0.0106362024818,      0.00981093510475, 0.00916507201128, 0.00887198986058, 0.0073929115025,      0.00659077291791, 0.00716191546131, 0.00942304513143, 0.0106886280007,      0.0123527175979, 0.0171022290546, 0.0142909490656, 0.0157642220699,      0.0265140538974, 0.0194395354708, 0.0241685144124, 0.0229897123662,      0.017921889568, 0.0155115839714, 0.0145263157895, 0.017609281127,      0.0157671315949, 0.0190258751903, 0.0138453217956, 0.00958058335108,      0.0122924304507, 0.00929741151611, 0.00885235535884, 0.00509319462505,      0.0061314863177, 0.0063104189044, 0.00729117134253, 0.010843373494,      0.0217755443886, 0.0181687353841, 0.0155402963498, 0.017310022503,      0.0214746959003, 0.026357827476, 0.0194751217195, 0.0196820590462,      0.0184317400812, 0.0130208333333, 0.0128666035951, 0.0120045731707,      0.0122374253228, 0.00874940561103, 0.0114368092263, 0.00922893718369,      0.00479041916168, 0.00644107774653, 0.00775830595108, 0.00829578041786,      0.00681348095875, 0.00573782551125, 0.00772002058672, 0.0112488083889,      0.00908907291456, 0.0157722638969, 0.00994270306707, 0.0134179772039,      0.0126050420168, 0.0113648781554, 0.0153894803415, 0.0126959699913,      0.0116655865198, 0.0112065745237, 0.0122006737686, 0.010251878038,      0.010891174691, 0.0148273273273, 0.0138516532618, 0.0136552722011,      0.00986993819758, 0.0097852677358, 0.00889011089726, 0.00816723383568,      0.00917641660931, 0.00884466556108, 0.0182179529646, 0.0183156760639,      0.0217806648835, 0.0171099125907, 0.0186579938377, 0.019360390076,      0.0144603654529, 0.0177730696798, 0.0153226598566, 0.0134016909516,      0.0126480805202, 0.0115501519757, 0.0127156322248, 0.0124326204138,      0.0240245215806, 0.0130234933606, 0.0144222706691, 0.00854005693371,      0.0053560967445, 0.00504132231405, 0.00288778877888, 0.00593526847816,      0.00455653279644, 0.00433014040152, 0.00535770564135, 0.0131095962244,      0.0126319758673, 0.0154982879798, 0.0125940464508, 0.0169948745616,      0.0257535512184, 0.0256175663312, 0.0265191262043, 0.0228974403622,      0.0193122555411, 0.0165794768612, 0.015658837248, 0.0168208578638,      0.0129912843282, 0.0119498443154, 0.0112663755459, 0.00838112042347,      0.00925767186696, 0.0113408269771, 0.0210861519924, 0.0156036134684,      0.0121687119728, 0.011006497812, 0.0107891491985, 0.0134615384615,      0.0147229755909, 0.015756893641, 0.0176257128046, 0.016776075857,      0.0169553999263, 0.0179193118984, 0.0190055672874, 0.0183088625509,      0.0155489923558, 0.0152507401094, 0.0160748342567, 0.0161532350605,      0.0139190952588, 0.0161469457497, 0.0118186629035, 0.0109259765092,      0.00950587391265, 0.00928986154533, 0.00815520645549, 0.00702576112412,      0.00709539362541, 0.00827287768869, 0.0104688211197, 0.0130375888927,      0.0160891089109, 0.0188415910677, 0.0203265044814, 0.0183175033921,      0.0139940353292, 0.0124648170487, 0.0131685758095, 0.00957428620277,      0.0119647893342, 0.00835800104475, 0.0101892285298, 0.00904207699194,      0.00772134522992, 0.00740740740741, 0.00776823249863, 0.00642254601227,      0.00484237572883, 0.00361539964823, 0.00414811817078, 0.00358072916667,      0.00433306007729, 0.00485008818342, 0.00905280804694, 0.00931847250137,      0.00779271381259, 0.00779912497622, 0.00908230842006, 0.0058152538582,      0.0102777777778, 0.00807537012113, 0.00648535564854, 0.0145492582731,      0.00694127317563, 0.00759878419453, 0.00789242911429, 0.00635050701629,      0.00785233530492, 0.00607964332759, 0.00531968282646, 0.00361944157187,      0.00305157155935, 0.00276327909119, 0.00318820364651, 0.00184464029514,      0.00412550211703, 0.00516567972786, 0.00463655399342, 0.00702897308418,      0.0100714154917, 0.00791168353266, 0.00959190791768, 0.00736,      0.00738007380074, 0.012573964497, 0.0117919562013, 0.00842919476398,      0.00778887565289, 0.00623967700496, 0.0062232955601, 0.00447815755803,      0.00511135450894, 0.00502557659517, 0.00330328263712), .Tsp = c(1,      15.9583333333333, 24), class = "ts")          T2 <- structure(c(0, 0, 0, 0, 0.000109673173942, 0, 0, 0, 0, 0, 0,      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9.66183574879e-05, 0,      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9.43930526713e-05,      0, 0, 0, 8.95255147717e-05, 0, 0, 0, 0, 0.000191699415317, 0.000207792207792,      0, 0, 0, 0.00019727756954, 0.000205338809035, 0.000205423171734,      0.000704225352113, 0.000450180072029, 0.000493218249075, 0.000120860526952,      0.000410846343468, 0.000384393619066, 0.000643264105863, 0.000189915487608,      0.000915499404925, 0.000185099490976, 0.000936568752661, 0.000451385754266,      0.000757217226692, 0.000273722627737, 0.000187020759304, 0.000211565585331,      0.000141823854772, 9.63948332369e-05, 0.000117536436295, 0.000287150035894,      0, 0, 0.000400320256205, 0.000388048117967, 0.000345721694036,      0.000296868042155, 0.000609533097647, 0.000424043252412, 0.000290360046458,      0.000546996079861, 0.000556534644282, 0.00036866359447, 0.000275077938749,      0.000964404699281, 0.00152310035539, 0.00113339145597, 0.00061570938517,      0.000362877619523, 0.000472634464505, 0.000102923013586, 0.000187511719482,      0.000294869274622, 0.00011522064754, 0.000248787162582, 0, 0.00035593521979,      0.000392233771328, 0.000551166636046, 0.000165727543918, 0.000143472022956,      0.00012030798845, 0.000438260107374, 0.000195713866327, 0.000184009568498,      0.000537297394108, 0.000365096750639, 0.000102480016397, 0.000452857531021,      0.000180848177955, 0.000770745910765, 0.00219818869252, 0.000357685773048,      0.000362023712553, 0.000660501981506, 0.000419709560984, 0.000488949735967,      0.00177758026886, 4e-04, 0.000475661962898, 0.000879816998064,      0.0014942099365, 0.00378173960022, 0.00274725274725, 0.00192545729611,      0.0016462841016, 0.00176238855484, 0.00260780478718, 0.00447289949132,      0.0034435261708, 0.00290522941294, 0.002694416055, 0.0041329904482,      0.00729244577412, 0.0296930503689, 0.00982375036117, 0.00453023439039,      0.00327031170158, 0.00221573169503, 0.00211237853823, 0.00108719286801,      0.00131815458358, 0.000983008004494, 0.00132253265002, 0.00227790432802,      0.00247054351957, 0.00307455803228, 0.0029314767314, 0.00222755311857,      0.00492610837438, 0.00454430699318, 0.00753880266075, 0.00671845475541,      0.00590490003108, 0.00288356368698, 0.00294736842105, 0.00248601615911,      0.00197089144936, 0.00326157860404, 0.00302866414278, 0.00202256759634,      0.00258788009489, 0.00169043845747, 0.00137000737696, 0.000433463372345,      0.000908368343363, 0.000805585392052, 0.00142653352354, 0.00189328743546,      0.00558347292016, 0.00161899622234, 0.00162631008312, 0.00276960360048,      0.00585673524553, 0.00519169329073, 0.0045125282033, 0.00562344544176,      0.00322815786733, 0.00330528846154, 0.00255439924314, 0.00285823170732,      0.00240894199268, 0.00218735140276, 0.00201826045171, 0.00168701002282,      0.000460617227084, 0.00127007166833, 0.00109529025192, 0.000819336337567,      0.00158170093685, 0.000588494924231, 0.00120089209127, 0.00305052430887,      0.00161583518481, 0.00211579149837, 0.0010111223458, 0.00346270379455,      0.00228091236495, 0.00207627581685, 0.00295140718878, 0.0022121765894,      0.00240718451995, 0.00224131490474, 0.0031867431485, 0.00176756517897,      0.00233382314807, 0.00178303303303, 0.00169794459339, 0.00162778079219,      0.000737939304492, 0.00135906496331, 0.000733205022454, 0.000875060768109,      0.00114705207616, 0.000967385295744, 0.00182179529646, 0.00359130903214,      0.00420328620558, 0.00446345545843, 0.00376583361862, 0.00659687365553,      0.00433810963586, 0.00353107344633, 0.00333955407131, 0.00341788091383,      0.0024939877082, 0.00538428137212, 0.00906989151698, 0.00773778473309,      0.0210421671775, 0.00859720803541, 0.00511487506289, 0.00406669377796,      0.00117164616286, 0.00206611570248, 0.00107260726073, 0.00148381711954,      0.000741761152909, 0.00104973100643, 0.00110305704381, 0.00209753539591,      0.00452488687783, 0.00486574157506, 0.00850507033039, 0.0101159967629,      0.0163991223005, 0.0150452373691, 0.0156443766097, 0.0112310639039,      0.00635593220339, 0.00627766599598, 0.00583041812427, 0.00622371740959,      0.00624897220852, 0.00420769166036, 0.00305676855895, 0.00291133656815,      0.00120006857535, 0.00501806503412, 0.00490575781048, 0.00593119810202,      0.00226874291018, 0.00304999336958, 0.00339087546239, 0.00541958041958,      0.00445563734986, 0.00431438754455, 0.0038016243304, 0.0037928519329,      0.00491460867428, 0.00460782305959, 0.00508734881935, 0.00300725278613,      0.00390896455872, 0.00367811967345, 0.00953591862683, 0.00529614264278,      0.00243584167029, 0.00427167876976, 0.00291056623743, 0.00227624510607,      0.00439422473321, 0.00232246538633, 0.00317623830372, 0.00263466042155,      0.00180200473026, 0.00190912562047, 0.0034896070399, 0.00338638672536,      0.00548090523338, 0.00697836706211, 0.00720230473752, 0.00746268656716,      0.00367056664373, 0.0032167269803, 0.00523135203391, 0.00299196443837,      0.00299119733356, 0.00287306285913, 0.00154657933042, 0.00214861235452,      0.00163006177076, 0.00157407407407, 0.00137086455858, 0.00124616564417,      0.000790591955727, 0.00107484854407, 0.00121408336706, 0.00108506944444,      0.00105398758637, 0.000881834215168, 0.00184409052808, 0.00237529691211,      0.0013637249172, 0.00190222560396, 0.00264900662252, 0.00156564526951,      0.00263888888889, 0.00183531139117, 0.00303347280335, 0.0120768352986,      0.00365330167139, 0.00351443768997, 0.00263080970476, 0.0029703984431,      0.00265143789517, 0.0014185834431, 0.00150557061126, 0.00144777662875,      0.00111890957176, 0.000716405690308, 0.000797050911627, 0.000512400081984,      0.000868526761481, 0.00113392969636, 0.00134609632067, 0.00240013715069,      0.00128181651712, 0.00110395584177, 0.00156958493198, 0.00208,      0.00184501845018, 0.00110946745562, 0.000736997262582, 0.00208250694169,      0.00229084578026, 0.00137639933933, 0.00111462010032, 0.000822518735149,      0.00200803212851, 0.000987166831194, 0.00041291032964), .Tsp = c(1,      15.9583333333333, 24), class = "ts")          T3 <- structure(c(0.00192287148809, 0.00149812734082, 0.00192410475681,      0.00151122625216, 0.00120640491336, 0.00167845582065, 0.00121261115602,      0.000802568218299, 0.00109170305677, 0.00250626566416, 0.00273597811218,      0.00242854474127, 0.00160915430002, 0.00124571784491, 0.00192943770673,      0.00329388800781, 0.00191032700303, 0.00156168662155, 0.00174753289474,      0.0014917951268, 0.00143639464943, 0.000543773790103, 0.000929525097178,      0.00141560496294, 0.000966183574879, 0.000719359769805, 0.00190740419629,      0.00137804317869, 0.00197177251972, 0.001443001443, 0.00203399680372,      0.00158954433063, 0.00256562068285, 0.00228310502283, 0.00302053966975,      0.00227352221056, 0.00263239393001, 0.00202608585539, 0.00272386789241,      0.00269206875129, 0.0027045300879, 0.00276480122033, 0.00405890126487,      0.00341070582662, 0.00351591413768, 0.00336004135436, 0.00358102059087,      0.00257289879931, 0.00235733228563, 0.00239624269146, 0.00136103801833,      0.000862647368926, 0.00145454545455, 0.00168959691045, 0.00246305418719,      0.0020964360587, 0.00335371868219, 0.00390143737166, 0.00349219391947,      0.00334507042254, 0.00255102040816, 0.00332922318126, 0.00386753686246,      0.00246507806081, 0.00432442821449, 0.00312442565705, 0.00408318298357,      0.00375354756019, 0.00416473854697, 0.00263942103023, 0.0028888688273,      0.00321817321344, 0.00310218978102, 0.002150738732, 0.00296191819464,      0.00134732662034, 0.00221708116445, 0.00152797367184, 0.00157932519742,      0.00220077873709, 0.00207100591716, 0.00260208166533, 0.00310438494373,      0.00311149524633, 0.00385928454802, 0.00292575886871, 0.00222622707516,      0.00329074719319, 0.00282614641262, 0.00287542899545, 0.00221198156682,      0.00311754997249, 0.00315623356128, 0.00287696733796, 0.00296425457716,      0.00263875450787, 0.00208654631226, 0.00179601096512, 0.00164676821737,      0.00206262891431, 0.00235895419697, 0.00241963359834, 0.0028610523697,      0.00516910352976, 0.00160170848905, 0.00254951951363, 0.00275583318023,      0.00298309579052, 0.00286944045911, 0.00288739172281, 0.00394434096636,      0.00254428026226, 0.00285214831171, 0.0034924330617, 0.00246440306681,      0.00266448042632, 0.00389457476678, 0.00253187449136, 0.00171276869059,      0.00184647850171, 0.00134132164893, 0.00153860077835, 0.000990752972259,      0.00117518677075, 0.00312927831019, 0.00188867903566, 0.0024,      0.00269541778976, 0.00263945099419, 0.00242809114681, 0.00378173960022,      0.00274725274725, 0.00165039196809, 0.00211665098777, 0.00290275761974,      0.00149017416411, 0.00105244693913, 0.00309917355372, 0.00240432779002,      0.00297314875035, 0.0015613519471, 0.00196335078534, 0.00227707441479,      0.00279302706347, 0.00295450068938, 0.00316811446091, 0.00211501661799,      0.00168990283059, 0.00195694716243, 0.00131815458358, 0.00112343771942,      0.00214911555629, 0.00157701068863, 0.00171037628278, 0.00230591852421,      0.00183217295713, 0.00102810143934, 0.00130396986381, 0.00151476899773,      0.00188470066519, 0.00220449296662, 0.00238267895991, 0.00238639753406,      0.00147368421053, 0.00113942407292, 0.0018192844148, 0.00152207001522,      0.00151433207139, 0.00117096018735, 0.000862626698296, 0.00095087163233,      0.00137000737696, 0.00119202427395, 0.00170319064381, 0.000805585392052,      0.0012680297987, 0.00189328743546, 0.00186115764005, 0.000719553876597,      0.000903505601735, 0.000865501125151, 0.00210241778045, 0.00146432374867,      0.00130625816411, 0.0011895749973, 0.00135374362178, 0.00120192307692,      0.00160832544939, 0.0015243902439, 0.00240894199268, 0.00218735140276,      0.00230658337338, 0.00188548179022, 0.0016582220175, 0.00263086274154,      0.00155166119022, 0.00204834084392, 0.00194670884536, 0.00308959835221,      0.00154400411734, 0.00152526215443, 0.00343364976772, 0.00269282554337,      0.00235928547354, 0.00230846919636, 0.00300120048019, 0.00327833023713,      0.00347844418678, 0.00259690295277, 0.00157392833997, 0.00345536047815,      0.00336884275699, 0.0023862129916, 0.00216094735932, 0.00478603603604,      0.00330652368186, 0.00551636824019, 0.00313624204409, 0.00253692126484,      0.00201631381175, 0.00243072435586, 0.00229410415233, 0.00386954118297,      0.00298111957602, 0.00305261267732, 0.0038211692778, 0.00334759159383,      0.00479287915098, 0.0045891294995, 0.00525831471014, 0.00800376647834,      0.0076613299283, 0.00638604065479, 0.00587868531219, 0.00633955709944,      0.00453494575849, 0.00617283950617, 0.00314804075884, 0.00425604358189,      0.00536642629549, 0.00422936152908, 0.00234329232572, 0.00454545454545,      0.00305280528053, 0.00389501993879, 0.0040267034015, 0.00275554389188,      0.00409706901986, 0.00506904387345, 0.0065987933635, 0.00594701748063,      0.00343473994112, 0.00579983814405, 0.00750664048966, 0.00365965233303,      0.00467423447486, 0.00348250043531, 0.00464471968709, 0.00603621730382,      0.00358154256205, 0.00445752733389, 0.00501562243052, 0.0035344609947,      0.00410480349345, 0.00467578297309, 0.00265729470255, 0.00210758731433,      0.00223771408899, 0.00218998083767, 0.00309374033206, 0.00291738496221,      0.00184956843403, 0.00297202797203, 0.00329329717164, 0.00318889514162,      0.00397442543632, 0.00481400437637, 0.002580169554, 0.00440303092361,      0.00335956997504, 0.00318415000884, 0.00269284225156, 0.00242217637032,      0.00381436745073, 0.00238326418925, 0.0037407568508, 0.00290474156343,      0.00335156112189, 0.00227624510607, 0.00376647834275, 0.00223313979455,      0.00197441840501, 0.00214676034348, 0.00225250591283, 0.00140002545501,      0.0034896070399, 0.00220115137149, 0.002828854314, 0.00418702023726,      0.00176056338028, 0.00393487109905, 0.00217939894471, 0.00331724969843,      0.00234508884279, 0.00282099504189, 0.00239295786685, 0.00269893783737,      0.00263828238719, 0.00250671441361, 0.00231640356898, 0.00231481481481,      0.00127947358801, 0.0017254601227, 0.00207530388378, 0.00185655657612,      0.00131525698098, 0.00227864583333, 0.0018737557091, 0.00220458553792,      0.00184409052808, 0.00109629088251, 0.00253263198909, 0.00228267072475,      0.00170293282876, 0.00134198165958, 0.000833333333333, 0.00269179004038,      0.00198744769874, 0.00209205020921, 0.00146132066855, 0.00113981762918,      0.00185131053298, 0.00194612311789, 0.00203956761167, 0.00111460127673,      0.00170631335943, 0.00186142709411, 0.00183094293561, 0.00194452973084,      0.0014944704593, 0.00153720024595, 0.00184561936815, 0.00151190626181,      0.000897397547113, 0.00222869878279, 0.00201428309833, 0.00202391904324,      0.00244157656087, 0.00256, 0.00184501845018, 0.00160256410256,      0.00115813855549, 0.0016858389528, 0.001741042793, 0.0026610387227,      0.00167193015047, 0.00201060135259, 0.00219058050383, 0.00233330341919,      0.000963457435827), .Tsp = c(1, 15.9583333333333, 24), class = "ts")      I know that T1 and T2 are correlated and consider them as ground truth so any distance metric should tell me that (T1, T2) are closer than (T2, T3) and (T1, T3). However, when using `dtw` in R, I am getting the following:               > dtw(T1, T2, k = TRUE)$distance; dtw(T1, T3, k = TRUE)$distance; dtw(T3, T2, k = TRUE)$distance     [1] 1.107791     [1] 1.568011     [1] 0.4102962      Can someone please explain how to use Dynamic Time Warping for nearest- neighbor queries?


title: The cross validation (CV) and the generalized cross validation (GCV) statistics

text: I have found possibly conflicting definitions for the cross validation (CV) statistic and for the generalized cross validation (GCV) statistic associated with a linear model $Y = X\boldsymbol\beta + \boldsymbol\varepsilon$ (with a normal, homoscedastic error vector $\boldsymbol\varepsilon$). On the one hand, Golub, Heath & Wahba define the GCV estimate $\hat{\lambda}$ as (p. 216) > the minimizer of $V\left(\lambda\right)$ given by $$ V\left(\lambda\right) = > \frac{\frac{1}{n} \left\|\left(I - > A\left(\lambda\right)\right)y\right\|^2}{\left(\frac{1}{n} > \mathrm{tr}\left(I - A\left(\lambda\right)\right)\right)^2} $$ where > $A\left(\lambda\right) = X\left(X^T X + n\lambda I\right)^{-1} X^T$ On the other hand, Efron defines the same concept as $V\left(0\right)$ (p. 24), yet he attributes the introduction of this concept to Craven & Wahba, where its definition (p. 377) is essentially the same as Golub, Heath & Wahba's above mentioned definition. Does this mean that $0$ minimizes $V\left(\lambda\right)$? Similarly, Golub, Heath & Wahba define the CV estimate of $\lambda$ (p. 217) as the minimizer of $$ P\left(\lambda\right) = \frac{1}{n}\sum_{k=1}^n \left(\left[X \beta^{(k)}\left(\lambda\right)\right]_k - y_k\right)^2 $$ where $\beta^{\left(k\right)}\left(\lambda\right)$ is the estimate $$ \hat{\beta}\left(\lambda\right) = \left(X^T X + n \lambda I\right)^{-1} X^T y $$ of $\beta$ with the $k$th data point $y_i$ omitted. The authors attribute the introduction of the CV estimate (also called the PRESS estimate) to Allen ("Allen's PRESS", ibid.) Yet in Allen's paper, the PRESS estimate is defined (p. 126) as $n P\left(0\right)$ (in Efron's article it is defined as $P\left(0\right)$ (p. 24)). Again, does this mean that $0$ minimizes $P\left(\lambda\right)$? * * *   1. Allen, David M. The Relationship Between Variable Selection and Data Agumentation and a Method for Prediction. Technometrics, Vol. 16, No. 1 (February, 1974), pp. 125-127   2. Craven, Peter and Wahba, Grace. Smoothing Noisy Data with Spline Functions. Numerische Mathematik 31, (1979), pp. 377-403   3. Efron, Bradley. How Biased Is the Apparent Error Rate of a Logistic Regression? Technical report no. 232. Department of Statistics, Stanford University (April 1985)   4. Golub, Gene H., Heath and Grace Wahba. Generalized Cross-Validation as a Method for Choosing a Good Ridge Parameter. Technometrics, Vol. 21, No. 2 (May, 1979), pp. 215-223


title: What is shallow architecture in machine learning?

text: What is a precise definition of shallow architecture in machine learning?


title: Deviance vs Pearson goodness-of-fit

text: I am trying to come up with a model by using negative binomial regression (negative binomial GLM). I have a relatively small sample size (greater than 300), and the data are not scaled. I noticed that there are two ways to measure goodness of fit - one is deviance and the other is the Pearson statistic. How can I determine which goodness-of-fit measure to use? Are there some criteria that I can take a look at in selecting the goodness-of-fit measure?


title: Power of $\chi^2$ tests for homogeneity and independence

text: I recently studied these topics and I was wondering what their power is. My opinion is that since we are relying on an approximate chi-square distribution it cannot be much and so we have to be extra cautious. Is that right? Thanks.


title: How to choose the most appropriate distribution for a given vector in R?

text: When creating a Monte Carlo simulation model for a variable, a critical step is to choose the distribution that best fits the variable's probability density. I generally do this by looking at the density plot and determining what distribution best fits the density shape. For (a very lame) example, this…               x <- rnorm(1000)     plot(density(x))      … appears to be a normal distribution (but only because it's a random sample from the normal distribution). However, when dealing with real world data, it's difficult to know which of the 17 built-in distributions best represent the shape of the data. For example, this data…               data <- c(6.515, 0.243, 0.725, 2.276, 1.456, 4.047, 0.766, 0.29, 2.368,      0.543, 2.223, 0.488, 0.47, 3.511, 0.544, 4.191, 0.414, 0.704,      4.917, 0.434, 0.773, 0.477, 3.257, 0.415, 1.921, 0.278, 3.159,      4.193, 0.132, 1.109, 1.538, 4.088, 0.468, 0.047, 2.204, 3.765,      0.168, 2.231, 0.164, 0.371, 2.33, 4.458, 0.046, 1.195, 1.714,      1.046, 1.915, 2.66, 5.409, 0.466)          plot(density(data))      ![Plot](http://i.stack.imgur.com/l9SWm.png) … seems like it could be best modeled with the chi-squared distribution, but it could also be a gamma distribution. The only way I've found to fit the best type of model is to overlay a bunch of different possible distributions until I see one that visually matches (or comes close). But surely there's a more numerical, official way to do that, right? Is there a systematic, non-visual (and automated) way to find the best distribution for a given variable? Is there a function in some R function that runs through different distributions to check their goodness of fit, or is that terribly inefficient?


title: What is the relationship between profile likelihood and confidence intervals?

text: To make this chart I generated random samples of different size from a normal distribution with mean=0 and sd=1. Confidence intervals were then calculated using alpha cutoffs ranging from .001 to .999 (red line) with the t.test() function, the profile likelihood was calculated using the code below which I found in lecture notes put on line (I can't find the link at the moment **Edit:**Found it), this is shown by the blue lines. Green lines show the normalized density using the R density() function and the data is shown by the boxplots at the bottom of each chart. On the right is a caterpillar plot of the 95% confidence intervals (red) and 1/20th of max likelihood intervals (blue). **R Code used for profile likelihood:**                 #mn=mean(dat)       muVals <- seq(low,high, length = 1000)       likVals <- sapply(muVals,                         function(mu){                           (sum((dat - mu)^2) /                              sum((dat - mn)^2)) ^ (-n/2)                         }       )      ![enter image description here](http://i.stack.imgur.com/silT5.png) My specific question is whether there is a known relationship between these two types of intervals and why the confidence interval appears to be more conservative for all cases except when n=3. Comments/answers about whether my calculations are valid (and a better way to do this) and the general relationship between these two types of intervals are also desired. **R code:**               samp.size=c(3,4,5,10,20,1000)     cnt2<-1     ints=matrix(nrow=length(samp.size),ncol=4)     layout(matrix(c(1,2,7,3,4,7,5,6,7),nrow=3,ncol=3, byrow=T))     par(mar=c(5.1,4.1,4.1,4.1))     for(j in samp.size){                 #set.seed(200)       dat<-rnorm(j,0,1)       vals<-seq(.001,.999, by=.001)       cis<-matrix(nrow=length(vals),ncol=3)       cnt<-1       for(ci in vals){         x<-t.test(dat,conf.level=ci)$conf.int[1:2]         cis[cnt,]<-cbind(ci,x[1],x[2])         cnt<-cnt+1       }                 mn=mean(dat)       n=length(dat)       high<-max(c(dat,cis[970,3]), na.rm=T)       low<-min(c(dat,cis[970,2]), na.rm=T)       #high<-max(abs(c(dat,cis[970,2],cis[970,3])), na.rm=T)       #low<--high                 muVals <- seq(low,high, length = 1000)       likVals <- sapply(muVals,                         function(mu){                           (sum((dat - mu)^2) /                              sum((dat - mn)^2)) ^ (-n/2)                         }       )                 plot(muVals, likVals, type = "l", lwd=3, col="Blue", xlim=c(low,high),            ylim=c(-.1,1), ylab="Likelihood/Alpha", xlab="Values",            main=c(paste("n=",n),                    "True Mean=0 True sd=1",                    paste("Sample Mean=", round(mn,2), "Sample sd=", round(sd(dat),2)))       )       axis(side=4,at=seq(0,1,length=6),            labels=round(seq(0,max(density(dat)$y),length=6),2))       mtext(4, text="Density", line=2.2,cex=.8)            lines(density(dat)$x,density(dat)$y/max(density(dat)$y), lwd=2, col="Green")       lines(range(muVals[likVals>1/20]), c(1/20,1/20), col="Blue", lwd=4)       lines(cis[,2],1-cis[,1], lwd=3, col="Red")       lines(cis[,3],1-cis[,1], lwd=3, col="Red")       lines(cis[which(round(cis[,1],3)==.95),2:3],rep(.05,2),              lty=3, lwd=4, col="Red")       abline(v=mn, lty=2, lwd=2)       #abline(h=.05, lty=3, lwd=4, col="Red")       abline(h=0, lty=1, lwd=3)       abline(v=0, lty=3, lwd=1)            boxplot(dat,at=-.1,add=T, horizontal=T, boxwex=.1, col="Green")       stripchart(dat,at=-.1,add=T, pch=16, cex=1.1)            legend("topleft", legend=c("Likelihood"," Confidence Interval", "Sample Density"),              col=c("Blue","Red", "Green"), lwd=3,bty="n")            ints[cnt2,]<-cbind(range(muVals[likVals>1/20])[1],range(muVals[likVals>1/20])[2],                          cis[which(round(cis[,1],3)==.95),2],cis[which(round(cis[,1],3)==.95),3])       cnt2<-cnt2+1     }     par(mar=c(5.1,4.1,4.1,2.1))               plot(0,0, type="n", ylim=c(1,nrow(ints)+.5), xlim=c(min(ints),max(ints)),           yaxt="n", ylab="Sample Size", xlab="Values")     for(i in 1:nrow(ints)){       segments(ints[i,1],i+.2,ints[i,2],i+.2, lwd=3, col="Blue")       segments(ints[i,3],i+.3,ints[i,4],i+.3, lwd=3, col="Red")     }     axis(side=2, at=seq(1.25,nrow(ints)+.25,by=1), samp.size)


title: Correlation of distribution standard deviation with subtraction of its mean from mean of normal distribution

text: The following question might have relations to this question: Given $\mu$ and $\sigma$ (mean and standard deviation) of normal distribution, find the set of all distributions with their $\mu_1$ and $\sigma_1$ that $\mu-\mu_1$ is any increasing function of $\sigma_1$. Thank you. please at least give me a clue or examples of such distributions. My question might be wrong or I could put it in a better way.


I am familiar with the Bayesian score, which is used to compare competing structures of a BN. However, I have difficulty in understanding how the Bayesian score formula below is derived. **1)** What is an intuitive explanation of the Bayesian score? (e.g. how do the fractions within the product operators equate to $P(D|B_S)$) **2)** What type of structures does it favour? (e.g. more/less connected, higher/lower number of parents, etc.) **Background:** The Bayesian score of a BN structure, $B_S$, for a data set, $D$, is used to calculate the joint probability of $P(B_S,D)$. For a discrete network (after some assumptions*), it is given as: $$P(B_S,D) = P(B_S)\prod_{i=0}^n\prod_{j=1}^{qi} \frac{\Gamma(N'_{ij})}{\Gamma(N'_{ij} + N_{ij})} \prod_{k=1}^{ri} \frac{\Gamma(N'_{ijk} + N_{ijk})}{\Gamma(N'_{ijk})} $$   * $\Gamma()$ is the gamma function. Above, since all counts are positive integers, $\Gamma(int) = (int-1)!$   * $P(B_S)$ is the prior on the network structure ( taken to be uniform/uninformative by many versions of the Bayesian score)   * $xi$ is a node in the network with $n$ nodes;    * $ri (1 ≤ i ≤ n)$ is the cardinality of $xi$;    * $qi$ denotes the cardinality of the parent set of $xi$ , that is, the number of different values to which the parents of xi, $pa(xi)$, can be instantiated.   * $N_{ij} (1 ≤ i ≤ n, 1 ≤ j ≤ qi)$ denotes the number of records in $D$ for which $pa(xi)$ takes its $j$th value.    * $N_{ijk} (1 ≤ i ≤ n, 1 ≤ j ≤ qi, 1 ≤ k ≤ ri)$ denotes the number of records in $D$ for which $pa(xi)$ takes its $j$th value and for which $xi$ takes its $k$th value. So, $N_{ij} = \sum_{k=1}^{ri} N_{ijk}$.   * $N′_{ij}$ and $N′_{ijk}$ represent choices of priors on counts restricted by $N'_{ij} = \sum_{k=1}^{ri} N'_{ijk}$. *List of asssumptions: 1 - Multinomial Sample; 2- Parameter Independence; 3- Parameter Modularity; 4 - Dirichlet; 5- Complete Data


title: Calculating Regression Coefficients for Very Large Observation Matrices

text: I am trying to run a regression from a 11300x21500 observation matrix (where there are 11300 observations and 21500 independent variables). However, when I try to implement the usual $(X^T X)^{-1}$ formula in C++, I can't even initialize the observation matrix $X$ due to memory limitations (The application closes itself when the memory usage is around 1 GB. I have 8 GB memory but I think there is an OS command to limit memory usage of each application). Is there any other technique which is more suitable for such a regression? Can there be any way to overcome this situation? Thanks in advance


title: Minimum number of observations to determine long range dependency

text: What is the minimum number of observations to be considered in order to determine long range dependency? I'm trying to estimate Hurst parameter using R/S method. I've used SELFIS as well as R. Logically, since we are interested in _long range_ dependency, we should have _enough_ observations. SELFIS requires at least 64 values before it could analyze. So, is that a standard value (or is there any standard value)? Also, for the `rsFit` method of `fArma` package in R, does the `levels` parameter indicate that minimum value?


title: Clustering a dataset with both discrete and continuous variables

text: I have a dataset X which has 10 dimensions, 4 of which are discrete values. In fact, those 4 discrete variables are ordinal, i.e. a higher value implies a higher/better semantic. 2 of these discrete variables are categorical in the sense that for each of these variables, the distance e.g. from 11 to 12 is not the same as the distance from 5 to 6. While a higher variable value implies a higher in reality, the scale is not necessarily linear (in fact, it is not really defined). My question is:   * Is it a good idea to apply a common clustering algorithm (e.g. K-Means and then Gaussian Mixture (GMM)) to this dataset which contains both discrete and continuous variables? If not:   * Should I remove the discrete variables and focus only on the continuous ones?   * Should I better discretize the continuous ones and use a clustering algorithm for discrete data?


title: Parameter learning of Markov random field

text: Given a Markov random field $\mathcal{G} = (\mathcal{V},\mathcal{E})$, the corresponding density function to which is expressed by $f(x) \propto \prod_{x\in\mathcal{V}} \psi_u(x) \prod_{(x_i,x_j)\in\mathcal{E}} \psi_p(x_i,x_j; \theta)$ where $\psi_u$ is the unary potential and $\psi_p$ is the pairwise potential. Now assume $\psi_u$ is known, the structure of the MRF and the functional form of $\psi_p$ (which is Gaussian) are also given, is there any efficient method to estimate the common hyper-parameter $\theta$ for all the pairwise potentials? MCMC is a good option, however, it's too costly in my case. I prefer some methods of computing the point-estimate of $\theta$.


title: Is there a method for constructing decision trees that takes account of structured/hierarchical/multilevel predictors?

text: Is there a method for constructing decision trees that takes account of structured/hierarchical/multilevel predictors, that would allow me to impose domain knowledge or constraints on interactions for example?


title: What is empirical entropy?

text: In the definition of jointly typical sets (in "Elements of Information Theory", ch. 7.6, p. 195), we use $$-\frac{1}{n} \log{p(x^n)}$$ as the **empirical entropy** of an $n$-sequence with $p(x^n) = \prod_{i=1}^{n}{p(x_i)}$. I never came across this terminology before. It is not defined explicitly anywhere according to the index of the book. My question basically is: Why is empirical entropy not $-\sum_{x}{\hat p (x) \log(\hat p(x))}$ where $\hat p(x)$ is the empirical distribution? What are the most interesting differences and similarities between these two formulas? (in terms of properties they share/do not share).


title: How to compare dbscan clusters / choose epsilon parameter

text: I am currently trying to make a DBSCAN clustering using scikit learn in python. I would like to compare the different outputs when varying the epsilon parameter in order to choose the right epsilon parameter. I took as an example the iris dataset. In order to compare clusters I thought about trying to cluster with epsilon within a range (ex : 0.1, 0.2, ..., 1). Now, when I run a kmeans or a hierarchical clustering I can choose my k value by checking the gap statistic for example, or by looking at inertia and choosing a k for which there is an 'elbow' on the inertia vs k plot. My problem is that I assume this will not work anymore because the total number of points within all the clusters is not constant in DBSCAN. Indeed, depending on epsilon, the number of 'noisy sample' unclassified points will vary. As a result, I may have only a few points for a low epsilon, resulting on a very small Inertia which would be biased. I could consider gap statistic because I may be able to generate random samples with the right size each time. But I wonder if I will then leave the validity framework of the paper and I'm not sure the different clustering can still be compared. Has anyone an idea on how to compare different total sizes clustering, and more precisely the results of dbscan for different epsilon ? Would silhouette coefficient work or is it total size sensitive too ? Also, can someone create the 'dbscan' tag please ? I do not have the reputation. Thanks


title: Inferring dimension weight in a mapping from a triangle to a distribution over its vertices

text: I have a dataset $(y_i, \mathbf{X}_i)$, where $\mathbf{X}_i$ is a $3 \times n$ matrix of reals and $y_i$ takes a value in $\\{1, 2, 3\\}$. Essentially, $y_i$ represents a "selection" of the row vector $\mathbf{x}_{iy_i}$. (One way to think about this is as selecting a vertex of a triangle in $n$ dimensions.) I would like to find a general method for learning $f: \mathbb{N}^{m \times n} \rightarrow \\{1, \dots, m\\}$ (where in my case $m = 3$) that also allows me to calculate feature importance (dimensions weight) for each of the $n$ features represented in the columns of $\mathbf{X}_i$. This last constraint precludes me from, for instance, calculating a Minkowski $p$-distance between each vertex and using those distances as features. My initial stab at this problem was to map $\mathbf{X}_i$ to a matrix of naturals with some precision $p$, $\mathbf{X}^{(p)}_i$, s.t. $x^{(p)}_{ijk} = \lfloor(10^p)x_{ijk}\rfloor$, and then to a unique $n$-dimensional vector $\mathbf{z}_i$, e.g. by applying the Cantor pairing function down the columns of $\mathbf{X}^{(p)}_i$ such that $z_{ij} = \langle\langle x^{(p)}_{ij1}, x^{(p)}_{ij2}\rangle, x^{(p)}_{ij3}\rangle$. I then used an off-the-shelf supervised method, e.g. a random forest, to learn $g : \mathbb{N}^n \rightarrow \\{1, 2, 3\\}$. The essential idea here is to map each $\mathbf{X}^{(p)}_i \in \mathbb{N}^{3 \times n}$ to a unique point in $n$ dimensions, then map from those points to (probabilities distributions over) categories. EDIT: Another possible solution is to assume that there exists a fourth ($[m+1]^{th}$ for the general case) $n$ dimenional point $\mathbf{z}$ and that the parameter(s) $\phi$ of the distribution over the triangle's vertices are drawn from a distribution whose parameters are some function $d$ of the distance between each of the vertices $\mathbf{x}_{ij}$ and this fourth ($[m+1]^{th}$) point. The feature importances $\mathbf{w}$ can be incorporated into the distance parameter as weights applied to each dimension: for instance, if $d$ is euclidean distance $d(\mathbf{x}_{ij}, \mathbf{z}) = \sqrt{\sum_{k=1}^n(x_{ijk} - z_k)^2}$ then the $\mathbf{w}$-weighted distance is $d_{\mathbf{w}}(\mathbf{x}_{ij}, \mathbf{z}) = \sqrt{\sum_{k=1}^n w_k(x_{ijk} - z_k)^2}$ For concreteness, let $\phi_{k(i)} \sim \mathrm{Dirichlet}\left(d_{\mathbf{w}}(\mathbf{x}_{i1}, \mathbf{z}), \ldots, d_{\mathbf{w}}(\mathbf{x}_{im}, \mathbf{z})\right)$. Since each unique triangle triangle has some number of vertex selections associated with it, we draw a single $\phi_{k(i)}$ for all instances of triangle $k(i)$ (of which $\mathbf{X}_i$ is one). Let the vertex selections $\mathbf{y}_i \sim \mathrm{Categorical}(\phi_{k(i)})$. We place an exponential prior on each $w_l \sim \mathrm{Exponential}(\alpha)$ and a multivariate gaussian prior on $\mathbf{z}_{k(i)} \sim \mathcal{N}(\frac{1}{m}\sum_{j=1}^m \mathbf{x}_{ij}, \mathrm{Cov}(\mathbf{X}_i))$. We are interested in the posterior over the weights $p(\mathbf{w}\;|\;\mathbf{y}, \mathbf{X}; \alpha) \propto p(\mathbf{w};\alpha) \int \mathrm{d}\mathbf{z}\; p(\mathbf{z}\;|\;\mathbf{X}) \int \mathrm{d}\phi\; p(\mathbf{y}\;|\;\phi)p(\phi\;|\;\mathbf{z})$. Due to the Dirichlet- Multinomial conjugacy, the deepest integral can be analytically solved, but I suspect the integral over $\mathbf{z}$ will require brute force. Further clarity---and consequently, less bookkeeping---can be achieved by specifying this model in terms of a multinomial likelihood instead of a categorical likelihood (since we can remove the need for the indexing function $k()$).


title: Splitting an imbalanced dataset for training and testing

text: I have a highly imbalanced dataset. My question is how to split the dataset for training and testing? I want to have a separate training set and a separate testing set. One idea I had is just to shuffle the dataset and then take 70% for training and the rest for testing. Is that good strategy?


title: Modelling atemporal causal link between two (or three) indices

text: I want to construct two (or three) indices with two sets of variables (say around fifteen variables in each set) and then examine the causal link between these two indices. For example, to test the hypothesis that first indices affect negatively the second indices, how I can build such a model and test this hypothesis? Is Structural Equation Modelling a way of doing this. Please note that Granger causality may not be suitable because I am looking at atemporal causal links. Please suggest something.


title: Variational distance between samples of the same distribution

text: Suppose that I have $n$ independent samples from a probability distribution. Then, I divide them into two separate sets both with cardinality $n/2$. I compute the pdf from the first and the second group from the samples. As $n$ goes to infinity the variational distance between both pdfs should go to zero, and also both of them converge to the true pdf. Do you know if there is any result for finite $n$ that gives bounds on the variational distance between the pdfs of the two sets? (note: this is related to the previous question: Is it possible to test for a Gaussian distribution?)


title: Implications of using survival analysis with only observed data

text: I am doing survival analysis on a dataset where the number of right-censored individuals is much greater than the number of observed individuals. There are about 50,000 individuals and only about 5,000 of them are observed. So as you might imagine, the survival curves are not that informative because they never go below about 80%. So I was wondering what the implications would be by looking at only the observed data (i.e. the 5,000). I would imagine that this would artificially lower survival probabilities on the curve. However, if you are only interested in the relative difference in survival curves for different groups (I am using a regression-based SA model, where the different groups are defined by different values of the covariates) would it matter?


title: Goodness of fit to Poisson Distribution

text: What are some of the well known statistical tests to measure goodness of fit of observed random variables to a poisson distribution? I know the Kolmogorov- Smirnov test is one such, are there any others out there? Thanks!


title: All factors have low proportion variance explained. What could this mean about the data?

text: I'm running factor analysis on a survey data set. Across any number of factors, the % variance explained for all factors is very low. No factor ever has a % variance explained greater than 9%. In basic terms, this means that it's not possible to construct a factor, or principal component for that matter, that accounts for or "explains" very much of the variability in the data. Beyond that, what implications might we draw about the data? Does this suggest very high variances and/or covariances in the data? Anything else? (I'm not expecting that there is a single thing we can infer about the data given uniformly weak factors, but I'm looking for ideas so I can then explore the data.)


title: Measuring intra-observer correlation

text: We're far from what you'd call expert statisticians, so we'd like to ask some help in an ongoing medical research. We have 60 subjects. A researcher took at least 3 times the measurements of 8 of their values of interest, giving more than 60×8×3=1440 data points. As a first step we'd like to determine the intra-observer correlation (Edit: mistake, this should have been variability) of these measurements, but we hardly even know how to get started. Asking Google and Wikipedia gave us a few options: Pearson: this expects pairs of data points, but we have 3 or more measurements for each subject and each attribute. We could arrange them in pairs but we'd rather find a method that requires less manual work Intra-Class Correlation: The concept of a "group" is not so evident here. How would we group the measurements in this case? If the above isn't clear, we'd gladly elaborate more. Thank you all in advance!


title: K-means & BIC (to validate clusters) in R

text: I'm wondering if there is a good way to calculate the clustering criterion based on BIC formula, for a k-means output in R? I'm a bit confused as to how to calculate that BIC so that I can compare it with other clustering models. Currently I'm using the stats package implementation of k-means.


title: What is the difference between estimate and contrasts in SAS?

text: I am having difficulty determining the difference between the two statements.


title: Can the proof of the central limit theorem be expressed as the limit of a convolution with an operator

text: The central limit theorem, as I understand it (engineer, non-statistician), says that the distribution comprised of means of some other reasonably behaved distributions converges to a normal distribution in the limit as sample size goes to infinity. Convolution (1d) between scalar values is shown here. The mean can be considered an expectation operator. Can the central limit theorem be expressed as N(mu,sig) = conv(scalar values, operator) I am thinking to explore frequency domain deconvolution to learn about the nature of the distribution upon which the mean is being operated, but need help to determine if it can be expressed in these terms, and then to express it in a way that is useful. If so please give complete details.


title: Multiple Linear Regression with 2 predictors

text: I have a set of data and supposed to see if the linear model $$Y=\beta_0 + \beta_1 X_1 + \beta_2X_2 + \varepsilon$$ fits the data (mydata, using R) below is a scatter plot.               > mydata=read.table(....)     > pairs(mydata)      ![scatterplot using pairs\(\) function in R](http://i.stack.imgur.com/GLKMc.png) I am not very good at interpreting scatter plots(yet). From the output, how can i tell whether the model fits? thanks


title: Variogram modeling and interprtation

text: I have 23 stations and I have calculated the variogram as following: ![Variogram](http://i.stack.imgur.com/zs7RZ.png) The maximum distance between the stations is 2169 meter and the range in the variogram is 1650 meter (exponential effective range). I would like to group the stations into clusters in which in the presence of rain all the stations can observe it simultaneously. So, I have divided stations to 2 clusters where the maximum distance in each cluster is decreased to 200 meter. Will the rain be observed in all the stations simultaneously since the maximum distance is 200 m and range is 1650?


title: Re-scaling a confusion matrix after down sampling one class

text: Let's say I have a large, un-balanced binary classification problem (in reality nrow is more like 500k, and ncol is more like 500):               set.seed(42)     nrow <- 10000     ncol <- 50     X <- matrix(rnorm(nrow*ncol), ncol=ncol)     Y <- X %*% rnorm(ncol(X)) * sample(0:1, ncol(X), replace=TRUE) + rnorm(nrow(X))     Y <- Y-20     Y <- exp(Y)/(1+exp(Y))     Y <- round(Y, 0)     > sum(Y==1)/length(Y)     [1] 0.0027      Before modeling, I down-sampled the negative class. I don't have a strong theoretical justification for doing this, but it makes my models fit a lot faster, and they seem to be better too.               keep <- which(Y==1)     sample <- sample(which(Y==0), length(keep))     Xfull <- X     Yfull <- Y     X <- X[c(keep, sample),]     Y <- Y[c(keep, sample),]     > sum(Y==1)/length(Y)     [1] 0.5      Fitting a model to the down-sampled dataset is pretty quick:               library(caret)     Y <- factor(paste('X', Y, sep=''))     X <- as.data.frame(X)     model <- train(X, Y, method='glmnet',                     tuneGrid=expand.grid(.alpha=0:1, .lambda=0:30/10),                    trControl=trainControl(                      method='cv',                       summaryFunction=twoClassSummary,                       classProbs=TRUE))     plot(model)      And I can use the cross-validation folds to estimate some statistics about the model's predictive ability:               > max(model$results$ROC)     [1] 0.9777778     > confusionMatrix(model)          Cross-Validated (10 fold) Confusion Matrix           (entries are percentages of table totals)                    Reference     Prediction   X0   X1             X0 41.0  3.3             X1  9.0 46.7      However, I would like to estimate these statistics on the FULL dataset, preferably without cross-validating my model on the full dataset, which would be extremely slow. I was thinking of doing a naive re-scaling of the confusion matrix, like this:               scaling_factor <- 0.5/0.0027     CM <- confusionMatrix(model)$table * nrow(X)     CM[,1] <- CM[,1]*scaling_factor     > round(CM/sum(CM)*100, 2)               Reference     Prediction    X0    X1             X0 81.56  0.04             X1 17.90  0.50      Does this seem like a reasonable calculation? Is there a similar method I could use to re-scale AUC? Or do I expect AUC to stay the same? /edit: in response to B_Miner. I am fairly certain that fitting the downsampled model to the full dataset will overestimate its performance. It's easy to see why if we fit a random forest instead of a glmnet:               model <- train(X, Y, method='rf',                         trControl=trainControl(                          method='cv',                           summaryFunction=twoClassSummary,                           classProbs=TRUE))      And predict this model on the full dataset:               pred_full <- predict(model, Xfull, type='raw')     > table(pred_full, Yfull)              Yfull     pred_full    0    1            X0 8531    0            X1 1442   27      Because _every single positive instance_ was used to train the model, the model can perfectly predict these instances, even on the full dataset. /edit2: To clarify. I understand the the down-sampled model is biased. However, I suspect that the model's bias is predictable and consistent. I'm looking for a theoretical way to correct for this bias, under the assumption that the removed negative observations come from the same distribution as the negative observations in the training set.


title: Robust mean estimation with O(1) update efficiency

text: I am looking for a robust estimation of the mean that has a specific property. I have a set of elements for which I want to calculate this statistic. Then, I add new elements one at a time, and for each additional element I would like to recalculate the statistic (also known as an online algorithm). I would like this update calculation to be fast, preferably O(1), i.e. not dependent on the size of the list. The usual mean has this property that it can be updated efficiently, but is not robust to outliers. Typical robust estimators of the mean, like inter- quartile mean and trimmed mean cannot be updated efficiently (since they require maintaining a sorted list). I would appreciate any suggestions for robust statistics that can be calculated/updated efficiently.


title: How is the signed rank statistic defined in Matlab?

text: 1. Given a sequence of sample points $X_1, \dots, X_n$, I am trying to find a Matlab function to      * compute the signed rank statistic defined as $$ \sum_{i=1}^n sign(X_i) r_i $$ where $r_i$ is the rank of $|X_i|$ in $|X_1|, \dots, |X_n|$, and $sign(X_i)$ returns $1$ if $X_i >0$ and $-1$ if $X_i <0$.     * and also to compute the p-value of this signed rank statistic.   2. The function I suspect to be relevant is `signrank` in Matlab. Matlab signrank function takes the minimum between the sum of the ranks of positive differences and the sum of the ranks of negative differences. In the implementation of signrank function              % Compute signed rank statistic (most extreme version)     w = sum(tierank(neg));     w = min(w, n*(n+1)/2-w);      **Questions:**   1. I was wondering why Matlab signrank function takes the minimum between the sum of the ranks of positive differences and the sum of the ranks of negative differences? I checked the following references listed in the source file of signrank.m, but haven't found them taking the minimum as the signed rank statistic.              %   References:     %      [1] Hollander, M. and D. A. Wolfe.  Nonparametric Statistical     %          Methods. Wiley, 1973.     %      [2] Gibbons, J.D.  Nonparametric Statistical Inference,     %          2nd ed.  M. Dekker, 1985.        2. Is the p value provided by Matlab signrank function, same as the p value I want? Thanks and regards!


title: Algorithim to determine if point is "too far from the average"

text: Long story short I have a collection of about 30 scripts that manipulate data sets and place them in a database. These scripts report their running times as well as any errors that occur to a separate database. I wrote another script that goes through this database daily and for each script determines if an error occurred. It also checks the running time of each script 30 days back and averages them. I take the running time of the current script and see if it greater than 3 standard deviations than the average. If it is I report the running time is too far from average. Is this the correct method for performing such a task? I feel as if I get far too many "running time too far from average" errors. Would increasing the sample size help, or does the 3 standard deviations rule not apply? I was under the assumption 99% of data lies within 3 standard deviations and a reliable way to detect outliers (a script that took a "hella long time" to run) would be to use this method.


title: Conditional Forest - Calculate Overfitting, OOB Error

text: My name is Abhi and I am trying to teach myself conditional forests. I have got a basic example working myModel <\- cforest(Survived ~ .,data = modelData[,c(-1)],controls=cforest_unbiased(ntree=2000, mtry=3)) print(myModel) But I have a few questions.   * How to calculate the OOB Error?    * Is there something similar to `importance(model)` function that exists for random forest? I would like to understand the contribution of individual predictors?    * What are the tradeoffs between random forest and conditional forests. I am using R-3.1 and the party package for conditional forests. Any help would be appreciated.


title: Power/sample size calculation for one group prospective study

text: Study prospectively enrolls patients who submit for stool testing at a hospital, stools of those who have nosocomial diarrhea are processed to detect 6 pathogens to determine cause. If prevalence of ND is 12%, and initial sample size with non response rate of 10% brings size to 180. How do I calculate the power? I don't have a mean and expected value. Not sure how to calculate effect size either. I think i'm missing something basic here.


title: What statistical test to use for A/B testing

text: I'm trying to understand how A/B testing works for a simple 2x2 contingency table. I have found two different types of solutions: 1) Calculate an aggregate probability of conversion assuming that there is no difference in the conversion rate between version A and version B. Then, perform a chi-squared goodness-of-fit test between the actual results in the 2x2 table vs the predicted results (assuming no difference). (see: Chi-square independence test for A/B split testing) 2) `((Prob_conversion(A) - Prob_conversion(B))/sqrt(Standard_deviation(Prob_conversion(A))^2 + Standard_deviation(Prob_conversion(B))^2)` should have a N(0,1) distribution. I feel like method 1) and 2) must be linked somehow. Can someone explain how? Thanks!


title: How to fit an ARMAX model with more than one exogenous time series?

text: I am trying to fit an ARMAX with two exogenous time series with the following code but it gives me an "computationally singular" error. I know it is about defining more than 2 time series for `xreg` because when I include only one exogenous it works!               library(forecast)          data<-read.csv("DATA.csv")          Y1=ts(data[,1], start=1978, frequency=12)          #exogenous time series     Y2<-matrix(0,360,2)     Y2[,1]<-cbind(data[,2])     Y2[,2]<-cbind(data[,3])          model<- arima (Y1, order=c(1, 0, 0), xreg=as.ts(Y3, start=1978, frequency=12))          predict (model, 10, newxreg=0)


title: What are good basic statistics to use for ordinal data?

text: I have some ordinal data gained from survey questions. In my case they are Likert style responses (Strongly Disagree-Disagree-Neutral-Agree-Strongly Agree). In my data they are coded as 1-5. I don't think means would mean much here, so what basic summary statistics are considered usefull?


title: Combining several Likert items from a Likert scale into one variable

text: Hi there I have a very important question for my statistical analysis which needs to be finished in a few days. First of all I have to say that I have not much knowledge of statistics and therefore have problems to understand most of the answers for similar questions. My question is: I have conducted a survey with nine questions grouped into 3 variables (3 questions for each variable). Now I would like to have the score my variable which results from the 3 questions in each variable. So now I'm wondering how can I get my scores for the single variables? Can I calculate the mean of the scores from every respondent and then use the mean? Would that be okay? I would be happy if I got a quick (and easy :D) answer. Thank you!!


title: Testing correlation between two variables when there are other variables involved

text: I want to see if there is any correlation between number of pages in a paper and the citation it receives. If I just plot number of pages vs number of citations of a paper, I might get a correlation, but that might be wrong, because there are other variables such as age of the paper, venue impact factor etc. How to get around this problem? If I do a multiple regression with some known variables and report the coefficient for my variable, is that sufficient?


title: Cross validation accuracy is the same as the fraction of negative labels - what does it mean?

text: I have a dataset for classification (binary - 1/0) that has around 4000 samples that I use to train the model (I'm using an SVM, if that's relevant). To check whether everything is working fine, I cross validate 5 folds. But I'm noticing that the cross validation accuracy that I get is the exact same fraction of the 0 labels in the set. For instance, when I pick the last 800 samples for testing and train on the first 3200 samples, the accuracy of prediction on the last 800 samples is the exact same number as the fraction of samples having label = 0 in those 800 data points. Does this mean that the classifier is not taking the features into consideration at all? Or is this completely normal behavior and I'm missing something?


title: Transformation to fit gamma distribution for glm

text: The data simulated below has a maximum value of 4 and is interestingly skewed. The maximum of 4 is a limitation imposed by the instrument used and the data is semi-discrete, i.e., there are a reasonably large number of numbers it could be between -4 and 4. Because of the shape of the data, I thought about transforming it so it would approximate a gamma distribution: _Edit to update for comments:_   It is limited to this range in this instance because it is a signal detection measure (d prime http://en.wikipedia.org/wiki/D%27) and the accuracy we have for this particular measure limits us to +-4. It is skewed like this because one population does not very often get false positives and will generally get more hits while the other populations often do get false positives and less hits.               set.seed(69)     g1<-rnorm(700,0,1); g2<-rnorm(100,-0.5,1.5); g3<-rnorm(100,-1,2.5)     gt<-data.frame(score=c(g1, g2, g3), fac1=factor(rep(c("a", "b", "c"), c(700, 100, 100))), fac2=ordered(rep(c(0,1,2), c(3,13,4))))     gt$score<-with(gt, ifelse(fac2 == 0, score, score-rnorm(1, 0.5, 2)))     gt$score<-with(gt, ifelse(fac2 == 2, score-rnorm(1, 0.5, 2), score))     gt$score<-round(with(gt, ifelse(score>0, score*-1, score)), 1)+4     gt$score<-with(gt, ifelse(score < -4, -4, score))     gt$cov1<-with(gt, score + rnorm(900, sd=40))/40     hist(gt$score)     gt$score2<-with(gt, 4-score+0.0000001) #Gamma distribution can't have 0s (and is positive skewed???)     hist(gt$score2)          glm1<-glm(score2~cov1+fac1*fac2, family="Gamma", data=gt)      This is quite new territory for me.   1\. Is this a reasonable thing to do?   2\. Are there other distributions I might try and compare (exponential perhaps)? _Update:_   After some comments below, I investigated beta regression using the _betareg_ package in R. It gave me skewed residuals:               gt$scorer<-with(gt, (score--4)/(4--4))     gt$scorer<-with(gt, (scorer*(length(scorer)-1)+0.5)/length(scorer))     b1 <- betareg(scorer ~ cov1 + fac1 * fac2, data=gt)     plot(density(resid(b1))) #Strange residuals, even straight lm looks better      So I had a look at a quasibinomial regression and it gave me smaller and better looking residuals:               glm2 <- glm(scorer~cov1 + fac1 * fac2, data=gt, family="quasibinomial")     plot(density(resid(g1))) #Better residuals      Are the residuals good enough to go on in this case?   Or is the fact that d', while based upon T/F, is not a binary variable, a serious issue? _Edit 3: d' clarification_ The below is an example of my d' scores, with the rough distributional qualities and similar raw scores for hits and false positives.               hitrate<-sample(0:16, 100, replace=T, prob=c(rep(0.02,11), 0.025, 0.05, 0.1, 0.2, 0.3, 0.2))/16     hitrate<-ifelse(hitrate==1, 31/32,hitrate); hitrate<-ifelse(hitrate==0, 1/32,hitrate)     farate<-sample(0:32,100, replace=T, prob=c(0.7,0.1,0.05,0.05,0.05,0.02,rep(0.001, 27)))/32     farate<-ifelse(farate==0, 1/64,farate); farate<-ifelse(farate==1, 63/64,farate)          dprime<-round(qnorm(hitrate) - qnorm(farate),1)     plot(density(dprime))


title: Use ARIMA equation outside R

text: I'm using R together with the `forecast` package to set up a ARIMA model, that will be used to predict a energy related variable. I used `auto.arima()` to fit different models (according to geographic region), and I need to put the model coefficients in our database, so that the IT folks can automate things. That's exactly the problem: I simply don't know how set up the equations by looking at the model:               ARIMA(1,0,1)(2,0,1)[12] with non-zero mean           Coefficients:                 ar1     ma1    sar1    sar2     sma1   intercept    prec0    prec1          0.3561  0.3290  0.6857  0.2855  -0.7079  11333.240   15.5291  28.0817          s.e. 0.2079  0.1845  0.2764  0.2251   0.3887   2211.302    6.2147   6.0906      I have 2 regressor variables (prec0 and prec1). Given the residuals, the ARIMA vector `ARIMA(1,0,1)(2,0,1)[12]`, the time series up to period $t$, the number $h$ of forecasting periods and the regressor matrix reg, how can I set a function to return the forecast values? I.e:               do.forecast = function(residuals, ARIMA, timeSeries, h, regMatrix)     {       p = ARIMA[1]       q = ARIMA[3]            ## arima equations here...     }      Thanks! PS: I know this is a possible duplicate of Reproducing ARIMA model outside R, but my model seems very different, and I really don't know how to start with.


title: Next steps after performing a principal component analysis?

text: A few months ago, I developed a questionnaire using a principal component analysis (PCA) and tested the questionnaire for split-half reliability (using a sample which I will call `sample #1`). I am in the process of writing a manuscript to submit for publication, which utilizes the questionnaire and its relationship to depression, anxiety, and stress. I have never developed a questionnaire before, and I don't know what further analyses to conduct with the new sample (`sample #2`). It should be noted that sample #2 is different from sample #1 but not terribly so (both samples are a composite of undergraduate students and individuals who attend exercise classes). I already did the split-half reliability for the questionnaire with sample #2 and it was above .8, but is there anything else that should be done? Do I need to somehow confirm the factor structure of the questionnaire? I did a PCA with sample #2 just to see what would happen, but the items are not lining up with the factors that were specified when I first did the PCA. Frankly, I want to know if I made a mistake in analyzing the questionnaire data for sample #2 with a PCA. Should I have done the second principal components analysis, or is another analysis "better"? Or, should I just stick with the split-half reliability?


title: Taylor Series and Multivariate Delta Method

text: I asked this question on http://math.stackexchange.com/ but did not get any answer. Sorry for cross posting. I'm trying to understand delta method for matrices and vectors to find the variance-covariance matrices for the functions of matrices and vectors. Please see my attempt below. I'm not sure it is right or wrong and especially confused whether it should be $\mathbb{V}\left\\{ \mathbf{\textrm{vec}\left(\mathbf{X}\right)}\right\\} $ or $\mathbb{V}\left\\{ \mathbf{\textrm{vec}\left(\mathbf{X}^{\prime}\right)}\right\\} $. I would highly appreciate experts help in this regard. Thanks If $\mathbf{X}$ is a matrix then a Taylor series expansion of $f\left(\mathbf{X}\right)$ about $\mathbf{X}_{0}$ is given by: \begin{align*} f\left(\mathbf{X}\right) & =f\left(\mathbf{X}_{0}\right)+\left[\left.\frac{\partial f\left(\mathbf{X}\right)}{\partial\mathbf{X}}\right|_{\mathbf{X}_{0}}\right]^{\prime}\left(\mathbf{X}-\mathbf{X}_{0}\right)+\left(\begin{array}{c} \textrm{2-nd and higher}\\\ \textrm{order term} \end{array}\right)\\\ f\left(\mathbf{X}\right) & \approxeq f\left(\mathbf{X}_{0}\right)+\left[\left.\frac{\partial f\left(\mathbf{X}\right)}{\partial\mathbf{X}}\right|_{\mathbf{X}_{0}}\right]^{\prime}\left(\mathbf{X}-\mathbf{X}_{0}\right). \end{align*} The approximated variance of $f\left(\mathbf{X}\right)$ about $\mathbf{X}_{0}$ can be computed as \begin{align*} \mathbb{V}\left[f\left(\mathbf{X}\right)\right] & \approxeq\left[\left.\frac{\partial f\left(\mathbf{X}\right)}{\partial\mathbf{X}}\right|_{\mathbf{X}_{0}}\right]^{\prime}\mathbb{V}\left\\{ \mathbf{\textrm{vec}\left(\mathbf{X}\right)}\right\\} \left[\left.\frac{\partial f\left(\mathbf{X}\right)}{\partial\mathbf{X}}\right|_{\mathbf{X}_{0}}\right] \end{align*} **Example** If $f\left(\mathbf{X}\right)=\mathbf{X}^{-1}$ then \begin{align*} \frac{\partial f\left(\mathbf{X}\right)}{\partial\mathbf{X}} & =\frac{\partial\mathbf{X}^{-1}}{\partial\mathbf{X}}=-\left(\mathbf{X}^{\prime}\otimes\mathbf{X}\right)^{-1} \end{align*} Therefore the approximated variance of $f\left(\mathbf{X}\right)=\mathbf{X}^{-1}$ about $\mathbf{X}_{0}$ is \begin{align*} \mathbb{V}\left(\mathbf{X}^{-1}\right) & \approxeq\left[-\left(\mathbf{X}^{\prime}\otimes\mathbf{X}\right)^{-1}\right]^{\prime}\mathbb{V}\left\\{ \mathbf{\textrm{vec}\left(\mathbf{X}\right)}\right\\} \left[-\left(\mathbf{X}^{\prime}\otimes\mathbf{X}\right)^{-1}\right]\\\ & =\left[\left(\mathbf{X}^{\prime}\otimes\mathbf{X}\right)^{-1}\right]^{\prime}\mathbb{V}\left\\{ \mathbf{\textrm{vec}\left(\mathbf{X}\right)}\right\\} \left[\left(\mathbf{X}^{\prime}\otimes\mathbf{X}\right)^{-1}\right]. \end{align*} If $\mathbf{X}$ is a matrix and $\mathbf{y}$ is a column vector then a Taylor series expansion of $f\left(\mathbf{X},\mathbf{y}\right)$ about $\left(\mathbf{X}_{0},\mathbf{y}_{0}\right)$ is given by: \begin{align*} f\left(\mathbf{X},\mathbf{y}\right) & =f\left(\mathbf{X}_{0},\mathbf{y}_{0}\right)+\left[\left.\frac{\partial f\left(\mathbf{X},\mathbf{y}\right)}{\partial\mathbf{X}}\right|_{\left(\mathbf{X}_{0},\mathbf{y}_{0}\right)}\right]^{\prime}\left(\mathbf{X}-\mathbf{X}_{0}\right)+\left[\left.\frac{\partial f\left(\mathbf{X},\mathbf{y}\right)}{\partial\mathbf{y}}\right|_{\left(\mathbf{X}_{0},\mathbf{y}_{0}\right)}\right]^{\prime}\left(\mathbf{y}-\mathbf{y}_{0}\right)+\left(\begin{array}{c} \textrm{2-nd and higher}\\\ \textrm{order term} \end{array}\right)\\\ f\left(\mathbf{X},\mathbf{y}\right) & \approxeq f\left(\mathbf{X}_{0},\mathbf{y}_{0}\right)+\left[\left.\frac{\partial f\left(\mathbf{X},\mathbf{y}\right)}{\partial\mathbf{X}}\right|_{\left(\mathbf{X}_{0},\mathbf{y}_{0}\right)}\right]^{\prime}\left(\mathbf{X}-\mathbf{X}_{0}\right)+\left[\left.\frac{\partial f\left(\mathbf{X},\mathbf{y}\right)}{\partial\mathbf{y}}\right|_{\left(\mathbf{X}_{0},\mathbf{y}_{0}\right)}\right]^{\prime}\left(\mathbf{y}-\mathbf{y}_{0}\right) \end{align*} The approximated variance of $f\left(\mathbf{X},\mathbf{y}\right)$ about $\left(\mathbf{X}_{0},\mathbf{y}_{0}\right)$ can be computed as \begin{align*} \mathbb{V}\left[f\left(\mathbf{X},\mathbf{y}\right)\right] & \approxeq\left[\left.\frac{\partial f\left(\mathbf{X},\mathbf{y}\right)}{\partial\mathbf{X}}\right|_{\left(\mathbf{X}_{0},\mathbf{y}_{0}\right)}\right]^{\prime}\mathbb{V}\left\\{ \mathbf{\textrm{vec}\left(\mathbf{X}\right)}\right\\} \left[\left.\frac{\partial f\left(\mathbf{X},\mathbf{y}\right)}{\partial\mathbf{X}}\right|_{\left(\mathbf{X}_{0},\mathbf{y}_{0}\right)}\right]+\left[\left.\frac{\partial f\left(\mathbf{X},\mathbf{y}\right)}{\partial\mathbf{y}}\right|_{\left(\mathbf{X}_{0},\mathbf{y}_{0}\right)}\right]^{\prime}\mathbf{\mathbb{V}\left(\mathbf{y}\right)}\left[\left.\frac{\partial f\left(\mathbf{X},\mathbf{y}\right)}{\partial\mathbf{y}}\right|_{\left(\mathbf{X}_{0},\mathbf{y}_{0}\right)}\right] \end{align*} **Example** If $f\left(\mathbf{X},\mathbf{y}\right)=\mathbf{X}\mathbf{y}$ then \begin{align*} \frac{\partial f\left(\mathbf{X},\mathbf{y}\right)}{\partial\mathbf{X}} & =\frac{\partial\mathbf{X}\mathbf{y}}{\partial\mathbf{X}}=\frac{\partial\textrm{vec}\left(\mathbf{X}\mathbf{y}\right)}{\partial\textrm{vec}\left(\mathbf{X}\right)}=\frac{\partial\textrm{vec}\left(\mathbf{I}\mathbf{X}\mathbf{y}\right)}{\partial\textrm{vec}\left(\mathbf{X}\right)}\\\ & =\frac{\partial\left(\mathbf{y}^{\prime}\otimes\mathbf{I}\right)\textrm{vec}\left(\mathbf{X}\right)}{\partial\textrm{vec}\left(\mathbf{X}\right)}=\left(\mathbf{y}^{\prime}\otimes\mathbf{I}\right)\frac{\partial\textrm{vec}\left(\mathbf{X}\right)}{\partial\textrm{vec}\left(\mathbf{X}\right)}\\\ & =\left(\mathbf{y}^{\prime}\otimes\mathbf{I}\right)\qquad\textrm{and}\\\ \frac{\partial f\left(\mathbf{X},\mathbf{y}\right)}{\partial\mathbf{y}} & =\frac{\partial\mathbf{X}\mathbf{y}}{\partial\mathbf{y}}=\mathbf{X}. \end{align*} Therefore the approximated variance of $f\left(\mathbf{X},\mathbf{y}\right)=\mathbf{X}\mathbf{y}$ about$\left(\mathbf{X}_{0},\mathbf{y}_{0}\right)$ is \begin{align*} \mathbb{V}\left(\mathbf{X}\mathbf{y}\right) & \approxeq\left(\mathbf{\mathbf{y}_{0}}^{\prime}\otimes\mathbf{I}\right)^{\prime}\mathbb{V}\left\\{ \mathbf{\textrm{vec}\left(\mathbf{X}\right)}\right\\} \left(\mathbf{\mathbf{y}_{0}}^{\prime}\otimes\mathbf{I}\right)+\mathbf{\mathbf{X}}_{0}^{\prime}\mathbf{\mathbb{V}\left(\mathbf{y}\right)}\mathbf{\mathbf{X}}_{0}. \end{align*} **Edited** Let $\mathbf{A}$ be $m\times n$ matrix with $\mathbf{a}_{i}(i=1,\ldots,n)$ as the \emph{$i$}-th column vector. The vectorization operation , $\textrm{vec}\left(.\right)$, is an operation from $\mathcal{R}^{m\times n}$ to $\mathcal{R}^{mn}$, with $ \textrm{vec}(\mathbf{A})=\textrm{vec}\left(\begin{bmatrix}a_{11} & \ldots & a_{1n}\\\ \vdots & \ddots & \vdots\\\ a_{m1} & \ldots & a_{mn} \end{bmatrix}\right)=\textrm{vec}\left(\left[\begin{array}{ccc} \mathbf{a}_{1} & \cdots & \mathbf{a}_{n}\end{array}\right]\right)=\left[\begin{array}{c} \mathbf{a}_{1}\\\ \vdots\\\ \mathbf{a}_{n} \end{array}\right]. $


title: Central moments of a gaussian mixture density?

text: Given the pdf $f(x) = \sum_i \omega_i \mathcal{N}(x; \mu_i, C_i )$ of a gaussian mixture density, where the $i$-th component has mean $\mu_i$ and covariance matrix $C_i$ and the weights $\omega_i$ sum to 1, is there a formula to compute the moments and central moments of this density from the given $\omega_i$'s, $\mu_i$'s and $C_i$'s?


title: Choice of default bin in categorical dummy data changes coefficients in Logistic Regression

text: I have 8 mutually exclusive categorical predictors (i.e., only one can have a value of 1 at a given time) and I built a logistic regression model using 7 (i.e. n-1) of these. When I leave out the last bin, I get both positive & negative coefficients. If I leave out the first bin, all my coefficients become negative. Most of them continue to remain statistically significant. Here is my Q: 1\. If all coefficients are negative, does it say something about the model? It is like saying there are negative odds for the event in all the Xs! 2\. How come changing the default bin (which is left out) changes the coefficients so radically? Is it because of residuals, or am I missing something. Please note that this is a rare event regression (1700 events in 187k observations)


title: Finding overall survival in competing risk model

text: I am dealing with time until event data and there are three types of events that I am tracking Type A Death Type B Death Event C. Type B Death cannot happen unless C happened. At any given time every person is at risk to type A and type B death. When somebody dies, it is difficult to say what type of death they had: they are at risk to both and have a nonzero probability of dying because of B but it may turn out that A kills them first and we are never able to see the death due to B I am interested in finding the overall survival $P($Survive$)$ which I have thought of as $(1 - P(A))(1 - P(B))$. I wonder if I should think of $P($Survive$)$ as 1 - P(any type of death) and P(any type of death) = P(A) + P(B) - P(A and B) If they are competing risks then does this quantity have to be 0? If A and B aren't physical deaths, then just events, can this quantity be nonzero? I have two separate CoxPh models to model the hazard of A and B separately. To complicate things further, what if A can be followed by B? Normally we would have two models where the model for A would treat every observation as censored if it experiences B. Would it better if, for the times that B follows A, we treat that individual as having both events A and B? Since this person is counted twice we can deal with that by subtracting out P(A and B)?


title: Stata's -margins- shows negative predicting value for count model

text: I am modeling "Innovation" counts over "HHI" and "squared HHI" with fixed effect poisson model. (My sample is unbalanced, small & wide one.) I am aware of that I can use Stata's -margins- command to obtain the estimated innovation count based on observable values. Below are my codes and results:                   xtpoisson innovation hhi c.hhi#c.hhi fee age c.age#c.age size c.size#c.size, fe         margins, dydx(*)         margins, at(hhi=0.0001(0.005)0.100)      The estimated coefficient for HHI is 32.80 and for squared HHI is -131.31. The "Average Marginal Effect" of HHi is 25.37. However, when I used the third line to get the "marginal effect at representative values," all predicted values are negative. I am surprised because I thought I was modeling some positive count number and I don't understand why I got negative predicted values, e.g.,-3.xx innovations, which is generated from the sample mean of all dependent variables. My questions may be three-fold: 1) Do I understand -margins- results correctly? Or I missed something? 2) Is there any coding mistakes that may result to what I got now? 3) If there is no coding mistakes, does the result imply that my fixed effect poisson model have some problems? A bit more information: when I plot the observed innovation over HHI fit with quadratic line, the curve is u-shape and above zero. However, when I use the -marginscontplot- command to plot the predicted innovation over HHI, it's an inverted u-shape and below zero. Thank you for your help!


title: interaction of categorical and continuous variables

text: I have a dependent variable that is continuous and I have two independent variables: one continuous and one categorical (with 2 categories) The interaction between the independent variables is significant. Which statistical analysis should I use (in R) to proceed with the analysis and document the interaction? (Should I simply analyze each of the two categories separately using simple linear regression?)


title: Nested Anova Question:

text: I've a dataset for measurement of protein marker abundance for two groups of individuals, Treatment and Control. Within these two groups there are several subjects and for each subject there are three replicates. The number of subjects in these Treatment and Control groups are unequal with 4 individuals for treatment and 10 for Control. So overall,               :Control and Treatment groups(Denoted by 1 and 2 respectively)         :10 individuals in Control(A to J), 4 individuals in Treatment(K to N)             :3 replicates for each indiv. in Control, 3 replicates for each indiv. in Treatment      I would like to perform a nested ANOVA on this data however I'm new to this kind of data and would like to know the followings:   * Is the 'Replicate' column is redundant here and can I simply drop/ignore it?   * Can I use the following ANOVA formula in R? Is there anything special I need to take care of because of unequal numbers of control and treated individuals? protein.aov <\- aov(abundance ˜ treatgroup + Error(treatgroup %in% individual), data=protein.df)   * Can I do a post hoc pairwise comparison like a TukeyHSD ? The data is as Follows:               treatgroup  individual  replicate   abundance     1   A   R1  -0.709258936     1   A   R2  -0.54767131     1   A   R3  -0.907607661     1   B   R1  3.729646649     1   B   R2  -0.650402382     1   B   R3  0.884222978     1   C   R1  -4.443417184     1   C   R2  3.709624162     1   C   R3  3.076829126     1   D   R1  2.109771383     1   D   R2  4.950350294     1   D   R3  -1.162741304     1   E   R1  0.799105402     1   E   R2  2.929226412     1   E   R3  2.95692962     1   F   R1  2.011646397     1   F   R2  -3.4757793     1   F   R3  4.615843439     1   G   R1  7.324129703     1   G   R2  -6.56365647     1   G   R3  5.848340873     1   H   R1  7.375089916     1   H   R2  -0.709544581     1   H   R3  -1.715528803     1   I   R1  6.70860394     1   I   R2  -4.325520039     1   I   R3  7.623999717     1   J   R1  1.959268861     1   J   R2  3.794791979     1   J   R3  -0.443267523     2   K   R1  4.489545974     2   K   R2  -0.93677524     2   K   R3  8.030252255     2   L   R1  -0.133320899     2   L   R2  3.802649555     2   L   R3  1.118932954     2   M   R1  2.054437925     2   M   R2  -3.872643548     2   M   R3  5.695342112     2   N   R1  -4.913796298     2   N   R2  4.647048982     2   N   R3  6.729868259


title: Analyzing the shape of distribution on a histogram

text: I have three different histograms which are generated from one sample. In each of the histograms, both variables are the same. Although binwidths for each histogram is different. By looking at each of them, we can have different interpretation of the data that made these histogram. One of them shows a bimodal distribution shape, another shows a skew distribution. I want to have the opinion of other people for how do they analyze the distribution of each of histograms by their shape. ![Dist](http://oi60.tinypic.com/210j2c1.jpg)


title: Why my MA plot shows a pattern/trend while the normalized expression values not?

text: I am doing microarray analysis using bioconductor. I first performed the quality checks which was ok according to recommendations. I created an MA plot after normalizing the AffyBatch data. My MA plot looks as follows:               mva.pairs(exprs(myEset),col=2:7)       ![MA plot after normalization](http://i.stack.imgur.com/V7oGn.png) Then after, I performed fit using the following commands.               library(limma)     Group<- factor(pData(ABatch)[,'group'] , levels = levels(pData(ABatch)[,'group']))     design<- model.matrix(~Group)     fit <-lmFit(myEset,design)     fit<-eBayes(fit)       top.all <- topTable(fit,n=nrow(exprs(myEset)),adjust="BH",coef=2)      Now again I am trying to create MA plot from `top.all` define above using following commands               M <- top.all$logFC #log fold-change         A <- top.all$AveExpr # ave. expr. level     mat <- cbind(A,M)                plotMA(mat,main="MA plot")                        abline(h=3,col="red")      I get the following plot ![MA plot after fit](http://i.stack.imgur.com/jeLUA.png) Now, I am not sure why the second plot is showing a pattern (fold-change increases as the avg. expression increases). Is there any problem in the analysis here?


title: Post-hoc chi-squared multiple comparison tests using complex samples module in SPSS v. 20

text: This questions refers to computing cross-tabulations and chi-squared tests of independence through the SPSS complex samples module available in SPSS versions 19 and up. I am able to compute the overall chi-squared through the test of independence option, but I'd like to compute post-hoc tests of independence between all cells. Does anyone know if this is even possible (i.e., post-hoc tests of independence through the complex samples module of SPSS)?


title: Which machine learning algorithm should I use?

text: I have one dependent binary categorical variable, and one independent continuous variable. There is a lot of randomness deciding the result of the dependent variable. The relationship between the independent variable and the dependent variable is linear. I have 2,000 data points to train data on. Some possibilities are:   * Logistic regression - simplest option   * SVM (support vector machines)   * Naive bayes   * Random forests - I see this does well on kaggle, but I have a simple one variable linear relationship, so it seems random trees isn't necessary here.


title: Random forest and prediction

text: I am trying to understand how Random Forest works. I have a grasp about how trees are build but can not understand how Random Forest make predictions on out of bag sample. Could anyone give me a simple explanation, please?:)


title: Smart way to search through a very large parameter space

text: I have a system whose performance is based on a rather large parameter set (200 parameters, lets say, of which each can take a very wide range of values). There are tests to evaluate the performance using a given parameter set, but it is expensive computationally, and takes time. Considering that, there may or may not be patterns relating how the parameters change and how the performance changes, is there a smart way to search this rather huge parameter space for optimum values that maximize the performance of the system, and minimize the computational overhead? A trend MAY be guessed for some of those parameters, but not with much confidence. A test range may be specified by the user, for each of the parameters. Is there any other way that is definitely better than the brute force method of doing an exhaustive search? Something that looks for trends on-the-fly and avoids searching in a direction, maybe? But, it needs good statistical support to ignore a branch. A subset of the parameter space may contain dependent variables, but we do not know any information about that. So, no assumptions there.


title: Path diagram for a second order factor model

text: I want to fit a second order factor model to some data. Do I need to add an error term to the endogenous factor? I'm thinking about e10, here. I omitted it at first, but AMOS was surprised about that and questioned my model. What's the usual way to set up a two-level factor model? ![A path diagram from AMOS](http://i.stack.imgur.com/vTzyq.jpg)


title: Consistency of the extremum estimators

text: I am trying to understand intuitively the assumptions mentioned in the following theorem: _Theorem._ If there is a function $Q_0(\theta)$ s.t. (i) $Q_0(\theta)$ is uniquely maximized at $\theta_0$; (ii) $\Theta$ is compact; (iii) $Q_0(\theta)$ is continuous; (iv) $\hat{Q}_n(\theta)$ converges uniformly in probability to $Q_0(\theta)$, then $\hat{\theta} \rightarrow \theta_0$ in probability. I understand that the last assumption says that $\hat{Q}_n(\theta)$ is between $Q_0(\theta) - \varepsilon$ and $Q_0(\theta) + \varepsilon$. for some $\varepsilon$. Thus, given assumption (i) we can find $\theta_L < \theta_0 < \theta_U$ s.t. $\hat{\theta} \in (\theta_L, \theta_U)$. Does assumption (i) guarantee that $|\theta_0| \neq \infty$? Why intuitively assumptions (ii) and (iii) are needed?


title: Is lm.morantest valid on residuals from geeglm fit of binary data?

text: Is it valid to apply lm.morantest (in the package spdep) to test for spatial autocorrelation among residuals from a Generalized Explanatory Equation (geeglm in geepack) model fit of binary data? If not, is there another solution to test for residual spatial autocorrelation since Moran's is not appropriate for the binary raw response data either?


title: Open, multi-class , active learning classfiers

text: I am trying to classify text documents using a huge corpora. Thats a huge tagset (more than 1000 tags). Corpus will have 1000 samples for each tag. But the tagset is not closed. New tags can be added in future as per the text content. There will be human intervention in between where I can get new tags if the situation arise. There can be multiple tags per text document. nltk book talks about it in chapter 6. They referred to it as open-class classifiers in section 6.1 after the three bullet points. Is it possible via some already implemented toolkit in python/ R/ C/ Java? If not I guess it is a problem of semi-supervised learning. I am familiar with nltk as well as scikit classifiers. Can those functions be used with other reinforcement learning toolkit? I have no experience with reinforcement/ semi-supervised learning. I hope I am clear in the question. I will be thankful for your suggestions. Thanks.


title: Is it wrong to use line plots for discrete data?

text: I've often seen discrete datasets plotted as line plots, but it occurs to me that the line infers a value at a point between the measurement intervals which is meaningless for discrete datasets. Is it therefore the case the using line plots for discrete data is wrong? As an example, take two time-series datasets, one continuous (my weight, measured daily in the morning) and one discrete (the number of donuts I eat per day). It makes sense for the first dataset to be a line plot, as it is reasonable to infer that my weight in any given afternoon will be related to my weight the preceding and following mornings. However, if the number of donuts is represented as a line graph the lines between the dots no meaning can be inferred from that line. EDIT Here's another example: the Federal Hourly Minimum Wage Since Its Inception plot at http://mste.illinois.edu/courses/ci330ms/youtsey/lineinfo.html Unless I'm mistaken, the minimum wage changes are discrete, and hence it's not possible to look up some arbitrarily selected time and establish the minimum wage at the point using the line interconnecting the dots.


title: How to handle with missing values in order to prepare data for feature selection with LASSO?

text: My situation:   * small sample size: 116    * binary outcome variable   * long list of explanatory variables: 44   * explanatory variables did not come from the top of my head; their choice was based on the literature.   * most cases in the sample and most variables have missing values. Approach to feature selection chosen: LASSO R's glmnet package won't let me run the glmnet routine, apparently due to the existence of missing values in my data set. There seems to be various methods for handling missing data, so I would like to know:   * Does LASSO impose any restriction in terms of the method of imputation that I can use?   * What would be the best bet for imputation method? Ideally, I need a method that I could run on SPSS (preferably) or R. UPDATE1: It became clear from some of the answers below that I have do deal with more basic issues before considering imputation methods. I would like to add here new questions regarding that. On the the answer suggesting the coding as constant value and the creation of a new variable in order to deal with 'not applicable' values and the usage of group lasso:   * Would you say that if I use group LASSO, I would be able to use the approach suggested to continuous predictors also to categorical predictors? If so, I assume it would be equivalent to creating a new category - I am wary that this may introduce bias.   * Does anyone know if R's glmnet package supports group LASSO? If not, would anyone suggest another one that does that in combination with logistic regression? Several options mentioning group LASSO can be found in CRAN repository, any suggestions of the most appropriate for my case? Maybe SGL? This is a follow-up on a previous question of mine (How to select a subset of variables from my original long list in order to perform logistic regression analysis?). OBS: I am not a statistician.


title: How to conduct the two sample test with extra data to eliminate between-subjects variability

text: I am stuck with the correlated and independent data combined in one study. Here's my dilemma: Say X is a drug(explanatory variable) and Y is a gene expression(response variable). Normally, you would give out drugs to a half of your group(randomly chosen) and placebo to the rest, measure the gene expression for each person in the group, and conduct the differential expression analysis between two conditions using a standard package such as deseq. However, let's say between-subjects variability is huge, so you also measured the gene expression prior to taking the pill as well. In other words, each subject has two data points(before and after the taking the pill). How do you incorporate this information into the analysis? I tried to subtract the 'before' value from after 'value' for each patient, and conduct t-test on those statistics between conditions(drug & placebo), but I am not sure if I am doing it right. I am curious if there is a standard, already established way to do this? You can assume Y follows Poisson or negative binomial.


title: Logistic Regression doesn't converge if I use a particular baseline

text: Apologies if the question is too broad, then please close the thread. Or maybe this belongs in StackOverflow? I observed something strange and wonder if someone more educated can shed some light on the issue. I don't think it's necessary to go into a ton of detail (i.e., providing code and data). I'm running a logistic regression in R using glm. I have some dataset where about 90% of my dependendent variables = 1. There are also 15 independent variables. About 10 of the variables are factors with 3 levels: "Yes", "No", "No Data" (not NA, to prevent exclusion). Overwhelmingly most of the data (80%+) for these factor variables is "Yes", but the default baseline for those independent variables was "No Data." With this configuration, the regression converges. Also, if I use "No" as the baseline, the regression converges. If I change the baseline for those variables to "Yes," the regression no longer converges. Any particular reason for this? Is this some property of logistic regression or a result of R code? Is it because my dependent variable has low variance or because the independent variables? I'm curious because ideally we would like faster convergance, and if I knew some properties of logistic regression that prevent convergance, then I could build better regressions and be more confident about the results. For example, I had a similar dataset, and made the same changes (releveling). This regression happened to converge, but it took much more iterations, and some of the results were actually different than before I made the change. So, it seems like important knowledge. Thanks!


title: How to find out some particular distribution given the grouped data and a polynomial fitted to the data

text: I have to analyse a set of grouped data.The data is divided into groups by some categories for example: BP(<=60), BP((60,80]); Pulse(<75), Pulse((75,90]) & Pulse(>90) etc having many more factors affecting health.Each of these categories have no. of people as their frequency.[Example has been given below].I have managed to identify the distribution of some categories (for example for BP(60,80] -Normal(10.8,4.43) fits well), but some of the categories I am not able to identify and fit any statistical/Standard distributions to them (for example: BP(<=60)). I have fitted polynomials to all the categories individually and the R-squared value of the polynomial fitted is more than 0.95 for almost all the categories. How can I use these polynomials to find out about the distributions of the data? Also is there any other way to find out the distributions even if I don't use these polynomials? Is there any function in R which can help me to analyse the data using these good fitted polynomials? ### Interval[BP(<=60) Freq Relative Freq Interval BP((60,80]) Freq Relative Freq   1. 6.5                      0              0                              7                   0            0   2. 8.5                      8          0.296296296               9                   7       0.23333333   3. 10.5                    9          0.333333333              11                10       0.333333333   4. 12.5                    4          0.148148148              13                 7        0.233333333   5. 14.5                    3          0.111111111              15                 6        0.2   6. 16.5                    3          0.111111111              17                 0             0   7. 18.5                    0            0                             Sum:              30           1   8. Sum:                 27            1


title: Speech recognition, words out of dictionary

text: I'm performing word recognition by using a tradional procedure. I'm extracting MFCC features. Then I'm creating a code book in order to do vector quantization. After that, I train discrete HMM for two words: 1stWrod, 2dWord. Until now, I have been performing the classification like this: I estimate the probability in the two trained models of a new audio segment with the proper feature extraction and quantization. I say the audio correspond to the class with the highest probability. That gives me good results. But any audio segment is classified as any of those tho words, when sometimes it's not. I don't not how to say that is not corresponding to any class. I'm not sure if I could solve that by training another model with all of the other data, because it's very different and I think the model wouldn't be enough.


title: Why is one parameter estimate so high in logistic regression?

text: I am doing logistic regression on a model with a dependent variable of 4 different sizes of fish. I originally tried to do ordinal logistic but I ended up binning responses into "small" and "large" for various reasons. I am now running binomial logistic. The problem is that the parameter estimate of one of my independent variables with values that range from .16 to 2.0 is VERY high (7.5 +- 2.7). I have this problem whether I am doing ordinal or binomial, with the same independent variable. What could be causing this?


title: I want to show$X : ((0,1],\mathcal B,\lambda)\to \mathbb R$ is random variable

text: Let $F:\mathbb R \to [0,1]$ be a distribution function of a probability measure $P$ $(i.e.,F(x)=P((-\infty,x])) $. Then show that There is a random variable $X : ((0,1],\mathcal B,\lambda)\to \mathbb R$,(where $\mathcal B$ is the borel $\sigma $-algebra and $\lambda$ is Lebesgue measure) ,such that $P_{X}=P$


title: How to show $P_{\hat X}=P_{X}$, where $P_{X}$ is a known distribution?

text: Let $X$ be a random variable on the probability space $(\Omega,\mathcal B,P)$, with distribution $P_{X}$. Consider the random variable $\hat X$ on the probability space $(\mathbb R,\mathcal B_{\mathbb R},P_{X})$,defined by $\hat X(x)=x$ . Then $P_{\hat X}=P_{X}$. $\mathcal B $ is $\sigma$-algebra.