title: How to use profile likelihood?

text: I am new to profile likelihood and do not really understand what advantages it may have. Lets say I have the following results estimating the means of three groups. What can I say about them? ![enter image description here](http://i.stack.imgur.com/1USWx.png) **R code:**               profile.likelihood<-function(dat, muVals){       likVals <- sapply(muVals,                         function(mu){                           (sum((dat - mu)^2) /                              sum((dat - mean(dat))^2)) ^ (-length(dat)/2)                         }       )          return(cbind(muVals,likVals))     }                 muVals <- seq(0,20, length = 10000)     dat1<-runif(10,0,20)     a<-profile.likelihood(dat1,muVals)     dat2<-runif(10,0,20)     b<-profile.likelihood(dat2,muVals)     dat3<-runif(10,8,20)     c<-profile.likelihood(dat3,muVals)          plot(a, type="l", lwd=4,     ylab="Likelihood", xlab="Score"     )     lines(b, lty=2, lwd=4)     lines(c, lty=3, lwd=4)      **EDIT:** Here is data more similar to my actual. 1) data can take any value from 0 to 20 2)I know from a relatively large number of previous results that the distribution is not normal and looks like the top plot. My data is shown in the boxplots in the lower panel. I wish to estimate the means of each group and compare them using likelihoods or know why I should not do this. ![enter image description here](http://i.stack.imgur.com/65kEj.png) **dput() of new data:**               structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2,      2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 13.6986301369863,      16.1643835616438, 12.0547945205479, 12.1722113502935, 9.74559686888454,      0.430528375733855, 11.3502935420744, 10.6457925636008, 9.9412915851272,      10.7240704500978, 10.958904109589, 11.6242661448141, 17.9701232444495,      15.9326690901071, 7.98247314058244, 14.4031004607677, 13.5198221541941,      2.82421704847366, 16.114045586437, 19.328767512925, 3.74181577004492,      17.4085859861225, 19.2017483590171, 8.26946665905416, 10.0207103956491,      16.1689247898757, 13.9989542039111, 9.80047978740185, 17.8596440777183,      18.1706223106012, 18.1891529858112, 11.7567204562947), .Dim = c(32L,      2L))      **dput() of prior density**               structure(list(x = c(0, 0.0391389432485323, 0.0782778864970646,      0.117416829745597, 0.156555772994129, 0.195694716242661, 0.234833659491194,      0.273972602739726, 0.313111545988258, 0.352250489236791, 0.391389432485323,      0.430528375733855, 0.469667318982387, 0.50880626223092, 0.547945205479452,      0.587084148727984, 0.626223091976517, 0.665362035225049, 0.704500978473581,      0.743639921722113, 0.782778864970646, 0.821917808219178, 0.86105675146771,      0.900195694716243, 0.939334637964775, 0.978473581213307, 1.01761252446184,      1.05675146771037, 1.0958904109589, 1.13502935420744, 1.17416829745597,      1.2133072407045, 1.25244618395303, 1.29158512720157, 1.3307240704501,      1.36986301369863, 1.40900195694716, 1.44814090019569, 1.48727984344423,      1.52641878669276, 1.56555772994129, 1.60469667318982, 1.64383561643836,      1.68297455968689, 1.72211350293542, 1.76125244618395, 1.80039138943249,      1.83953033268102, 1.87866927592955, 1.91780821917808, 1.95694716242661,      1.99608610567515, 2.03522504892368, 2.07436399217221, 2.11350293542074,      2.15264187866928, 2.19178082191781, 2.23091976516634, 2.27005870841487,      2.30919765166341, 2.34833659491194, 2.38747553816047, 2.426614481409,      2.46575342465753, 2.50489236790607, 2.5440313111546, 2.58317025440313,      2.62230919765166, 2.6614481409002, 2.70058708414873, 2.73972602739726,      2.77886497064579, 2.81800391389432, 2.85714285714286, 2.89628180039139,      2.93542074363992, 2.97455968688845, 3.01369863013699, 3.05283757338552,      3.09197651663405, 3.13111545988258, 3.17025440313112, 3.20939334637965,      3.24853228962818, 3.28767123287671, 3.32681017612524, 3.36594911937378,      3.40508806262231, 3.44422700587084, 3.48336594911937, 3.52250489236791,      3.56164383561644, 3.60078277886497, 3.6399217221135, 3.67906066536204,      3.71819960861057, 3.7573385518591, 3.79647749510763, 3.83561643835616,      3.8747553816047, 3.91389432485323, 3.95303326810176, 3.99217221135029,      4.03131115459883, 4.07045009784736, 4.10958904109589, 4.14872798434442,      4.18786692759295, 4.22700587084149, 4.26614481409002, 4.30528375733855,      4.34442270058708, 4.38356164383562, 4.42270058708415, 4.46183953033268,      4.50097847358121, 4.54011741682975, 4.57925636007828, 4.61839530332681,      4.65753424657534, 4.69667318982387, 4.73581213307241, 4.77495107632094,      4.81409001956947, 4.853228962818, 4.89236790606654, 4.93150684931507,      4.9706457925636, 5.00978473581213, 5.04892367906066, 5.0880626223092,      5.12720156555773, 5.16634050880626, 5.20547945205479, 5.24461839530333,      5.28375733855186, 5.32289628180039, 5.36203522504892, 5.40117416829746,      5.44031311154599, 5.47945205479452, 5.51859099804305, 5.55772994129159,      5.59686888454012, 5.63600782778865, 5.67514677103718, 5.71428571428571,      5.75342465753425, 5.79256360078278, 5.83170254403131, 5.87084148727984,      5.90998043052838, 5.94911937377691, 5.98825831702544, 6.02739726027397,      6.0665362035225, 6.10567514677104, 6.14481409001957, 6.1839530332681,      6.22309197651663, 6.26223091976517, 6.3013698630137, 6.34050880626223,      6.37964774951076, 6.4187866927593, 6.45792563600783, 6.49706457925636,      6.53620352250489, 6.57534246575342, 6.61448140900196, 6.65362035225049,      6.69275929549902, 6.73189823874755, 6.77103718199609, 6.81017612524462,      6.84931506849315, 6.88845401174168, 6.92759295499021, 6.96673189823875,      7.00587084148728, 7.04500978473581, 7.08414872798434, 7.12328767123288,      7.16242661448141, 7.20156555772994, 7.24070450097847, 7.27984344422701,      7.31898238747554, 7.35812133072407, 7.3972602739726, 7.43639921722114,      7.47553816046967, 7.5146771037182, 7.55381604696673, 7.59295499021526,      7.6320939334638, 7.67123287671233, 7.71037181996086, 7.74951076320939,      7.78864970645793, 7.82778864970646, 7.86692759295499, 7.90606653620352,      7.94520547945205, 7.98434442270059, 8.02348336594912, 8.06262230919765,      8.10176125244618, 8.14090019569472, 8.18003913894325, 8.21917808219178,      8.25831702544031, 8.29745596868884, 8.33659491193738, 8.37573385518591,      8.41487279843444, 8.45401174168297, 8.49315068493151, 8.53228962818004,      8.57142857142857, 8.6105675146771, 8.64970645792564, 8.68884540117417,      8.7279843444227, 8.76712328767123, 8.80626223091977, 8.8454011741683,      8.88454011741683, 8.92367906066536, 8.96281800391389, 9.00195694716243,      9.04109589041096, 9.08023483365949, 9.11937377690802, 9.15851272015655,      9.19765166340509, 9.23679060665362, 9.27592954990215, 9.31506849315068,      9.35420743639922, 9.39334637964775, 9.43248532289628, 9.47162426614481,      9.51076320939335, 9.54990215264188, 9.58904109589041, 9.62818003913894,      9.66731898238747, 9.70645792563601, 9.74559686888454, 9.78473581213307,      9.8238747553816, 9.86301369863014, 9.90215264187867, 9.9412915851272,      9.98043052837573, 10.0195694716243, 10.0587084148728, 10.0978473581213,      10.1369863013699, 10.1761252446184, 10.2152641878669, 10.2544031311155,      10.293542074364, 10.3326810176125, 10.3718199608611, 10.4109589041096,      10.4500978473581, 10.4892367906067, 10.5283757338552, 10.5675146771037,      10.6066536203522, 10.6457925636008, 10.6849315068493, 10.7240704500978,      10.7632093933464, 10.8023483365949, 10.8414872798434, 10.880626223092,      10.9197651663405, 10.958904109589, 10.9980430528376, 11.0371819960861,      11.0763209393346, 11.1154598825832, 11.1545988258317, 11.1937377690802,      11.2328767123288, 11.2720156555773, 11.3111545988258, 11.3502935420744,      11.3894324853229, 11.4285714285714, 11.46771037182, 11.5068493150685,      11.545988258317, 11.5851272015656, 11.6242661448141, 11.6634050880626,      11.7025440313112, 11.7416829745597, 11.7808219178082, 11.8199608610568,      11.8590998043053, 11.8982387475538, 11.9373776908023, 11.9765166340509,      12.0156555772994, 12.0547945205479, 12.0939334637965, 12.133072407045,      12.1722113502935, 12.2113502935421, 12.2504892367906, 12.2896281800391,      12.3287671232877, 12.3679060665362, 12.4070450097847, 12.4461839530333,      12.4853228962818, 12.5244618395303, 12.5636007827789, 12.6027397260274,      12.6418786692759, 12.6810176125245, 12.720156555773, 12.7592954990215,      12.7984344422701, 12.8375733855186, 12.8767123287671, 12.9158512720157,      12.9549902152642, 12.9941291585127, 13.0332681017613, 13.0724070450098,      13.1115459882583, 13.1506849315068, 13.1898238747554, 13.2289628180039,      13.2681017612524, 13.307240704501, 13.3463796477495, 13.385518590998,      13.4246575342466, 13.4637964774951, 13.5029354207436, 13.5420743639922,      13.5812133072407, 13.6203522504892, 13.6594911937378, 13.6986301369863,      13.7377690802348, 13.7769080234834, 13.8160469667319, 13.8551859099804,      13.894324853229, 13.9334637964775, 13.972602739726, 14.0117416829746,      14.0508806262231, 14.0900195694716, 14.1291585127202, 14.1682974559687,      14.2074363992172, 14.2465753424658, 14.2857142857143, 14.3248532289628,      14.3639921722113, 14.4031311154599, 14.4422700587084, 14.4814090019569,      14.5205479452055, 14.559686888454, 14.5988258317025, 14.6379647749511,      14.6771037181996, 14.7162426614481, 14.7553816046967, 14.7945205479452,      14.8336594911937, 14.8727984344423, 14.9119373776908, 14.9510763209393,      14.9902152641879, 15.0293542074364, 15.0684931506849, 15.1076320939335,      15.146771037182, 15.1859099804305, 15.2250489236791, 15.2641878669276,      15.3033268101761, 15.3424657534247, 15.3816046966732, 15.4207436399217,      15.4598825831703, 15.4990215264188, 15.5381604696673, 15.5772994129159,      15.6164383561644, 15.6555772994129, 15.6947162426614, 15.73385518591,      15.7729941291585, 15.812133072407, 15.8512720156556, 15.8904109589041,      15.9295499021526, 15.9686888454012, 16.0078277886497, 16.0469667318982,      16.0861056751468, 16.1252446183953, 16.1643835616438, 16.2035225048924,      16.2426614481409, 16.2818003913894, 16.320939334638, 16.3600782778865,      16.399217221135, 16.4383561643836, 16.4774951076321, 16.5166340508806,      16.5557729941292, 16.5949119373777, 16.6340508806262, 16.6731898238748,      16.7123287671233, 16.7514677103718, 16.7906066536204, 16.8297455968689,      16.8688845401174, 16.9080234833659, 16.9471624266145, 16.986301369863,      17.0254403131115, 17.0645792563601, 17.1037181996086, 17.1428571428571,      17.1819960861057, 17.2211350293542, 17.2602739726027, 17.2994129158513,      17.3385518590998, 17.3776908023483, 17.4168297455969, 17.4559686888454,      17.4951076320939, 17.5342465753425, 17.573385518591, 17.6125244618395,      17.6516634050881, 17.6908023483366, 17.7299412915851, 17.7690802348337,      17.8082191780822, 17.8473581213307, 17.8864970645793, 17.9256360078278,      17.9647749510763, 18.0039138943249, 18.0430528375734, 18.0821917808219,      18.1213307240705, 18.160469667319, 18.1996086105675, 18.238747553816,      18.2778864970646, 18.3170254403131, 18.3561643835616, 18.3953033268102,      18.4344422700587, 18.4735812133072, 18.5127201565558, 18.5518590998043,      18.5909980430528, 18.6301369863014, 18.6692759295499, 18.7084148727984,      18.747553816047, 18.7866927592955, 18.825831702544, 18.8649706457926,      18.9041095890411, 18.9432485322896, 18.9823874755382, 19.0215264187867,      19.0606653620352, 19.0998043052838, 19.1389432485323, 19.1780821917808,      19.2172211350294, 19.2563600782779, 19.2954990215264, 19.3346379647749,      19.3737769080235, 19.412915851272, 19.4520547945205, 19.4911937377691,      19.5303326810176, 19.5694716242661, 19.6086105675147, 19.6477495107632,      19.6868884540117, 19.7260273972603, 19.7651663405088, 19.8043052837573,      19.8434442270059, 19.8825831702544, 19.9217221135029, 19.9608610567515,      20), y = c(0.0217667400980039, 0.0227118960613568, 0.0236399642746827,      0.024547462294379, 0.0254284817500052, 0.0262697990497476, 0.0270761925341352,      0.0278447201687443, 0.0285725349842736, 0.0292431898522125, 0.0298599237216537,      0.0304262344415342, 0.0309405476728982, 0.0313967502030459, 0.0317805594200756,      0.0321087524689367, 0.0323813637754389, 0.0325985887432935, 0.0327456085503099,      0.0328338748111955, 0.032871094252979, 0.032858917719392, 0.0327952170079012,      0.0326741710852976, 0.0325130172841641, 0.0323142824586605, 0.0320805363813515,      0.0318066450090944, 0.0315032349150802, 0.0311762923391493, 0.0308283666713486,      0.0304606397730375, 0.0300742685852501, 0.0296767520098838, 0.0292699115474762,      0.0288554692880644, 0.0284340054317109, 0.0280086390955194, 0.0275806684358679,      0.0271507171781247, 0.0267191796498216, 0.0262860885684517, 0.0258516152315375,      0.025415436557538, 0.0249771598811644, 0.0245350861062508, 0.0240885560431356,      0.023637102870035, 0.023179973519991, 0.0227156094045636, 0.022241184162465,      0.0217581061784532, 0.0212659006410898, 0.0207641513896269, 0.020249460771043,      0.0197238756753159, 0.0191886405031699, 0.0186441419420248, 0.0180900824959544,      0.0175265997231419, 0.0169574756929679, 0.0163840336778012, 0.0158076842167135,      0.0152307532198712, 0.0146567117250279, 0.0140876768708212, 0.0135257717396179,      0.0129746875992836, 0.0124411735382691, 0.0119248408581489, 0.0114281180660412,      0.0109534391260541, 0.0105125675916208, 0.0101040127121235, 0.0097278039838056,      0.00938617839440894, 0.00908600262449913, 0.00883795563273938,      0.00863315495962393, 0.00847325872834633, 0.00835985054125376,      0.0083104707966061, 0.00831574579983596, 0.00837242006745159,      0.0084811812661746, 0.00864895114585229, 0.0088847277032945,      0.0091733846139949, 0.00951451815849311, 0.00990760864554539,      0.0103676944983387, 0.0108803185070754, 0.0114398691731338, 0.0120446169019066,      0.0126975354786556, 0.0134003348410125, 0.014138322350119, 0.0149085888122459,      0.0157081161739461, 0.0165397010226779, 0.0173908415086893, 0.0182558660392905,      0.0191308446396648, 0.0200115146141704, 0.0208896684310384, 0.0217595621107009,      0.0226168398271489, 0.0234571367239625, 0.0242658576551506, 0.0250415384796808,      0.025781743996249, 0.0264824465880188, 0.0271330271813973, 0.0277197602228349,      0.0282514207578494, 0.0287250463560853, 0.0291378145433322, 0.0294643336036007,      0.0297185497373367, 0.0299034233614394, 0.0300178272409036, 0.0300510165057044,      0.0299933390096322, 0.029864296171915, 0.0296646737354628, 0.0293954446422498,      0.02903609387537, 0.0286087954739586, 0.0281202389214091, 0.0275730852809592,      0.0269631864592354, 0.0262908755920268, 0.025574205859082, 0.0248169722677519,      0.024023046736675, 0.0231883685262104, 0.0223289370148913, 0.0214507109581384,      0.0205578894465608, 0.0196543161993343, 0.0187475783069642, 0.0178440120329414,      0.0169474708137715, 0.0160617456778675, 0.0151982413151739, 0.0143579764501946,      0.0135431035900457, 0.0127565541603801, 0.0120065112945572, 0.0113013352426194,      0.010635019980095, 0.0100094469431421, 0.00942639127892244, 0.0089039890191574,      0.00843070966487366, 0.00800519340240519, 0.00762818830262024,      0.00730797543721608, 0.00704980283664362, 0.00684141502281347,      0.00668265968608497, 0.00657330236985191, 0.00652994734551145,      0.00653626373538692, 0.00658916348197513, 0.00668767640137901,      0.00683740677506747, 0.00703949814277543, 0.00728169259995597,      0.00756244741730784, 0.00788017305139687, 0.00824447686134098,      0.00864128336177512, 0.00906733428314206, 0.00952075397149063,      0.0100031836945359, 0.010512665615019, 0.0110414546541446, 0.0115878311927959,      0.0121501195536992, 0.0127308183673382, 0.01332332344487, 0.0139262209784856,      0.0145387177488175, 0.015161575268352, 0.0157951749967284, 0.016438040193859,      0.0170909905985346, 0.0177550383014629, 0.0184378178562425, 0.019138615981402,      0.0198599761627206, 0.0206051085684539, 0.0213835326226107, 0.0222057991861808,      0.0230705544613072, 0.0239831651499172, 0.0249491866774483, 0.0260019858815712,      0.0271307566962957, 0.0283404432091247, 0.0296377725554787, 0.0310483777483443,      0.0325909750094088, 0.0342501536507845, 0.0360320721161026, 0.0379427304435736,      0.0400427505040893, 0.0422911469270196, 0.0446880555428848, 0.0472364098981542,      0.0499665903504324, 0.0528865496293085, 0.0559605500295509, 0.0591864336230866,      0.0625614652392627, 0.0661323901936574, 0.0698369942423732, 0.0736642649820459,      0.0776053170213216, 0.0816654985557755, 0.0858243890440125, 0.0900462454820139,      0.0943168565446659, 0.0986215759610166, 0.102937806780721, 0.107234576343045,      0.111494854705369, 0.115701725627134, 0.119818728151249, 0.123810263483841,      0.127676838260079, 0.131403340542788, 0.134975040188428, 0.138301149786118,      0.141429558148372, 0.144353043040705, 0.147063324586097, 0.149509011012514,      0.151678338013187, 0.153617609445566, 0.155327745461799, 0.156810642833852,      0.157987822962458, 0.15895731472489, 0.159731001895812, 0.160320192816835,      0.16070927729284, 0.160927335641176, 0.161026262422779, 0.161024035352265,      0.160939101818898, 0.160783676919756, 0.160610951636628, 0.160440813332408,      0.160292840283904, 0.160202196337883, 0.160205045445256, 0.160307691092289,      0.160525158759831, 0.160872099587403, 0.161427322298974, 0.162146341126922,      0.16303494787568, 0.1640978213177, 0.165374536988198, 0.166861439709314,      0.168517967752217, 0.170337253164731, 0.172312150290867, 0.174476220294316,      0.17675201562322, 0.17912237546253, 0.181569351447975, 0.184077063209181,      0.186601937519015, 0.189110996532675, 0.19158078566652, 0.193986444864657,      0.196248038237016, 0.198366671014338, 0.2003193238609, 0.202083414374645,      0.203583109911179, 0.204784727823111, 0.20571235693773, 0.206350205706743,      0.206678507619118, 0.206556713002678, 0.206095969055307, 0.205290936474901,      0.20413741111321, 0.202556572638439, 0.20056355298721, 0.198228613635869,      0.195558644146613, 0.192555647948946, 0.189126333152985, 0.185410492636506,      0.181423717899618, 0.177182210388674, 0.172658716028772, 0.167900183255391,      0.16296802586402, 0.157882368007903, 0.152661212543965, 0.147304635977954,      0.141882524743064, 0.136414358489409, 0.130919304274034, 0.125422474679471,      0.119957697079395, 0.11454052598248, 0.109185666974252, 0.103909993572189,      0.0987683404859644, 0.093741308688079, 0.0888376270971578, 0.0840654236882268,      0.0794666168513133, 0.0750427325158831, 0.0707727313029462, 0.0666591512693038,      0.0627091344513033, 0.0589774733641864, 0.0554044324211085, 0.0519885093076509,      0.0487279204345587, 0.0456550213639223, 0.0427535393057037, 0.0399945562439461,      0.0373743649778989, 0.0348937863392637, 0.0325878723180957, 0.0304033659564744,      0.0283361981437062, 0.0263823330082348, 0.0245619867105529, 0.0228595658836749,      0.0212541351296348, 0.0197423154392725, 0.0183242669772017, 0.0170220160814398,      0.0158005734119126, 0.0146575124941681, 0.0135905002875402, 0.0126143097870055,      0.0117191044664318, 0.0108918081116421, 0.0101309359870987, 0.0094379000631823,      0.00883019944826195, 0.00828380391399826, 0.00779783511360308,      0.00737145971429946, 0.00701812393006059, 0.0067311252542151,      0.00650097052493056, 0.00632717374230014, 0.00621198669319391,      0.00617143557767508, 0.0061855594050197, 0.00625405591658972,      0.00637663424018689, 0.00656643085439909, 0.00681747969859399,      0.00712155184640286, 0.00747843433306984, 0.0078907086989986,      0.00837340291572697, 0.00890806427672648, 0.009494513737272,      0.010132574128414, 0.010835143782014, 0.0115961949265227, 0.0124081127363209,      0.0132707231223347, 0.0141867435373067, 0.0151701717571193, 0.0162034512423126,      0.0172863373562776, 0.0184185690183941, 0.019612438019485, 0.0208613703536178,      0.0221579723104729, 0.0235016744919049, 0.0248946535533411, 0.0263479485897112,      0.0278446638511839, 0.0293835133691365, 0.0309631084740302, 0.0325912348690895,      0.0342591472113346, 0.0359590659180276, 0.0376882779600787, 0.0394452615206507,      0.0412293518464991, 0.043027543361404, 0.0448353514472685, 0.0466481043572906,      0.0484577337141087, 0.0502538712990142, 0.0520305447610213, 0.0537813149370389,      0.0554964047538021, 0.0571515262745735, 0.0587502027060215, 0.0602849490847391,      0.0617482500354599, 0.0631066082295335, 0.0643579690671984, 0.0655057363703737,      0.0665429490020974, 0.0674531640036809, 0.0681929810868339, 0.0687958988058021,      0.0692572514303207, 0.0695726946039624, 0.069695167044965, 0.0696433847901905,      0.0694350823030356, 0.0690697581056959, 0.068535116834938, 0.0677961078509799,      0.0669066143489608, 0.0658704054221769, 0.0646916508125923, 0.0633398550663818,      0.0618485040466386, 0.0602410483121883, 0.0585250659382392, 0.0567013372657685,      0.0547654519871085, 0.0527582945347469, 0.0506891518402187, 0.0485673654999008,      0.0463957642449994, 0.0441983057887947, 0.0419884488508511, 0.0397750260303709,      0.0375684836047041, 0.0353879440268247, 0.0332384672584751, 0.0311267921300081,      0.0290593901126105, 0.0270601872571371, 0.0251287992452735, 0.0232632081303122,      0.0214669543868477, 0.0197503127220639, 0.0181339971366657, 0.0165964380609598,      0.0151383742283497, 0.013760299683845, 0.012485591080821, 0.0112979221667829,      0.0101874296294859, 0.00915261833812058, 0.00819857079937494,      0.00733611052320258, 0.00654025398573884, 0.00580844808959148,      0.0051380838601794, 0.00454221827421074, 0.00400446040436533,      0.00351615254678973, 0.0030745243393893, 0.00268072850148871,      0.00233768602660161, 0.00202984861125548, 0.00175479227131742,      0.00151014383233181, 0.00130072566668639, 0.00111721604283242,      0.000955112980316734, 0.000812623159815952, 0.000689590092917738,      0.0005860912862362, 0.000495639552545513, 0.000417014922709681,      0.000349049284657548, 0.000292958836410902, 0.00024515801710463,      0.000204040052154968, 0.000168884867932852, 0.000139474416390484,      0.000115548025328669, 9.51646202245106e-05, 7.79159203801693e-05,      6.34177645960676e-05, 5.18748729574288e-05, 4.22967579760247e-05,      3.42684953287067e-05, 2.75891372479148e-05, 2.21737075433437e-05,      1.79085413159829e-05, 1.43641470627961e-05, 1.14430817192601e-05,      9.05508126397165e-06, 7.22026287729582e-06, 5.73639534240461e-06,      4.52392445545382e-06, 3.54207412160509e-06, 2.77052765004019e-06,      2.18172878629954e-06, 1.7042057736597e-06, 1.32082991902191e-06,      1.01596818604684e-06, 7.89867561242482e-07, 6.11506981747829e-07,      4.69391759719292e-07, 3.57363517027134e-07), bw = 0.715374250961732,          n = 35L, call = density.default(x = dat.prior, from = 0,              to = 20), data.name = "dat.prior", has.na = FALSE), .Names = c("x",      "y", "bw", "n", "call", "data.name", "has.na"), class = "density")


title: Relation of Mahalanobis Distance to Log Likelihood

text: The Wikipedia entry on Mahalanobis Distance contains this note: > Another intuitive description of Mahalanobis distance is that it is square > root of the negative log likelihood. That is, the exponential of the > negative square of the Mahalanobis distance will give you the likelihood of > your data point belonging to (a presumed normal) distribution of the sample > points you already have. How/where is this shown? What does it mean (I guess it's not so intuitive...)?


title: Exponential Distribution Expected Value

text: **Question 1 (part 1)** I was attempting to use the exponential distribution to solve the following problem: > the time between cars passing on a road can be modelled with an avg rate of > 3.5 cars per minute. Find the probability that the time between two > comsecutive cars passing the road is between 20 and 40 seconds I have set lamda to be 0.0583/sec. **Question 1 (part 2)** **While what I am going to ask has nothing to do with solving my question, it leaves me curious (since i am really sincere about learning stats):** If I did E(x) = 1/lamda = 17.15 May i say that **this value is the average cars passing per sec**? But it doesnt seem to make sense here * * * **Answer to Question 1 (part 1)** The answer to the original question is: Find P (20 < x < 40) Lamda = 0.0583 per sec P(x<40) = 1 - e ^ -(0.0583)(40) = 0.9029 P(x<20) = 1 - e ^ -(0.0583)(40) = 0.6884 Therefore P (20 < x < 40) = 0.9029 - 0.6884 = 0.2145 If theres a better way to do this, I will be glad to learn.


title: Empty significance in SPSS linear regression

text: I got a weird result with linear regression: there are empty significances. I use 'Stepwise' method, so it should drop the inadequate variables. (Of course the "000" Betas are confusing too.) What could be the reason, and what should be the solution? Thank you very much in advance. ![enter image description here](http://i.stack.imgur.com/Q3wJg.jpg)


title: Statistical tests to check quality of data before learning model

text: I have collected some gene expression data and I want to learn a Bayesian network out of it. Before that, I want to do some statistical analysis to test the quality of my data. Now I want to know that, which kind of statistical analysis I should do to assess the quality of my data?


title: Type I ANOVA tests not depending on the order of the factors

text: I have a dataset with two factors `A` and `B` and the following design (contigency table showing the number of individuals for each crossed level):               > xtabs(~A+B, data=dd2)         B     A    B1 B2 B3 B4       A1  8  8  8  8       A2  7  7  7  7       A3  8  8  8  8       A4  4  4  4  4       A5  4  4  4  4      I have observed that the Type I (sequential) ANOVA sum of squares do not depend on the order of `A` and `B` in the model. I was firstly surprised because the design is unbalanced and I have always believed that balanced design is the general condition required in order to have this property. Here the design is only balanced per each level of factor `A`. So finally what is the general condition under which the Type I sum of squares do not depend on the order of the factors ?


title: Pattern of events that leads to a goal

text: I am not a statistician, so go easy on me here... I have collected lots of data on the behaviour of customers in a shop. A customer will enter the shop and a series of events will then take place (eg look at product, look at advert in shop, talk to shop assistant, etc). Before they leave the shop, some customers will make a purchase. If my goal is to be able predict which future customers will make a purchase, what statistical analysis should I apply to the data? I basically want to see if there is a pattern of events that, when followed, will indicate that a customer is more likely to make a purchase. For example, it may be that customers who look at a particular advert in the shop, then talk to the shop assistant, then touch a product are very likely to go on and purchase.


title: Computing per cent correctness from a bell curve

text: I've been trying to remember my High School teachings and are falling short. I'm working on a project where I need to give a % of correctness for an integer (How close a given number is to the actual number within a 300% difference). For example if the number we want is 50, any number from -100 to 150 will return > 0% correctness. The problem is that we need a curve (log, or bell curve, or something similar) to return a non-linear % correctness (i.e. 100 is 50% correct in linear terms, but we would want maybe 66%? ... and 125 is 33%?) and I dont have a formula to get this response. Something like this (sorry, used mspaint quickly to try to explain) (please go to http://stackoverflow.com/questions/7759062/math-statistics-bell-curve- computing-correct-given-2-numbers-c to see the image, cant post it here due to me being a new user) Am I explaining this properly? Make sense? Any help? ;) I ran into standard deviation too, just its a bit complicated for me to process right now. If you understand it, can you throw me a quick formula?


title: Text classification with RapidMiner

text: I have a folder for each class that contains some text files related to that class and need to implement a classifier for these files in RapidMiner with KNN algorithm. Is there any sample?


title: Equivalence of Poisson and Weibull PH regression in a survival setting

text: Suppose we have the following dataset that records individual survival times (dur) and a covariate z:               id| dur | z     -------------     1 |  1  | -1     2 |  2  |  1     3 |  3  | -1     4 |  4  |  1      I want to model the duration as a function of z. I may specify dur ~ weibull() and parameterize the scale as a function of z. This model can be fitted easily, e.g. with phreg or survreg in R. If I am now interested to incorporate time-dependent covariates, I need to transform the dataset into something like this:               id | event  | time  | z  | orig.row     _____________________________________     1  | 1      | 1     | -1 | 1     2  | 0      | 1     | 1  | 2     2.1| 1      | 2     | 1  | 2     3  | 0      | 1     |-1  | 3     3.1| 0      | 2     |-1  | 3     3.2| 1      | 3     |-1  | 3     4  | 0      | 1     | 1  | 4     4.1| 0      | 2     | 1  | 4     4.2| 0      | 3     | 1  | 4     4.3| 1      | 4     | 1  | 4      Apparently, a to say event ~ poisson() with a "poisson mean" $$\log(u_i) = \beta_0 + \alpha \log(t_i) + \beta_1 z_i$$ is equivalent to a the Weibull model above (Lindsey 1995). But when I run the analysis in R, I get two largely different values for the scale:               # Generate the data:      library(eha)     enter <- rep(0, 4)     exit <- 1:4     event <- rep(1, 4)     z <- rep(c(-1, 1), 2)     dat <- data.frame(enter, exit, event, z)     binDat <- toBinary(dat)     binDat <- binDat[order(rownames(binDat)),]          # Run the model:      summary(phreg(Surv(enter,exit, event) ~ z, data = dat, dist="weibull"))     summary(glm(event ~ z + log(risktime), data = binDat, family = poisson("log")))      The shape parmeter (glm: log(risktime)) is sort of similar as well as the beta for z. But the scale parameter is different (glm: Intercept). What am I doing wrong?               ## Results     summary(phreg(Surv(enter,exit, event) ~ z, data = dat, dist="weibull"))     Call:     phreg(formula = Surv(enter, exit, event) ~ z, data = dat, dist = "weibull")          Covariate          W.mean      Coef Exp(Coef)  se(Coef)    Wald p     z                   0.200    -0.432     0.649     0.519     0.406           log(scale)                    1.020     2.774     0.197     0.000      log(shape)                    0.985     2.677     0.421     0.019           Events                    4      Total time at risk            10      Max. log. likelihood      -5.6567      LR test statistic         0.68      Degrees of freedom        1      Overall p-value           0.410367          summary(glm(event ~ z + log(risktime), data = binDat, family = poisson("log")))          Call:     glm(formula = event ~ z + log(risktime), family = poisson("log"),          data = binDat)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -1.0696  -0.7695  -0.5391   0.3482   1.0671            Coefficients:                   Estimate Std. Error z value Pr(>|z|)     (Intercept)    -1.6122     1.0019  -1.609    0.108     z              -0.3166     0.5149  -0.615    0.539     log(risktime)   1.0634     1.0928   0.973    0.330          (Dispersion parameter for poisson family taken to be 1)              Null deviance: 7.3303  on 9  degrees of freedom     Residual deviance: 6.1388  on 7  degrees of freedom     AIC: 20.139          Number of Fisher Scoring iterations: 5


title: Coverage probability

text: I know the definition, but what does it imply if one estimator has greater coverage probability than another estimator? Does it imply that it has converges to the parameter of interest faster?


title: How to keep time invariant variables in a fixed effects model

text: I have data on a large Italian firm's employees over ten years and I would like to see how the gender gap in male-female earnings has changed over time. For this purpose I run pooled OLS: $$ y_{it} = X'_{it}\beta + \delta {\rm male}_i + \sum^{10}_{t=1}\gamma_t d_t + \varepsilon_{it} $$ where $y$ is log earnings per year, $X_{it}$ includes covariates that differ by individual and time, $d_t$ are year dummies and ${\rm male}_i$ equals one if a worker is male and is zero otherwise. Now I have a concern that some of the covariates can be maybe correlated with unobserved fixed effects. But when I use the fixed effects (within) estimator or first differences I lose the gender dummy because this variable does not change over time. I don't want to use the random effects estimator because I often hear people saying that it puts assumptions that are very unrealistic and are unlikely to hold. Are there any ways for keep the gender dummy and control fixed effects at the same time? If there is a way, do I need to cluster or take care for other problems with the errors for hypothesis tests on the gender variable?


title: LMM models with two difference reference categories yield two different results

text: I have run LMM models with different reference categories and this yield different results:               > summary(lmer3)     Linear mixed model fit by maximum likelihood ['merModLmerTest']     Formula: v000001 ~ (1 | item) + (1 + color | speaker) + Language * color *      sex         Data: data1.frame                 AIC       BIC    logLik  deviance      16279.975 16377.355 -8119.988 16239.975           Random effects:      Groups   Name        Variance  Std.Dev.  Corr             speaker  (Intercept) 8.904e+05 9.436e+02                           colorblue   1.821e+05 4.267e+02 -0.35                     colorred    3.428e+05 5.855e+02 -0.44  1.00      item     (Intercept) 9.502e-06 3.083e-03                  Residual             1.067e+06 1.033e+03                 Number of obs: 962, groups: speaker, 53; item, 10          Fixed effects:                                       Estimate Std. Error       df t value Pr(>|t|)     (Intercept)                       10664.67     318.69    38.45  33.464   <2e-16     Languagel2_like                     391.48     421.40    42.13   0.917   0.3642     colorblue                          -179.31     211.02    44.50  -0.850   0.4000     colorred                            116.96     241.44    36.27   0.484   0.6310     sexmale                            -168.01     450.11    38.26  -0.373   0.7110     Languagel2_like:colorblue           758.22     301.01    54.20   2.519   0.0147     Languagel2_like:colorred            463.37     344.01    45.73   1.344   0.1857     Languagel2_like:sexmale            -811.49     607.85    43.49  -1.326   0.1917     colorblue:sexmale                   342.76     294.97    42.57   1.162   0.2517     colorred:sexmale                     13.25     337.44    34.81   0.039   0.9689     Languagel2_like:colorblue:sexmale  -721.37     438.78    54.19  -1.644   0.1059     Languagel2_like:colorred:sexmale   -605.76     497.75    45.29  -1.216   0.2304          (Intercept)                       ***     Languagel2_like                           colorblue                                 colorred                                 sexmale                                   Languagel2_like:colorblue         *       Languagel2_like:colorred                 Languagel2_like:sexmale                   colorblue:sexmale                         colorred:sexmale                          Languagel2_like:colorblue:sexmale         Languagel2_like:colorred:sexmale          ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      And this is the second model:               > summary(lmer43)     Linear mixed model fit by maximum likelihood ['merModLmerTest']     Formula: v000001 ~ (1 | item) + (1 + color3 | speaker) + Language * color3 *      sex         Data: data1.frame                 AIC       BIC    logLik  deviance      16279.975 16377.355 -8119.988 16239.975           Random effects:      Groups   Name        Variance  Std.Dev.  Corr             speaker  (Intercept) 7.945e+05 8.913e+02                           color3white 1.821e+05 4.268e+02 -0.11                     color3red   2.761e+04 1.661e+02 -0.24 -0.94      item     (Intercept) 4.961e-06 2.227e-03                  Residual             1.067e+06 1.033e+03                 Number of obs: 962, groups: speaker, 53; item, 10          Fixed effects:                                        Estimate Std. Error       df t value     (Intercept)                        10485.36     305.33    39.57  34.341     Languagel2_like                     1149.70     399.61    43.91   2.871     color3white                          179.31     211.03    44.50   0.850     color3red                            296.27     167.08   125.59   1.773     sexmale                              174.75     430.05    38.94   0.406     Languagel2_like:color3white         -758.22     301.01    54.10  -2.519     Languagel2_like:color3red           -294.85     244.46   159.09  -1.206     Languagel2_like:sexmale            -1532.85     577.70    44.74  -2.648     color3white:sexmale                 -342.76     294.98    42.57  -1.162     color3red:sexmale                   -329.51     228.57   113.99  -1.442     Languagel2_like:color3white:sexmale  721.36     438.78    54.10   1.644     Languagel2_like:color3red:sexmale    115.61     351.65   162.98   0.329                                        Pr(>|t|)         (Intercept)                         < 2e-16 ***     Languagel2_like                     0.00627 **      color3white                         0.40004         color3red                           0.07862 .       sexmale                             0.68671         Languagel2_like:color3white         0.01477 *       Languagel2_like:color3red           0.22953         Languagel2_like:sexmale             0.01114 *       color3white:sexmale                 0.25171         color3red:sexmale                   0.15215         Languagel2_like:color3white:sexmale 0.10602         Languagel2_like:color3red:sexmale   0.74275         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      Is it possible to report results from two models? (I know that very little of people would do this but these two models give different pictures). What should I do? Which one should I believe?


title: Automatic generation of crosswords

text: When solving crosswords today, I thought of if they can be generated automatically by any mathematical approach. My question is: If database of words and any grid with some cells already filled (solution of crossword) are given, how can ML (or maybe any optimization?) be applied to generate automatically at least one crossword?


title: Is repeated measures ANOVA the right test?

text: I am studying a parameter in several organs of mice and I have two categorial factors. I would like to study the effects of those two factors on my parameter in each organ (I don't care about differences of the parameter bewteen organs). I was previously doing a two-way ANOVA for each organ but this led to multiple ANOVA tests. I think this could be a problem since it increases the risk of having false positive results. To adress the problem of multiple comparisons, could I perform a single repeated measure ANOVA test, taking organs as the within-subject factor (even though I don't want to study this factor) and the two remaining factors as between-subject factors ? In my mind, repeated measures ANOVA is mainly used to measure a parameter at different time points so I'm not quite sure about this test in this case. Thank you !


title: What are the assumptions of negative binomial regression?

text: I'm working with a large data set (confidential, so I can't share too much), and came to the conclusion a negative binomial regression would be necessary. I've never done a glm regression before, and I can't find any clear information about what the assumptions are. Are they the same for MLR? Can I transform the variables the same way (I've already discovered transforming the dependent variable is a bad call since it needs to be a natural number)? I already determined that the negative binomial distribution would help with the over-dispersion in my data (variance is around 2000, the mean is 48). Thanks for the help!!


title: Using R online - without installing it

text: Is there a possibility to use R in a webinterface without the need to install it? I have only one small script which I like to run but I just want to give it a shot without a long installation procedure. Thank you.


title: Linear trend in SAS using contrast

text: I've been struggling with this problem for a couple of hours, and I could use some advice. I have a linear regression model with two continuous predictors and a categorical one (with 4 levels). I want to test for a linear trend in the categorical variable, with the lowest level as the reference. I've gotten the answer in Stata and R, but I can ONLY get the answer in SAS if I use its default reference level (which would be 4). If I change the reference level to 1, I get a different answer. Stata (F=18.56):               contrast {Age -3 -1 1 3}      R (F=18.56):               fit.contrast(model_A, "age_cat", c(-3, -1, 1, 3))      SAS (F=0.06):               proc glm data=data;     class age (ref='1');     model outcome=continuous1 continuous2 age / solution ss3;     contrast 'linear trend' age -3 -1 1 3;     run;      I figure that since I'm adjusting the reference category I have to adjust something else, but I cannot for the life of me figure out what it is. Thanks in advance for your help.


title: Obtaining an estimator via Rao-Blackwell theorem

text: Let $X_1,...,X_n$ be iid with pdf $$f(x|\theta) = exp(\theta -x) I(x)_{(\theta, \infty)}$$ It is asked to find an unbiased estimator for$ \theta $ that is function of a sufficient statistical for $\theta$ By factorization theorem, we show that $X_{(1)}$ is a sufficient statistical for $\theta.$ And, since $$E(X) = \theta +1$$ the estimator $ \bar{X} -1$ is unbiased. So, by the Rao Blackell theorem, $$W=E(\bar{X}-1|X_{(1)})$$ is an unbiased estimator that is function of the sufficient statistical. But it seems very complicated to evaluate the distribution of $\bar{X}-1|X_{(1)}$. How can I find this distribution? Is there a better unbiased estimador that I can use, in this case? Thanks in advance!


title: What kind of model can I try to fit in this plot?

text: I have a plot like this. ![enter image description here](http://i.stack.imgur.com/bru73.png) I wish to apply a model to this, however, I guess a linear regression model won't work on this. What I did was plot it on logarithm x and logarithm y axis as well but it came out to be of no use. With logarithm: ![https://www.dropbox.com/s/vyzk4u26vnludoq/1.png?dl=0](http://i.stack.imgur.com/blNav.png) I tried fitting in a model, but as expected, on plotting, the residual and fitted, it didn't turn out to be of much use: ![https://www.dropbox.com/s/0uaf53bxdc3wd1l/3.png?dl=0](http://i.stack.imgur.com/Dy7Zs.png) What else can I do? Anything that someone can suggest? Also, I wish to know how can I apply a non linear regression model?


title: Difficulty of testing linearity in regression

text: In Statistical Modeling: The Two Cultures Leo Breiman writes > Current applied practice is to check the data model fit using goodness-of- > fit tests and residual analysis. At one point, some years ago, I set up a > simulated regression problem in seven dimensions with a controlled amount of > nonlinearity. Standard tests of goodness-of-fit did not reject linearity > until the nonlinearity was extreme. Breiman doesn't give the details of his simulation. He references a paper that he says gives theoretical justification for his observation, but the paper is unpublished. Has anyone seen a published simulation result or theoretical paper to support Brieman's claim?


title: Prediction in simple and multiple ANOVA

text: In a (one or multi) way anova model, once a new individual is assigned to a treatment, the predicted value for him is calculated using the coefficients of the ANOVA model (simply assigning the treatment mean value to the individual). How should I construct a confidence (or prediction) interval for that predicted value? Can the `predict` function be used in R for an `aov` model? There must be something that I'm missing here, as I can't find anything like this in my books. Thank you.


title: Hot topics in mathematical statistics

text: What are some hot topics that mathematical statistics researchers are studying now?


title: Equivalent of Tukey test for chi-square?

text: As a follow up to my previous question, I learned that the chi square test can be used to compare expected and observed frequency distributions. But... Suppose I have the following data table in CSV, which surveys people on their fruit preferences in three cities (total people surveyed is different in each city):               apples,oranges,grapes,bananas          Acme Acres,5,15,5,15          Redmond,2,3,16,0          Cupertino,9,13,8,4      (1) I can run a chi-square test on each city, to see if my observed frequencies are significantly different from expected (in this case, equally distributed). But is there a test that I can run on _ALL THREE_ cities at once to see if they are different **from each other** , _AND_ also if they are individually different from the expected equal distribution? (2) If I can run a test that says there is a difference between cities, and/or if each one is different from the expected distribution, is there an equivalent of the Tukey test where I can see which specific fruits are the ones that are different? Thanks!


title: Is the case when response is discrete and ordinal classification?

text: When predicting the value of a response variable from the value of an explanatory variable, if I am correct,   * when the response is real-valued, the task is regression, and    * when the response takes values from a unordered discrete set, the task is classification. When the response takes its value from an ordered discrete set, is the task still classification? For example, let $Y$ be the response variable and $X$ be the explanatory variable. If the conditional distribution of $Y$ given $X=x$ is a Poisson distribution, is the task of predicting the value of $Y$ from the value of $X$ still classification? Is Poisson regression used for classification, just as logistic regression for classification? Thanks!


title: Distribution of the Rayleigh quotient

text: For a research project I need to find the expected value of the generalized Rayleigh quotient: $$E\,[w^T A w \ / \ w^T B w].$$ Here _A_ and _B_ are positive definite deterministic _p_ x _p_ covariance matrices, and _w_ follows a multivariate distribution with circular altitude lines (say, multivariate standard normal). The dimension _p_ is greater than 100. This problem is easy to solve using simulation; however, I was wondering if somebody might know how this problem could be solved (or approximated) analytically. My first idea was that possibly by the Lindeberg or Lyapunov central limit theorem both the numerator and the denominator are approximately normal distributed, which gives us a ratio of two (correlated) normal random variables, but simulation shows that that is not the case.


title: How does gbm determine the missing node if the training set is complete?

text: I have a training set which is complete - no missing values, but yet if I pretty print one of the trees in gbm I see that the algorithm has determined a strategy for missing data. This is smarter than I thought it would be, because in the future the data may not always be complete. I'm wondering if the tree just takes the mean of the target variable for the training observations that end up in the missing node The below output is from the R implementation of gbm. You'll see that the first split is done on Variable 1, and if it is missing it goes to row 9 which describes a terminal node with prediction -0.018                 SplitVar SplitCodePred LeftNode RightNode MissingNode ErrorReduction Weight  Prediction     0        1    0.00000000        1         8           9      4815.9182 199782 -0.01823429     1      117    1.00000000        2         3           7      1846.5907 184985 -0.05852970     2       -1   -0.12162550       -1        -1          -1         0.0000  69534 -0.12162550     3      115    2.00000000        4         5           6       402.3653 115451 -0.02052827     4       -1   -0.08380513       -1        -1          -1         0.0000  28280 -0.08380513     5       -1    0.00000000       -1        -1          -1         0.0000  87171  0.00000000     6       -1   -0.02052827       -1        -1          -1         0.0000 115451 -0.02052827     7       -1   -0.05852970       -1        -1          -1         0.0000 184985 -0.05852970     8       -1    0.48551969       -1        -1          -1         0.0000  14797  0.48551969     9       -1   -0.01823429       -1        -1          -1         0.0000 199782 -0.01823429


title: Label Switching in WinBugs/JAGS

text: I am using JAGS to estimate a Dirichlet Process Mixture of Normals. The code works well and the estimated density is accurate. However, I would like to know which component each observation is assigned to and the corresponding parameters for that component. This is hard due to the label switching problem in mixture models. Any suggestions? Putting an order restriction on the means is one recommended way to avoid this but I have not been able to implement it well in JAGS. Below is the R code followed by the JAGS model.               library(R2jags)     library(mixtools)          nobs=200     popmn=c(-6,-3,3); popstd=c(1.5,1.5,.7); popprob=c(.6,.3,.1)     x=rnormmix(n=nobs,lambda=popprob,mu=popmn,sigma=popstd)     xgrid=seq(-8,8,length=500)          ###true density     densTrue=rep(0,length(xgrid))     for(i in 1:length(xgrid))       densTrue[i]=sum(popprob*dnorm(xgrid[i],popmn,popstd))     densTrue=densTrue/sum(densTrue)          ###Call JAGS     K=10     dataList = list(y=x,n=200,K=K,grid=xgrid,n2=length(xgrid))     jags.inits <- function(){       list('mu'=rep(0,K),'tau'=rep(1,K),'prob'=rep(1/K,K))     }     parameters=c('f')     jags.fit=jags(data=dataList,inits=jags.inits,parameters,"jagsModel",n.chains=1,n.iter=800,n.burnin=300)     jags.fit     #contruct density     DPDens2=jags.fit$BUGSoutput$mean$f     DPDens2=DPDens2/sum(DPDens2)          ###Plot     plot(xgrid,trueDens,type='l')     lines(xgrid,DPDens2,col=2)      JAGS Model:               model{      for(i in 1:n){       y[i] ~ dnorm(mu[z[i]],tau[z[i]])       z[i] ~ dcat(prob[])       for(j in 1:K){        f.d[i,j] <- prob[j]*coi*exp(-0.5*pow(y[i]-mu[z[i]],2)*tau[z[i]])*sqrt(tau[z[i]])       }      }      prob[1:K] ~ ddirch(prob.par[])      for(i in 1:K){       repeat[i] <- step(1-dist[i])       prob.par[i] <- alpha/K       mu[i] ~ dnorm(0,.001)       tau[i] ~ dgamma(0.001,0.001)      }      # total is number of active "used" components      total <- sum(repeat[])      dist[1] <- 1      for(i in 2:K){       dist[i] <- sum(count[i,1:i])       count[i,i] <- 1       for(j in 1:i-1){        count[i,j] <- equals(z[i],z[j])       }      }      alpha <- 1     ###Density      for(i in 1:n2){       for(j in 1:K){        f.p[i,j] <- prob[j]*coi*exp(-0.5*pow(grid[i]-mu[j],2)*tau[j])*sqrt(tau[j])       }       # f[i] is predictive density over grid[i]       f[i] <- sum(f.p[i,1:K])       }       coi <- 0.3989422804     }


title: What are some basic ways to build a binary classifier based on questionaire data?

text: I apologize if this is a basic question - I'm a complete beginner with machine learning & stats. I am working on a website where signing up for a premium feature requires the user to fill out a questionnaire prior to adding the premium. So I have a dataset of:   1. Users who completed the questionnaire but decided not to sign up   2. Users who completed the questionnaire but did not sign up The questionnaire consists of multiple choice answers (i.e. categorical variable) and numerical inputs. You can think of the experience as similar to the questionnaire on Wealthfront (https://www.wealthfront.com/questions). The survey consists of 10 questions, 9 of which are multiple choice (i.e. categorical) & 1 which is a number. I have approx 10K filled out surveys and ~300 sign-ups (3% sign-up rate). What I need to find out is:   * Given a set of answers to the questionnaire, what is the probability that a user will sign up?   * What question/answer combinations are most likely to result in a sign up? I don't have any idea how to get started. What are some basic approaches to do this?


title: Using z scores in paired t-tests and regression analysis

text: I have one group of participants measured on a health outcome scale at two different time points. I first want to see whether there is a significant difference between time 1 and time 2 on this outcome variable, I then want to see if this outcome variable at T2 is predicted by another dimensional variable (BPD dimensional score). Seems straightforward using a paired samples t-test and then a regression analysis. However the health outcome variable has an adolescent version and a adult version. At time 1 most of the participants used the HONOSCA (adolescent version) and some used the HONOS (adult version). at the second time point the majority used the HONOS with only a few using the HONOSCA. while there have been numerous studies evidencing that the adult and adolescent versions can be compared, I have not found anything with the use of the different versions used for the same outcome variable. One scale has 13 items and the other 12 items, so total scores cannot be used. I have transformed the variables into z scores to make them comparable, but is it valid to use the z scores in further analysis? How can I compare T1 and T2 outcomes because the z-scores are not valid for the t-tests? Thanks


title: Mathematically modeling neural networks as graphical models

text: I am struggling to make the mathematical connection between a neural network and a graphical model. In graphical models the idea is simple: the probability distribution factorizes according to the cliques in the graph, with the potentials usually being of the exponential family. Is there an equivalent reasoning for a neural network? Can one express the probability distribution over the units (variables) in a Restricted Boltzmann machine or a CNN as a function of their energy, or the product of the energies between units? Also, is the probability distribution modelled by an RBM or Deep belief network (e.g. with CNNs) of the exponential family? I am hoping to find a text that formalizes the connection between these modern types of neural networks and statistics in the same way that Jordan & Wainwright did for graphical models with their Graphical Models, Exponential Families and Variational Inference. Any pointers would be great.


title: Bias term in support vector machine

text: In SVM, there is a bias term. But looks to me there are very few discussions on the physical meanings of this term. Why should we have that? How does this term affect the model?


title: Fit a smooth approximation line

text: I have a simulated data set and I want to fit a "smooth approximation line" like the image I have provided. ![http://i1024.photobucket.com/albums/y307/sarmadys/fitting_zps9d9651e7.jpg \(I did not have enough rep to post image\)](http://i.stack.imgur.com/aPE28.jpg) 1- Is it possible to do this in Excel or Matlab? May I have a pointer on how to do it? 2- How we find how much of the data is within a distance of 0.15 of the line? Or do we basically fit a line in a way that 90% of the points fall in 0.15 range?


title: Help modeling school consent rate changes after an intervention

text: I have # of students consenting to be vaccinated and # of students eligible to be vaccinated for 77 schools in three school years (2010/11, 2011/12, and 2012/13). Between 2011/12 and 2012/13 there was an intervention to increase consent rates. I want to be able to model (with a repeated measures model) the change in consent rates over time. First, I am having trouble deciding what distribution the outcome is since it is at the school level and not the individual-level. It's not a simple 1-consent, 0 - no consent as schools contribute multiple counts of #consents and #eligible. It's a rate but does not really follow a poisson distribution. Question 1 - is an aggregrate rate across schools, say in 2010/11, just the mean of individual school rates (to me this doesn't make sense) or the sum of #consents/sum of #eligible across all schools? If it's the latter, do I use a binomial or poisson distribution to calculate CIs, or does it matter? Question 2 - Since I don't have student level data, does it make sense to use GEE with the poisson family and log-link? This would be using count data of #consents per school and offset of #eligible by school. Schools in this sense would be the "subjects". I would model the pre and post-intervention time periods separately. Any insight is greatly appreciated!


title: Which Regression methods are suitable for binary valued features and continuous output?

text: I want to build a machine learning model to regression on continuous output given binary valued features(0,1). the dimension of my problem is around 200. which of the flowing methods seems suitable for this kind of problem ? SVR with different Kernels Regression random forest MARS Gradient boosting with regression tree Kernel regression (Nadya-Watson Kernel regression) LSR and LARS Stochastic gradient boostin'


title: Implementing a Kolmogorov-Smirnov permutation test

text: I would like some advice on how to implement a Kolmogorov-Smirnov test. **This is a condensed version of my previous question -- any further detail can be provided if needed.** I have data from _n_ experimental trials. Each set consists of two paired time-series signals, for which I am calculating an informative "overlap" measure (a single summary value), based on some threshold value. I calculate this summary value for all _n_ sets of data. I would like to determine whether the distribution of these overlap values is non-random. The distribution of the data is not normal, and having read some details about different comparative tests of distributions, the Kolmogorov- Smirnov test seems most appropriate. I'm hoping someone can give me advice on how to perform a permutation test using these data. My understanding is that this requires comparing the distribution of the actual data with those obtained from random combinations of the data, however, I'm uncertain of the procedure for implementing and interpreting the result. Any help would be much appreciated.


title: Joint distribution of sum of independent normals

text: Suppose we have three independent normally distributed random variables $$ X_0 \sim \mathcal{N}(\mu_0, \sigma_0^2), $$ $$ X_1 \sim \mathcal{N}(\mu_1, \sigma_1^2), $$ $$ X_2 \sim \mathcal{N}(\mu_2, \sigma_2^2).$$ Now, define two new random variables $Y_0 = X_0+X_1$ and $Y_1 = X_1+X_2$. Let $\vec{Y} = [Y_0 \;\;\; Y_1]^T$ What can we say about the distribution of $\vec{Y}$? Obviously, $Y_0$ and $Y_1$ are not independent. If they were, then $\vec{Y}$ would have been a multivariate normal variable. Any ideas?


title: What is this graph called?

text: ![http://i.stack.imgur.com/jGUjB.jpg](http://i.stack.imgur.com/4IGUc.jpg) Sorry for the game-image link, I've wondered what this diamond-shaped graph is called for some time, any idea what type of chart/graph this is?


title: How to determine if distribution in timeseries has shifted?

text: Is it possible to determine if distribution of time-series has shifted significantly without assuming any particular underlying data distribution? Furthermore, would it be possible to identify point of split between (possibly two) distributions? Any links to softwares or pointers to references would be greatly appreciated.


title: Is there a test for independence in a Bernoulli process?

text: I have a process with binary output. Is there a standard way to test if it is a Bernoulli process? The problem translates to checking if every trial is independent of the previous trials. I have observed some processes where a result "sticks" for a number of trials.


title: causal inference with correlated multivariate outcomes

text: I've been struggling with how to think about the causal estimate of a program on two outcomes, when one of the two outcomes affects the other outcome. It seems sort of like simultaneous equations, but reading up on SEM, I don't see immediately how it could apply. Would appreciate any ideas for how to approach this problem. Here is the setup: Two outcomes: $Y_1$ and $Y_2$ One treatment: $T$ Data is a panel, and we believe that $E[\epsilon|\alpha, T] = 0$ The true model is something like: $$ Y_{1,it} = \alpha_i + \delta T_{it} + \gamma Y_{2,it} + \beta (Y_{2,it}\times T_{it})+ \epsilon_{it}\\ Y_{2,it} = \alpha_i + \delta T_{it} + \eta_{it} $$ So, $T$ has some effect on $Y_2$, and both $T$ and $Y_2$ have an effect on $Y_1$. Importantly, $Y_2$ happens _before_ $Y_1$. If I were to estimate the first equation by OLS, $\delta$ would be the average effect of $T$ not accounting for $Y_2$. But that's not really what I'm after -- I want the average effect of $T$ on $Y_1$, accounting for both channels through which $T$ works. Estimating the second equation via OLS is valid (right?) Would it be valid to plug in fitted values of $\hat Y_2$ into the first equation? Probably not, because the standard errors wouldn't account for uncertainty in the second equation... One could sort of justify it algebraically, but I wouldn't know how to interpret the fitted model. This question falls under the heading of "bad controls" -- but perhaps one aspect of it that makes it a bit more difficult than typical "bad control" problems is that I really think that the interaction there is important, and because $Y_2$ isn't determined _entirely_ by $T$. How would I approach this problem? Is this something that can be solved by some well-known technique? If so, which one? Thanks in advance for any advice! **EDIT** : Thinking about it more, I'd get a causal estimate of the effect of $T$ alone by running $Y_{1,it} = \alpha_i + \delta T_{it} + \epsilon_{it}\\$. But I also want an estimate of $\gamma$ and particularly $\beta$. Given that the $\delta$ from the second equation is nonzero, is there some problem with estimating the first equation by OLS? Do the parameters have causal interpretation? If not, how do I get to causal interpretation?


title: Dealing with correlated regressors

text: In a multiple linear regression with highly correlated regressors, what is the best strategy to use? Is it a legitimate approach to add the product of all the correlated regressors?


title: Delta method for a two equation system

text: I want to estimate a nonlinear coefficient from a 2 equations system and use the delta method to calculate the standard errors for inference. \begin{align} y_{i1} &= \alpha_{0} + \alpha_{1}x_{i}+\alpha_{2}y_{i2}+\varepsilon_{i1} \\ y_{i2} &= \beta_{0} + \beta_{1}x_{i}+\beta_{2}z_i+\varepsilon_{i2} \end{align} $y_{i2}$ is endogenous in equation (1), $x_i$ (which appears in both equations) and $z_i$ are exogenous variables. They are used as instruments for equation $y_{i2}$ in equation (1). My interest lies in the combined coefficient : $\alpha_{comb}=\alpha_{1}/\beta_{1}+\alpha_{2}$ The standard errors should be calculated with the Delta method: \begin{eqnarray*} g(\hat{\theta}) \approx g(\theta) + \frac{\partial g(\theta)}{\partial \theta}(\hat{\theta}-\theta) \end{eqnarray*} \begin{eqnarray*} g(\hat{\theta}) \sim N(g(\theta) , \nabla g(\theta)' Var(\theta) \nabla g(\theta) ) \end{eqnarray*} with: $\theta = \begin{pmatrix} \alpha_{0} & \alpha_{1} & \alpha_{2} & \beta_{0} & \beta_{1} & \beta_{2} \end{pmatrix}'$ and: $g(\theta) = \alpha_{1} / \beta_{1} + \alpha_{2}$ \begin{eqnarray} \frac{\partial g(\theta)}{\partial \alpha_{1}} & = & \frac{1}{\beta_1} \\ \frac{\partial g(\theta)}{\partial \alpha_{2}} & = & 1 \\ \frac{\partial g(\theta)}{\partial \beta_{1}} & = & -\frac{\alpha_1}{(\beta_1)^2} \\ \frac{\partial g(\theta)}{\partial \alpha_{0}} & = & \frac{\partial g(\theta)}{\partial \beta_{0}} = \frac{\partial g(\theta)}{\partial \beta_{2}} = 0\\ \nabla g(\theta) & = & \begin{pmatrix} 0 & \frac{1}{\beta_1} & 1 & 0 & -\frac{\alpha_1}{(\beta_1)^2} & 0 \end{pmatrix}' \end{eqnarray} standard error of $g(\theta) = \alpha_{1} / \beta_{1} + \alpha_{2}$: \begin{eqnarray*} SE(g(\theta)) & = & \sqrt{\nabla g(\theta)' Var(\theta) \nabla g(\theta) } \end{eqnarray*} I think this should be correct, but Im not 100% sure. Do you think there is an (1) a mistake or (2) Delta method is not appropriate for this subject and maybe I should bootstrap the coefficient?


title: Good, authoritative recent book on factor analysis and principal component analysis

text: Questions regarding factor analysis and principal component analysis, both here and elsewhere, reveal considerable confusion of terms, disagreements about requirements, different recommendations for different issues and so on. Different programs label the same output differently. There is a fairly bewildering number of choices of extraction, rotation and so on. In addition, there is confusion about exploratory vs. confirmatory factor analysis and how both differ from principal component analysis. I studied factor analysis way back in the late 1990s. I've read a bunch of recent articles. I think I more or less know what I am doing. However, it would be good to know of some authoritative recent source; perhaps an edited collection. Some quick searching on Amazon didn't reveal any obvious candidates; are there any?


title: Probabilities from Graphical Model

text: I am trying to derive some probabilities from the following graphical model and I am having some difficulty. The model is actually a complete bipartite graph with two sets of nodes C and W. The graph is a directed graph where each $c_i$ connects to each $w_i$. That is all of the connections. I am interested in calculating the following probabilities: $$ P(w_i \vert w_j) $$ $$ P(w_i, c_k | w_j, c_l) $$ My approach to the first is: $$ P(w_i \vert w_j) = \sum_j P(w_i, c_j \vert w_j) = \sum_j P(w_i \vert c_j, w_j)P(c_j \vert w_j)) = P(w_i \vert c_j)P(c_j \vert w_j) $$ My approach to the second: $$ P(w_i, c_k \vert w_j, c_l) = P(w_i \vert c_l, c_k, w_j)P(c_k \vert w_j, c_l) = P(w_i \vert c_l, c_k, w_j)P(c_k \vert w_j) $$ First of all I am not sure if these are right. The second issue is I would like to eliminate the dependence on $w_j$ in the equation for the second case. Is this correct? Any suggestions? I am concerned about the step in the first equation when I remove the dependence on $w_j$.


title: List of sample arrangements at a given Wilcoxon / Mann-Whitney test statistic

text: Consider a Wilcoxon-Mann-Whitney test statistic, either in the form (i) "W = sum of ranks in sample A"; or in the form (ii) "sum of ranks in sample A - the minimum possible sum of ranks in sample A", or equivalently "U = the number of times a value in sample 1 exceeds a value in sample 2". **How would I generate a list of all possible sample arrangements (e.g. the specific ranks in sample A or some equivalent from which that could be computed) at a given value of the statistic (in either form)?** If I have $n_A=3, n_B=4$, I might specify either $W=10$ or $U = 4$ (the same thing, since $W = U + 6$ as the value I want to see all samples for; a solution for either is fine. As an example of what I seek, at $n_A=3, n_B=4$ the possible ranks in sample A are:                W      6: {1, 2, 3},      7: {1, 2, 4},       8: {1, 2, 5}, {1, 3, 4},       9: {1, 2, 6}, {1, 3, 5}, {2, 3, 4},      10: {1, 2, 7}, {1, 3, 6}, {1, 4, 5}, {2, 3, 5},      11: {1, 3, 7}, {1, 4, 6}, {2, 3, 6}, {2, 4, 5},      12: {1, 4, 7}, {1, 5, 6}, {2, 3, 7}, {2, 4, 6}, {3, 4, 5},      13: {1, 5, 7}, {2, 4, 7}, {2, 5, 6}, {3, 4, 6},      14: {1, 6, 7}, {2, 5, 7}, {3, 4, 7}, {3, 5, 6},      15: {2, 6, 7}, {3, 5, 7}, {4, 5, 6},      16: {3, 6, 7}, {4, 5, 7},      17: {4, 6, 7},      18: {5, 6, 7}      I'd like a way to specify "nA=3, nB=4, W=10", (or an equivalent), say and get   "`{1, 2, 7}, {1, 3, 6}, {1, 4, 5}, {2, 3, 5}`" or something equivalent to it back. I can come up with brute force approaches of course but I don't think that's feasible in the larger sample sizes. [For my particular application, it happens that $n_A=n_B$ and I'm mainly interested in cases between about $n_i=10$ to perhaps $n_i=40$ or so, so if either fact makes life significantly easier, feel free to use it. If it's doable at say samples of 25 but not at 40, that would still be useful.] R packages or code would be especially helpful, but I can code, so an algorithm would be fine and would probably be more generally useful. Indeed, even some pointers would help.


title: How to simplify coxme survival models?

text: I have some questions about specifying a `coxme` (mixed-effects Cox proportional hazards) model in R and then simplifying it after reading the `coxme` documentation. The situation is this: I glued 4 pieces of small coral fragments onto small ceramic tiles, which I placed at 4 distances east and west of 10 large coral colonies (i.e. site). Thus, each tile represents one distance-direction-site combination. I checked the small coral fragments daily to see which had died overnight and at the end of the experiment some were still alive (thus, censored). I therefore had 4 fragments per tile*4 distances*2 directions*10 sites = 320 small fragments. Distance and direction are fixed effects, while the tile that each fragment is on and the site are random effects. In addition, each large colony is a different size, so the size of the large colonies should be a random effect, too (`SiteSize`). The R model I wrote to express this is:               mefull<-coxme(Surv(death, censor) ~ Distance*Direction+(1|Site/Tile) +(1|SiteSize))        * First, can anyone tell me if this properly specifies the situation I described above? After running this model, I found that neither fixed effect nor their interaction was significant. Also, the standard deviation for `Site` and `SiteSize` are identical (~1.12), which seems strange to me. Is there a reason for that? The fact that they are both greater than 1 indicates to me that they contribute a lot of variation to survival. Is that correct?   * My next major question is how to simplify this model. My instinct (and based on reading the manuals and posts here) is to remove each random effect in turn and compare the AICs of the integrated log-likelihood of the resulting models; the model with higher AIC is preferred in this formulation. Is that correct?     * However, I'd also like to try to simplify the model through removal of the non-significant fixed effects, starting with their interaction. How can I do this while also removing random effects?      * What terms should I start with removing, or does the order not matter as long as I start with higher-order terms (i.e. interaction)?      * Can I try as many combinations as I like or do issues with multiple tests come into play? Some options are removing one of the random effects (me2) or removing the interaction between the fixed effects but keeping the random effects in place (`me3`).              me2<-coxme(surv ~ Distance*Direction+(1|Site/Tile))     me3<-coxme(surv ~ Distance+Direction+(1|Site/Tile)+(1|SiteSize))      When I run these and other combinations of factors, their AIC is always lower than that of the full model, which suggests to me that the full model is best. Any guidance on how to simplify this model would be greatly appreciated.   * Finally, to compare the model with random effects to one without, can I compare the NULL log-likelihood with the integrated log-likelihood? From my understanding of the coxme manual, the one closer to 0 is the better model, so if the integrated one is closer to 0 then the model with random effects is preferred over the one without random effects.  Thanks very much for your time and help.


title: Naive Bayes with unbalanced classes

text: As a part of a project for the university is should train a Naive Bayes classifier to classify question and answers in three different categories, the task should be easy since that the 3 classes are really different between each other. **Dataset** the dataset is a mixture of questions and answers from different domains (furniture and a C++ course) and they are in different language (Italian and English) so at the first sight it should be easy to classify them. The only problem with the dataset is that is really unbalanced, like C++ course 2700 training instance. furniture English 200 training instance. furniture Italian 60 training instance. **Feature Extraction** The feature are simply n-gram counts i don't remove stopword because i have 2 different language to work with. Using TF-IDF features and stemmed token i obtained lower results. **The NB algorithm** I have implemented the naive bayes by myself but it obtains the same result of the scikit learn one. i have trained it with per class prior and a smoothing using alpha=.5 **The results** The result at the end was in some sense good                             precision    recall  f1-score   support                0           0.92      1.00      0.96       318           1           1.00      0.50      0.67        44           2           1.00      0.33      0.50         6          avg / total       0.93      0.93      0.92       368      But the only draw back is that the recall on the 1 and 2 class is low, and the reason is simple (we have unbalanced classes). There is a way to actually obtain lower but more balanced results having the constraint of using a single model trained with Naive Bayes? Cheers * * * **Edit:** My needs of having more balanced results is due to the fact that I should use this classifier as part of an Question Answering task.


title: Get rid of the influence of a predictor

text: I am currently working on a data mining project during an in internship, and I use among others methods decision trees. My problem is that I have a categorical predictor very influent on classification, but I want to ignore this influence to better analyze the behaviour of other predictors. I could of course delete this predictor from the algorithm but the tree would be biased by this invisible parameter. I have done a Naive Bayes model of the influence of this predictor but I don't know how to use this to debiased the tree. Any idea would be appreciated Thanks in advance, Pierrick


title: How to sample from a multivariate normal given the $P^T LDL^T P$ decomposition of $\Sigma$?

text: How to sample from a multivariate normal given the $P^TLDL^TP$ decomposition of $\Sigma$, and $\mu$? Here, $P$ is a permutation matrix, $L$ is lower-triangular, and $D$ is diagonal Given the Cholesky decomposition $LL^T$ of $\Sigma$, and $\mu$, we can sample by doing: $sample = \mu + L * z$ ... where $z$ is a vector of univariate standard normals. I imagine that there must be some way of tweaking this formula to work with the $P^TLDL^TP$ decomposition, but my maths isn't quite good enough to see how to do that?


title: Exact p-values for real world outcomes

text: I had an interesting exam question yesterday, in which you had the outcomes of number people killed in traffic accidents in two subsequent years (the second being relatively higher than the first). Now you should calculate the **exact p-value** for the hypothesis that the second outcome was due to a chance event. The idea was that you should assume that the number of fatalities in traffic accidents was Poisson distributed then you should test for the parameter-ratio being 1. But is this really an "exact p-value" for this hypothesis. Since the Poisson distribution is just an approximation this is neither an exact p-value in theory nor in practice. How do you more experienced statisticians look at this?


title: Conditional inverse Gaussian distribution

text: I am considering the Inverse Gaussian distribution as the hitting time distribution for a Wiener process, $W(t)$, with drift parameter $\nu$ and variance parameter $\sigma$. Define the hitting time as the first time the process reaches the level $a > 0$: $T = \min_{t>0}\{ W(t) \geq a\}.$ The hitting time, $T$, of the level $a$ is then Inverse Gaussian. Let $\mu = a/\nu$ and $\lambda = a^2/\sigma^2$: $f_T(t) = \sqrt{\frac{\lambda}{2 \pi t^3}} \exp\Big\{- \frac{\lambda(t-\mu)^2}{2\mu^2t}\Big\}$ $P(T \leq t) = F(t) = \Phi\Big(\sqrt{\frac{\lambda}{t}}\Big(\frac{t}{\mu} - 1\Big)\Big) - \exp\Big(\frac{2\lambda}{\mu}\Big)\Phi\Big(- \sqrt{\frac{\lambda}{t}}\Big(\frac{t}{\mu} + 1\Big)\Big),$ where $\Phi(\cdot)$ is the standard normal distribution CDF. I am looking to find the conditional Inverse Gaussian distribution, conditioned on the endpoint value. Say the endpoint value, $W(t)$ = $r < a$ is known, I want to find the probability that the threshold has been crossed: $P(T \leq t \big| W(t) = r < a)$ What I have done leaves me with quite nasty expressions. Any help on approximations as well as insight on simplifications of the expressions is appreciated. What I have done is: for $r < a$ $P(T \leq t \big| W(t) = r) = \frac{P(T \leq t, W(t) = r)}{P(W(t) = r)} $ Then, for the sake of notation, consider the problem "discretely". I say that for each $0 < s < t$, the joint probability can be found by considering the hitting time $s$ and the increment from $s$ to $t$, $W(t) - W(s),$ which must be $r - a$: $ P(T = s) P(W(t) - W(s) = r - a)$, so I must sum over $s$. In continuous notation: $P(T \leq t, W(t) = r) = \int_{0}^t f_T(s) f_{W(t) - W(s)}(r-a) ds$ Note that $W(t) - W(s)$ is normally distributed with mean $(t-s)\nu$ and variance $(t-s)\sigma^2$. Finally in my problem I must sum over all $0 < r < a$ which leaves me with a double integral, to integrate out $W(t) = r$ on which it was conditioned in the beginning. $\displaystyle\int_{r=0}^a\frac{\int_{0}^t f_T(s) f_{W(t) - W(s)}(r-a) ds}{f_{W(t) - W(0)}(r)} dr.$ How can I simplify this?


title: Check if sample is representative of a larger sample

text: From a larger sample of tabular data, I have picked certain rows that meet a certain condition (this condition is unrelated to the actual data in the rows). Now, I want to know if the distribution of this subset that I have created is similar to the distribution of the original, larger sample. What test(s) can I use for this purpose? Thanks! I appreciate the help.


title: Variance-covariance matrix as the sum of variance covariance matrices

text: I have a variance-covariance matrix, $\mathrm{V}$. This allows me to take a vector, $x$ of independent random variables drawn from a known distribution, and induce a required variance-covariance structure to give a correlated random variable $y$ from: $y = \mathrm{V}^{\frac{1}{2}} x$ Because $\mathrm{V}$ is positive definite, I can find $\mathrm{V}^{\frac{1}{2}}$ using the Cholesky decomposition. However, I know that the variance-covariance structure represented by $\mathrm{V}$ actually comprises three separate components - say $\mathrm{V}_1$, $\mathrm{V}_2$ and $\mathrm{V}_3$. It would be highly useful to show the impact of these three separate components. Is there any way to do this?


title: Should I use independent samples t-test?

text: I need to analyse if there's a difference between poor and rich European countries (using a GDP per capita index variable) and look if it has an effect on the variable environment. The GDP per capita variable is a ratio variable going from 44-..., while the environment variable is a nominal variable with two categories (environment is important and environment isn't important). I was thinking about using the independent t-test since there's a scale variable and a nominal variable with two categories. What's your opinion? If I read the question right, I need to know the effect of GDP (scale variable) on environment (nominal variable). This environment variable is a recoded variable (before it was what is the most important thing for your country: economics, crime, health ...) so before I recoded it, it had 15 categories or so. So I decided to recode it in 2 categories those who thought that environment was the most important and those who didn't. There's also a Likert scale variable that I could use instead of this 'environment' variable, but I think it would be too difficult. This is the GDP variable descriptive:               GDP_Capita_Index     ----------------                 Frequency Percent Valid Percent Cumulative Percent Valid       44,00        1.003         6,1         6,1          6,1         51,00        1.000         6,1         6,1         12,2         74,00        1.159         7,1         7,1         19,3         80,00        1.428         8,7         8,7         27,9         85,00        1.082         6,6         6,6         34,5        100,00        2.560        15,6        15,6         50,1        112,00          928         5,6         5,6         55,8        115,00        1.211         7,4         7,4         63,1        118,00        1.407         8,6         8,6         71,7        119,00        1.142         7,0         7,0         78,7        123,00        1.181         7,2         7,2         85,9        126,00        1.019         6,2         6,2         92,1        127,00        1.305         7,9         7,9        100,0        Total        16.425        100,0      100,0


title: Interpretation of result in reference group

text: If I am looking at sports data of a whole league over a year:               Dependent Variable = Goals Scored          Independent Variables:     X1 = Traditional Home Uniforms     X2 = Retro Home Uniforms     X3 = Secondary Home Uniforms           Controls:     A whole bunch      My objective is to see what the effect of different uniforms are on the number of goals scored. More specifically, the effect of Retro versus Secondary Uniforms. When I run an OLS I get these results:               Model 1: X1 = Reference Group     X2 = -15.9 (p = 0.01)     X3 = -7.6 (p = 0.01)          Model 2: X3 = Reference Group     X1 = 7.6 (p = 0.01)     X2 = -6.2 (p = 0.34)      Is it correct to say (from Model 1) that compared to Traditional Uniforms that wearing Retro Uniforms is better for a team than Secondary Uniforms, ceteris paribus. It would be even better if it came through in Model 2 because that is a direct comparison of Retro Uniforms versus Secondary Uniforms. All models use the same controls. My questions: 1) Is it correct to say what I am saying? 2) Why are the results 'better' for what I am saying when I use Model 1 over Model 2?


title: Strange result of a fixed effect model

text: If I am looking at sports data of dozens of soccer leagues over a year:               Dependent Variable = Goals Scored          Independent Variables:     X1 = Traditional Home Uniforms     X2 = Retro Home Uniforms     X3 = Secondary Home Uniforms           Controls:     A whole bunch     Most importantly--Every team has a rank (1-50) from the previous year       When I run my regressions clustering my errors on rank the results are fine, but when I add a dummy for each rank my output becomes skewed:               Number of obs =    972     F( 18,   241) =       .     Prob > F      =       .     R-squared     =  0.2478      Any ideas why my F becomes a '.'?


title: Self Study on Simple Linear Regression (Not So Simple Afterall)

text: I have been practicing some simple linear regression exercises lately and have been deeply challenged by the following question. While the part (A) is something I can concur, the proceeding parts baffled me in many ways. In general, what I am challenged with is the (Value, Value) expressions. Appreciate if I will be able to be enlighten by what (B) means and shed some light on the rationale behind the provided answer, so that I may move on to the other parts. Thank you so much once again **Question** ![enter image description here](http://i.stack.imgur.com/UoKKw.jpg) **Data Background** ![enter image description here](http://i.stack.imgur.com/oJfQe.jpg) **Model Answer** ![enter image description here](http://i.stack.imgur.com/OXyq3.jpg)


title: spss principal axis factoring output problem

text: I am trying to run a factor analysis (SPSS) using principal axis factoring with an oblique (or promax) rotation because the variables are highly correlated. However, the output for "total variance explained" does not give any information on "extraction sum of squares" and "rotation sums" as it does if I use principal components (the table only shows eigenvalues). Has anyone encountered this issue? Does it have to do with PAF itself, the program, or something else that I am missing? I tried removing variables and running with fewer, but didn't change anything, tried without a rotation and didn't work again. Is it possible to use a general principal components to extract the factors, and then for individual factors to use PAF? Thank you :)


title: What is the meaning of operators in regression or anova formulas in R

text: I was wondering, what is the meaning of operators in anova or regression formulas in R For example   * " **+** " aov <- aov(x~time+sample, data=data) -> repeated mesures anova?   * " ***** " aov <- aov(x~time*sample, data=data) -> two way anova?   * " **/** " aov <- aov(x~time/sample, data=data) -> ?   * " **:** " aov <- aov(x~time:sample, data=data) -> ? And also are there more operators for this kind of formulas?


title: Arbitrariness of Euler's number in exponential of log-normal distribution

text: The Wikipedia article of the log-normal distribution says > If $X$ is a random variable with a normal distribution, then $Y = \exp(X)$ > has a log-normal distribution; likewise, if $Y$ is log-normally distributed, > then $X = \log(Y)$ is normally distributed. (This is true regardless of the > base of the logarithmic function: if $\log_a(Y)$ is normally distributed, > then so is $\log_b(Y)$, for any two positive numbers $a, b ≠ 1$.) (The general form of this has described the log-normal distribution on Wikipedia since at least 2006.) I have two questions. 1) is it the case that the logarithm can be taken to any base (i.e., it need not be the natural log), while the exponential must have Euler's number $e$ as the number which is raised to a power? Or, is it the case that I could use other bases to describe the log normally distributed variable, $X = a^{\mu+Z\sigma}$, where $Z$ is a standard normal and $a$ is any positive, real number? And, conditional on the answer to #1, 2) why the asymmetry?


title: Normalization of dummy variables

text: My data consists of several continuous measurements and some dummy variables representing the years the measurements have been made. Now, I want to learn a neural network with the data. Therefore, I am zScore-normalizing all variables, including the dummy variables. However, I wonder if this is a reasonable approach, because normalizing the dummy variables alters their ranges, which I guess makes them less comparable if their distributions differ. On the other hand, not normalizing the dummy variables might also be questionable, because without normalization their influence on the networks output might be suboptimal. What is the best approach to deal with dummy variables, normalizing them (zScore) or just leaving them as they are?


What is the difference between linear regression and logistic regression? When would you use each?


title: Estimating two break points in a broken stick model with random effects in R

text: Five months ago, jbowman posted a very useful answer to estimate the break point in a broken stick model with random effects in R. I never use "computing" like `ifelse` and I would like to estimate two break points. I should write two other functions like `b1` and `b2` but I don't know how. Can someone please tell me how to do that in R? Thanks! jbowman's code:               library(lme4)     str(sleepstudy)          #Basis functions     bp = 4     b1 <- function(x, bp) ifelse(x < bp, bp - x, 0)     b2 <- function(x, bp) ifelse(x < bp, 0, x - bp)          #Wrapper for Mixed effects model with variable break point     foo <- function(bp)     {       mod <- lmer(Reaction ~ b1(Days, bp) + b2(Days, bp) + (b1(Days, bp) + b2(Days, bp) | Subject), data = sleepstudy)       deviance(mod)     }          search.range <- c(min(sleepstudy$Days)+0.5,max(sleepstudy$Days)-0.5)     foo.opt <- optimize(foo, interval = search.range)     bp <- foo.opt$minimum     bp     [1] 6.071932     mod <- lmer(Reaction ~ b1(Days, bp) + b2(Days, bp) + (b1(Days, bp) + b2(Days, bp) | Subject), data = sleepstudy)


title: Split dataset randomly

text: I have a database with 500 records. I want to split these records to 75% and 25% *randomly*in order to use the different datasets for training and testing to machine learning algorithms. Does anyone knows how to do that? For example using an sql query


title: How to correctly model noise?

text: Assume a linear mixing model $x = As$, where $x = (x_{0}, ..., x_{n})^T$ are linear mixtures of $s = (s_{0}, ..., s_{n})^T$, and $A$ is the mixing matrix. Now, if I introduce additive noise to this system, that is, $x = As + e$, where $e = (e_{0}, ..., e_{n})$, what are the relations between $e_{0}, ..., e_{n}$? For example, assume I have three microphones ($x_{0}, x_{1}, x_{2}$) in a room. There are also three speakers ($s_0, s_1, s_2$) and some background noise ($e$). If we assume the noise is white gaussian (or laplacian, or uniform), and that the noise level is static over the room, do I get the same noise into each of the microphones, that is $e_0 = e_1 = e_2$? If not, why?


title: What to do with a markovian process with unknown path length?

text: I have a markovian process with an initial state distribution $p(Z_0)$; i want to predict $Z_i$ for a fixed $i$ and i observe $Z_{i+j}$. However, the catch is that $j$ is also a hidden variable with a distribution $p(j)$. How can i effectively compute $p(Z_i| Z_{i+j} )$? Some additional information about the problem:   * The process is deterministic but complex   * $p(j)$ is only non-zero for a (small) finite number of values   * $p(j)$ follows a uniform distribution


title: Is it possible to construct an empirical "distribution" from empirical data?

text: I am trying to compare two (univariate) data sets using the `HellingerDist` function described here. However, this package only allows the a `numeric` type with one of the distribution classes. Is it possible (and does it even make sense) to coerce one of my data sets into a distribution so I can use that function to compute the distance?


title: Measuring distance between two empirical distributions

text: This is similar to the question Measuring the "distance" between two multivariate distributions, except that I want to measure distance between two data sets. I can imagine simply computing the mean, standard deviation, and other summary statistics and then somehow aggregating the differences for each distribution, but that does not seem very theoretically sound. I wonder if there is a generally accepted method for doing this measurement.


title: R or python implementation of sparse PCA for p>n

text: According to this paper, there are 2 algorithms to perform sparse PCA. One is better if $p>n$. I need to run SPCA on a $2000\times12000$ matrix so I am looking for an implementation of this algorithm. There is an sparse PCA implementation in sklearn and one in the R package `elasticnet` but I think they are the $n>p$ version. Where can I find an implementation of this algorithm?


I am testing two path analysis models. Both models have the same predictors. The outcome variable is all that differs   * Model 1: Outcome variable is intention to use a face-to-face support service   * Model 2: Outcome variable is Intention to use a telephone support service When I tested the models, the fit indices for both indicated a very good fit to the data. On the second model, however, for one of the variables, no significant paths are evident (to the outcome variable or any other IV). I would have expected the relationships to be significant based on correlations. When I remove this predictor variable from the model, the fit is no longer acceptable. **Why does the fit get worse when I remove a non-significant predictor from my path model?**


title: True standard deviation (hypothesis testing)

text: I took two statistics classes at college and I remember "hypothesis testing" to be a very powerful tool if only I could apply it correctly. It seems that while I was in college solving problems, I was always given the population's standard deviation to solve the question. However, what if I want to apply it in my job (i.e. track a particular process)? How can I obtain the true standard deviation of a process? I am interested in normally distributed data. Let's say that I work at an industry that designs sprinkles, and I was interested to know if the product meets the specifications. The sprinkles were designed so that the average activation response time should be 25 (in a fire prevention system). Let's say I analyzed the times until 10 sprinkler samples activate in response to heat, and I got 27,41,22,27,23,35,30,33,24,27 (in seconds). Now I made the hypothesis test if the mean is equal to 25 or not equal to 25.


title: Steps to calculate the MGF of (1+sum of gamma distribution) over gamma distributions

text: What is the best way or steps to calculate the moment generating function (mgf) of (1+sum of gamma distribution) over gamma distributions? Is it the same way to derive beta distributions? I am trying to find the mgf of $$\exp[-(1+\sum X_j)/Y]$$ where $X$ are gamma distributions that have the same scale factor, say $m,$ and $Y$ also is a gamma distribution with scale factor $n.$ What I did is split it into $$\exp[-1/Y]\exp[\sum X_j /Y].$$ I managed to find the mgf of $\exp[-1/Y]$ which is $\frac{2}{\Gamma\left(m\right)}\left(ms\right)^{m/2}K\left(2\sqrt{ms}\right)$, but I am not sure what steps to do to find the mgf of the sum of $X$ over $Y$ ratio.


title: What type of sampling is this?

text: I have some confusion about the kind of sampling I'm using. I started by choosing the manufacturing sector and within it, chose 4 industries. Later I distributed questionnaires within those 4 industries as per my convenience (convenience sampling). What kind of sampling am I using?


title: Normal Bayesian Model: Marginal distribution of $\bar X$ with unknown mean and unknown variance

text: For $i=1, \ldots, K$ and $j=1, \ldots,n$, assume the following model. \begin{align} X_{ij} \mid \mu_i, \sigma^2 & \stackrel{_\text{ind}}{\sim} N(\mu_i, \sigma_i^2) \\ \mu_i & \stackrel{_\text{iid}}{\sim} N(\mu, \tau^2) \\ \ln(\sigma_i^2) & \stackrel{_\text{iid}}{\sim} N(\mu_v, \tau_v^2) \end{align} Is it possible to obtain the marginal distribution of $\bar X_i$ under this model? Thank you very much.


title: Comprehensive dataset for documents classification

text: In document classification (document categorization) field, researchers recommend always a few standard dataset, such as Reuters-21578, RCV1(Reuters Corpus Volume 1). However, these datasets only contain documents related to specific domain, for instance, **newswire stories**. However, if we really want to classify any given documents (1 million) in a more wide range of topics, like literature, sports, politics, science, etc. How can we accomplish this task? If the classification (supervised learning) doesn't work, could anyone tell me some more advanced methods to automatically categorize any topic documents ?


title: Simulating distributions

text: I am working on a Capacity Planning assignment and I have read some books. This is specifically about distributions. I use R.   1. What is the recommended approach to identify what my data distribution is ? I have this diagram but are there statistical methods to identify it ? ![PROBABILISTIC APPROACHES: SCENARIO ANALYSIS, DECISION TREES AND SIMULATIONS](http://i.stack.imgur.com/NYb5t.jpg)   2. What are simulations approaches available using R ? Here I want to generate data for a certain distribution like exponential. Is r-java the right approach if I want to integrate it with Java ?    3. Is there a way to predict what distribution the effect( CPU usage etc. ) will have when I pipe data for a particular distribution ? What are the different effects of sending certain distributions of data ? Please consider these as beginner's questions. Are there books or material that deal with these types of simulations ? ## Notes THe diagram is from the end of the paper http://people.stern.nyu.edu/adamodar/pdfiles/papers/probabilistic.pdf. ### Goodness of fit techniques I have come across Assessment of goodness-of-fit   1. Chi-squared   2. Kolmogorov-Smirnov,   3. Anderson-Darling statistics density, cdf, P-P and Q-Q plots I am not sure what the interpretation or next steps should be if I find that my distribution is normal or exponential etc. What does it allow me to do? Prediction? Hope this question is clear. Exponential delays will induce queue fluctuatuions as per my Capacity Planning book by Neil Gunther. So I know that one point.


title: How to calculate cumulative distribution in R?

text: I need to calculate the cumulative distribution function of a data sample. **Is there something similar to hist() in R that measure the cumulative density function?** I have tries ecdf() but i can't understand the logic.


title: Binomial probability homework questions

text: 38% claim to have at least one tattoo. If you asked 15 people, what is the probability that none of them have a tattoo? What is the probability that at least 10 of them have a tattoo?