title: Everything else held constant, what is the relationship between alpha and beta?

text: **I want the mathematical formula the statistical programs are using to compute alpha and beta.** I typically use G-Power to estimate the necessary sample size for a given power (1-beta) for an a-priori analysis. The input it requires for this calculation is: Type of analysis (t-test, difference between two independent samples) the number of tails (2) alpha (0.05) power (0.8) effect size (0.2) = estimated necessary sample size (788) I have looked everywhere I know to find out what the function/formula is that is being used to compute the sample size value from these variables. I know that:   * alpha=type I error=FP/(FP+TN)    * beta=type II error= FN/(TP+FN)    * where FP=false positive    * TP=true negative    * FN=false negative    * TN=true negative But I don't know how to compute false positives or true negatives a-priori. Unless Stevens (2009) is wrong, power (the inverse of beta) is a function of alpha therefore β ∝ α (e.g. if, {α = .1, β = .37, power = .63}, {α = .05, β = .52, power = .48}, {α = .01, β = .78, power = .22}, see p.4). I am simply interested in finding the function that describes the relationship between alpha and beta/power. My question: Do you know where I can find the function to describe this relationship? The reason? To actually understand what I am doing every time I run a statistical analysis. when I perused the source code for an R package pwr I found: What looks like the n is obtained numerically using a root finding algorithm computed with the highlighted part of the function below (uniroot looks like a general function to find where a supplied function equals 0).               > pwr.t.test     function (n = NULL, d = NULL, sig.level = 0.05, power = NULL,          type = c("two.sample", "one.sample", "paired"), alternative = c("two.sided",              "less", "greater"))      {         if (sum(sapply(list(n, d, power, sig.level), is.null)) !=              1)              stop("exactly one of n, d, power, and sig.level must be NULL")         if (!is.null(sig.level) && !is.numeric(sig.level) || any(0 >              sig.level | sig.level > 1))              stop(sQuote("sig.level"), " must be numeric in [0, 1]")         if (!is.null(power) && !is.numeric(power) || any(0 > power |              power > 1))              stop(sQuote("power"), " must be numeric in [0, 1]")         type <- match.arg(type)         alternative <- match.arg(alternative)         tsample <- switch(type, one.sample = 1, two.sample = 2, paired = 1)         ttside <- switch(alternative, less = 1, two.sided = 2, greater = 3)         tside <- switch(alternative, less = 1, two.sided = 2, greater = 1)         if (tside == 2 && !is.null(d))              d <- abs(d)         if (ttside == 1) {             p.body <- quote({                 nu <- (n - 1) * tsample                 pt(qt(sig.level/tside, nu, lower = TRUE), nu, ncp = sqrt(n/tsample) *                      d, lower = TRUE)             })         }         if (ttside == 2) {             p.body <- quote({                 nu <- (n - 1) * tsample                 qu <- qt(sig.level/tside, nu, lower = FALSE)                 pt(qu, nu, ncp = sqrt(n/tsample) * d, lower = FALSE) +                      pt(-qu, nu, ncp = sqrt(n/tsample) * d, lower = TRUE)             })         }         if (ttside == 3) {             p.body <- quote({                 nu <- (n - 1) * tsample                 pt(qt(sig.level/tside, nu, lower = FALSE), nu, ncp = sqrt(n/tsample) *                      d, lower = FALSE)             })         }         if (is.null(power))              power <- eval(p.body)         else if (is.null(n))              n <- uniroot(function(n) eval(p.body) - power, c(2 +                  1e-10, 1e+07))$root             else if (is.null(d)) {                 if (ttside == 2) {                     d <- uniroot(function(d) eval(p.body) - power, c(1e-07,                          10))$root                 }                 if (ttside == 1) {                     d <- uniroot(function(d) eval(p.body) - power, c(-10,                          5))$root                 }                 if (ttside == 3) {                     d <- uniroot(function(d) eval(p.body) - power, c(-5,                          10))$root                 }             }             else if (is.null(sig.level))                  sig.level <- uniroot(function(sig.level) eval(p.body) -                      power, c(1e-10, 1 - 1e-10))$root         else stop("internal error")         NOTE <- switch(type, paired = "n is number of *pairs*", two.sample = "n is number in *each* group",              NULL)         METHOD <- paste(switch(type, one.sample = "One-sample", two.sample = "Two-sample",              paired = "Paired"), "t test power calculation")         structure(list(n = n, d = d, sig.level = sig.level, power = power,              alternative = alternative, note = NOTE, method = METHOD),              class = "power.htest")     }     <environment: namespace:pwr>      I'm just not savvy enough with the notation to transform this into something prettier e.g. f(x)=...


title: isn't it possible giving weight to attributes in Clementine?

text: I use SPSS Clementine. I need to give a weight to one of my features to be considered much more than others as it is more important. I can't find any option to do this for training a C5 tree. the only option I saw is weightening some records, and not the attributes. does any one know how can I do it? regards


title: Time series forecasting using R

text: I have many time series(retail data). Some with trends, some seasonal, and some with neither. With period day, week or month. I need to make forecast, for each time serie. I'm looking for the most efficient methods for forecasting in R ? Which significant things should I know for it? Maybe someone has experience with random forest forecasting and would share with me? Any help would be truly appreciated. UPDATE 1: For example, one of mine time series is x:                  > dput(x)      c(1.07328072153326, 1.07385697538101, 1.10947204968944, 1.10501567398119,      1.08808510638298, 1.07468423942889, 1.06658878504673, 1.10157194679565,      1.10297619047619, 1.09510682288077, 1.07372549019608, 1.08457943925234,      1.09101316542645, 1.10577472841624, 1.08926553672316, 1.0929326655537,      1.08484848484848, 1.09699769053118, 1.10987124463519, 1.08726673984632,      1.09157959434542, 1.10070384407147, 1.08625486922649, 1.11432506887052,      1.0828313253012, 1.08040626322471, 1.07157157157157, 1.08369098712446,      1.08045977011494, 1.10748560460653, 1.11616161616162, 1.08371040723982,      1.10213414634146, 1.06835306781485, 1.07926829268293, 1.08721886999451,      1.10216718266254, 1.1241610738255, 1.08231707317073, 1.07698961937716,      1.08569953536396, 1.09771181199753, 1.07181984175289, 1.07288828337875,      1.07820419985518, 1.07210031347962, 1.07450628366248, 1.06662870159453,      1.07235494880546, 1.0979020979021, 1.08494690818239, 1.06716417910448,      1.08305369127517, 1.08023307933662, 1.07635746606335, 1.07701786814541,      1.08310249307479, 1.0768253968254, 1.096, 1.06787687450671, 1.07353535353535,      1.11226993865031, 1.07641196013289, 1.08066298342541, 1.09431605246721,      1.06678539626002, 1.06646525679758, 1.09977728285078, 1.07646420824295,      1.0973341599504, 1.0906432748538, 1.09831824062096, 1.09302325581395,      1.08199121522694, 1.073753605274, 1.0616937745373, 1.07997481108312,      1.08239202657807, 1.08798283261803, 1.07748776508972, 1.0552611657835,      1.0817746846455, 1.08978032473734, 1.08414985590778, 1.08205756276791,      1.11405835543767, 1.11866969009826, 1.07441154138193, 1.09642703400775,      1.07393209200438, 1.08049535603715, 1.09371428571429, 1.09732824427481,      1.10526315789474, 1.11575091575092, 1.08680994521702, 1.10028929604629,      1.09176340519624, 1.07464266807835, 1.10190664036818, 1.08295281582953,      1.08928571428571, 1.09341998375305, 1.0958605664488, 1.07885714285714,      1.07466814159292, 1.09463722397476, 1.07281903388609, 1.0812324929972,      1.08226102941176, 1.07101616628176, 1.08390410958904, 1.08528528528529,      1.09333333333333, 1.08073929961089, 1.09380234505863, 1.08012968967114,      1.07717391304348, 1.07066508313539, 1.06838106370544, 1.07199032062916,      1.08235294117647, 1.08157524613221, 1.11082474226804, 1.08620689655172,      1.08299477655252, 1.10016420361248, 1.10140093395597, 1.08766485647789,      1.10094850948509, 1.13925191527715, 1.11293859649123, 1.12204234122042,      1.10141364474493, 1.11103495544894, 1.09365558912387, 1.10044313146233,      1.11116279069767, 1.11053240740741, 1.09810671256454, 1.09899823217443,      1.10986101919259, 1.09649805447471, 1.08765778401122, 1.09922928709056,      1.07868303571429, 1.07439104674128, 1.08457374830852, 1.09739714525609,      1.0873440285205, 1.07574536663981, 1.10498812351544, 1.11056105610561,      1.09443402126329, 1.09200240529164, 1.1076573161486, 1.10090237899918,      1.09986225895317, 1.10569105691057, 1.09090909090909, 1.10409356725146,      1.107, 1.15349143610013, 1.08992562542258, 1.09016393442623,      1.08549783549784, 1.07950780880265, 1.08859223300971, 1.06225680933852,      1.08606557377049, 1.07929176289453, 1.09641873278237, 1.07554585152838,      1.05761316872428, 1.08054085831864, 1.09245172615565, 1.09028727770178,      1.06859756097561, 1.08278388278388, 1.06620808254514, 1.07001522070015,      1.06319485078994, 1.06764705882353, 1.08654416123296, 1.09310113864702,      1.06369008535785, 1.13811922753988, 1.12487100103199, 1.14294330518697,      1.15353181552831, 1.14426229508197, 1.1380042462845, 1.16727806309611,      1.09280544912729, 1.10660426417057, 1.13093858632677, 1.12244897959184,      1.09134045077106, 1.10821382007823, 1.09921875, 1.12583967756382,      1.0998268897865, 1.10657894736842, 1.12752114508783, 1.08413001912046,      1.14484272128749, 1.0859167404783, 1.09041501976285, 1.0887537993921,      1.05695364238411, 1.04765146358067, 1.04174820613177, 1.05854800936768,      1.04042904290429, 1.07479752262982, 1.07179197286603, 1.05997624703088,      1.06460369163952, 1.07920193470375, 1.081811541271, 1.08351810790835,      1.0703933747412, 1.07135523613963, 1.0532319391635, 1.05964730290456,      1.07206703910615, 1.07498383968972, 1.05938566552901, 1.08185840707965,      1.06121372031662, 1.05117647058824, 1.0734494015234, 1.05576208178439,      1.08180628272251, 1.06072555205047, 1.09534671532847, 1.08269794721408,      1.0863453815261, 1.07660577489688, 1.11460957178841, 1.09818731117825,      1.06873428331936, 1.08247925817472, 1.06818181818182, 1.09494725152693,      1.11903160726295, 1.10917361637604, 1.09464701318852, 1.10445468509985,      1.08333333333333, 1.06683804627249, 1.06380575945793, 1.07498766650222,      1.07160253287871, 1.07565588773642, 1.05174927113703, 1.07279344858963,      1.06560283687943, 1.06727037516171, 1.05085682697623, 1.06547285954113,      1.08014705882353, 1.0575296108291, 1.05748725081131, 1.04852071005917,      1.05421686746988, 1.05314846909301, 1.0538885486834, 1.04618937644342,      1.04105344694036, 1.06053604436229, 1.06058788242352, 1.04755700325733,      1.04994511525796, 1.05405405405405, 1.06622516556291, 1.07163323782235,      1.07538994800693, 1.06018957345972, 1.07800751879699, 1.07815198618307,      1.07247665629169, 1.07490217998882, 1.06998939554613, 1.05968331303289,      1.05139565795304, 1.07414104882459, 1.09087423312883, 1.06742556917688,      1.06096361848574, 1.07464929859719, 1.09754281459419, 1.10085400569337,      1.08974358974359, 1.09106168694922, 1.09333865177503, 1.08897569444444,      1.07627737226277, 1.14392723381487, 1.06422018348624, 1.07022471910112,      1.07848837209302, 1.06617647058824, 1.0828331332533, 1.08257858284497,      1.07761904761905, 1.06547619047619, 1.07017543859649, 1.06287069988138,      1.09431751611013, 1.09341500765697, 1.06916019760056, 1.06135831381733,      1.06491326245104, 1.06208955223881, 1.06825232678387, 1.06939409632315,      1.05837912087912)      * * *                 x<-ts(x, frequency=7)      When I try to:                 plot(forecast(ets(x),h=60))       plot(forecast(x,h=60))      ![plot](http://i.stack.imgur.com/qubPZ.png) I get the same results. Maybe someone could explain, why exponential smoothing in this case makes no difference? Also I have tryed to use                > plot(forecast(auto.arima(x),h=60))           Warning message:     In auto.arima(x) :       Unable to fit final model using maximum likelihood. AIC value approximated


title: Post hoc tests from Anova() from car-package

text: I have a small sample unbalanced RCBD, where I have reason to suspect unequal variances, hence I wish to fit a ANOVA with type-3 SS with heterscedasticity consistent esimators. The simplest way to do this is to utilize Anova() function from the package CAR, however, I cannot for the life of me understand how to run post hoc tests after using this function. Multcomp() doesnt let me use the object created in Anova(), neither will it let be use the similar (although without hce) drop1() objects. Can anyone help me with this? I know I can do this in SAS but I much rather learn to do it in R.


title: How to calculate loadings and use them to regress in SPSS

text: I am trying to build a model with 15 independent variables. However, because the parameters have been showing collinearity, I have been advised to use the Principal Component regression method. However, I am not sure how to calculate the principal components that can be used an inputs in regression. I am currently using SPSS and again am not sure if it is viable to use the factor loadings calculated when performing a PCA as principal components (new variables) for regression. I will be glad to have your suggestions.


title: Cross-correlation significance in SPSS

text: How do I test the significance of cross-correlation values in SPSS? It seems there is no option for this.


title: Segmenting an interval sensibly

text: Is there a canonical/recommended approach to or algorithm for splitting up an interval with the intent of minimizing the number of segments while keeping a high accuracy? It is essentially an optimization problem and one possible solution would be to apply classification trees as I describe below. But this is a typical problem faced in Machine Learning when it seems sensible to feed a feature as a category instead of as a float. So I would like to learn something about what is a good way to tackle this and whether there is a name for this type of problem. Especially because the case gets trickier when more classes come into play and the objective function is not a count but a float resulting from a function on other features. ![enter image description here](http://i.stack.imgur.com/aqWFL.png)               library(party)     D <- expand.grid(type = c("A","B"), val = 1:8)     D$num <- c(3,4,2,5,7,4,8,3,3,4,9,2,1,8,2,4)     D <- D[rep(1:nrow(D),D[,"num"]),-3]     plot(ctree(type~val,D,controls=ctree_control(mincriterion=.41)))      ![enter image description here](http://i.stack.imgur.com/kx6U3.png)


title: Forecasting hourly time series with daily, weekly & annual periodicity

text: **Major edit: I would like to say big thanks to Dave & Nick so far for their responses. The good news is that I got the loop to work (principle borrowed from Prof. Hydnman's post on batch forecasting). To consolidate the outstanding queries:** **a) How do I increase the maximum number of iterations for auto.arima - it seems that with a large number of exogenous variables auto.arima is hitting maximum iterations before converging on a final model. Please correct me if I'm misunderstanding this.** **b) One answer, from Nick, highlights that my predictions for hourly intervals are derived only from those hourly intervals and are not influenced by occurrences earlier on in the day. My instincts, from dealing with this data, tell me that this shouldn't often cause a significant issue but I am open to suggestions as to how to deal with this.** **c) Dave has pointed out that I require a much more sophisticated approach to identifying lead/lag times surrounding my predictor variables. Does anyone have any experience with a programmatic approach to this in R? I of course expect there shall be limitations but I would like to take this project as far as I can, and I don't doubt that this must be of use to others here as well.** **d) New query but fully related to the task at hand - does auto.arima consider the regressors when selecting orders?** I am trying to forecast visits to a store. I require the ability to account for moving holidays, leap years and sporadic events (essentially outliers); on this basis I gather that ARIMAX is my best bet, using exogenous variables to try and model the multiple seasonality as well as the aforementioned factors. The data is recorded 24 hours at hourly intervals. This is proving to be problematic because of the amount of zeros in my data, especially at times of the day which see very low volumes of visits, sometimes none at all when the store has just opened. Also, the opening hours are relatively erratic. Also, computational time is huge when forecasting as one complete time series with 3 years+ of historical data. I figured that it would make it faster by computing each hour of the day as separate time series, and when testing this at busier hours of the day seems to yield higher accuracy but again proves to become a problem with early/later hours that don't consistently receive visits. I believe the process would benefit from using auto.arima but it doesn't seem to be able to converge on a model before reaching the maximum number of iterations (hence using a manual fit and the maxit clause). I have tried to handle 'missing' data by creating an exogenous variable for when visits = 0. Again, this works great for busier times of day when the only time there are no visits is when the store is closed for the day; in these instances, the exogenous variable successfully seems to handle this for forecasting forward and not including the effect of the the day previously being closed. However, I'm not sure how to use this principle in regards to predicting the quieter hours where the store is open but doesn't always receive visits. With the help of the post by Professor Hyndman about batch forecasting in R, I am trying to set up a loop to forecasting the 24 series but it doesn't seem to want to forecast for 1 pm onwards and can't figure out why. I get "Error in optim(init[mask], armafn, method = optim.method, hessian = TRUE, : non-finite finite-difference value [1]" but as all the series are of equal length and I am essentially using the same matrix, I don't understand why this is happening. This means the matrix is not of full rank, no? How can I avoid this in this approach? https://www.dropbox.com/s/26ov3xp4ayig4ws/Data.zip               date()          #Read input files     INPUT <- read.csv("Input.csv")     XREGFDATA <- read.csv("xreg.csv")          #Subset time series data from the input file     TS <- ts(INPUT[,2:25], f=7)               fcast <- matrix(0, nrow=nrow(XREGFDATA),ncol=ncol(TS))          #Create matrix of exogenous variables for forecasting.     xregf <- (cbind(Weekday=model.matrix(~as.factor(XREGFDATA$WEEKDAY)),                         Month=model.matrix(~as.factor(XREGFDATA$MONTH)),                     Week=model.matrix(~as.factor(XREGFDATA$WEEK)),                         Nodata=XREGFDATA$NoData,                     NewYearsDay=XREGFDATA$NewYearsDay,                         GoodFriday=XREGFDATA$GoodFriday,                     EasterWeekend=XREGFDATA$EasterWeekend,                         EasterMonday=XREGFDATA$EasterMonday,                     MayDay=XREGFDATA$MayDay,                         SpringBH=XREGFDATA$SpringBH,                     SummerBH=XREGFDATA$SummerBH,                         Christmas=XREGFDATA$Christmas,                     BoxingDay=XREGFDATA$BoxingDay))     #Remove intercepts     xregf <- xregf[,c(-1,-8,-20)]          NoFcast <- 0          for(i in 1:24) {            if(max(INPUT[,i+1])>0) {            #The exogenous variables used to fit are the same for all series except for the       #'Nodata' variable. This is to handle missing data for each series        xreg <- (cbind(Weekday=model.matrix(~as.factor(INPUT$WEEKDAY)),                          Month=model.matrix(~as.factor(INPUT$MONTH)),                      Week=model.matrix(~as.factor(INPUT$WEEK)),                          Nodata=ifelse(INPUT[,i+1] < 1,1,0),                          NewYearsDay=INPUT$NewYearsDay,                      GoodFriday=INPUT$GoodFriday,                          EasterWeekend=INPUT$EasterWeekend,                      EasterMonday=INPUT$EasterMonday,                          MayDay=INPUT$MayDay,                      SpringBH=INPUT$SpringBH,                          SummerBH=INPUT$SummerBH,                      Christmas=INPUT$Christmas,                          BoxingDay=INPUT$BoxingDay))       xreg <- xreg[,c(-1,-8,-20)]            ARIMAXfit <- Arima(TS[,i],                           order=c(0,1,8), seasonal=c(0,1,0),                          include.drift=TRUE,                          xreg=xreg,                          lambda=BoxCox.lambda(TS[,i])                          ,optim.control = list(maxit=1500), method="ML")                   fcast[,i] <- forecast(ARIMAXfit, xreg=xregf)$mean           } else{       NoFcast <- NoFcast +1      }     }          #Save the forecasts to .csv     write(t(fcast),file="fcasts.csv",sep=",",ncol=ncol(fcast))               date()      I would fully appreciate constructive criticism of the way I am going about this and any help towards getting this script working. I am aware that there is other software available but I am strictly limited to the use of R and/or SPSS here... Also, I am very new to these forums - I have tried to deliver as full an explanation as possible, demonstrate the prior research I have done and also provide a reproducible example; I hope this is sufficient but please let me know if there is anything else I can provide to improve upon my post. EDIT: Nick suggested that I use daily totals first. I should add that I have tested this and the exogenous variables do produce forecasts that capture daily, weekly & annual seasonality. This was one of the other reasons I thought to forecast each hour as a separate series though, as Nick also mentioned, my forecast for 4pm on any given day will not be influenced by previous hours in the day. EDIT: 09/08/13, the problem with the loop was simply to do with the original orders I had used for testing. I should have spotted this sooner and places more urgency on trying to auto.arima to work with this data - see point a) & d) above.


title: Issues with ordinary kriging

text: I was following this wiki article related to ordinary kriging ![enter image description here](http://i.stack.imgur.com/UFQMB.png) Now my covariance matrix looks like this, for 4 variables               1   0.740818220681718   0.548811636094027   0.406569659740599     0.740818220681718   1   0.740818220681718   0.548811636094027     0.548811636094027   0.740818220681718   1   0.740818220681718     0.406569659740599   0.548811636094027   0.740818220681718   1      Well the relation between semvariogram and variogram is given by $\gamma(h)/(C0) = 1 - C(h)/C(0)$ So, I calculated the $\gamma(h)$ as well. Now when I try to calculate the weights as               A = 1.0000    0.7408    0.5488    1.0000         0.7408    1.0000    0.7408    1.0000         0.5488    0.7408    1.0000    1.0000         1.0000    1.0000    1.0000         0           B =  0.4066         0.5488         0.7408         1.0000      I am taking the fourth variable as missing                [W;mu] = inv(A)*B =  0.1148                           0.0297                           0.8555                          -0.1997      The above was by using covariance. Now using semi variance I had               A = 0         0.2592    0.4512    1.0000         0.2592         0    0.2592    1.0000         0.4512    0.2592         0    1.0000         1.0000    1.0000    1.0000         0          B = 0.5934         0.4512         0.2592         1.0000               inv(A)*B =  0.1148                 0.0297                 0.8555                 0.1997      As you can see the last terms are not equal. When according to the derivation they are equated or said to be equal. Any clarifications?


title: Denormalize value after prediction

text: I don't know how to denormalize (0-1 normalized) data after prediction. I have 1 output and several input values. It's clear for 1 input I must use min and max value previously used for normalization. But what should I do if for normalizing, each input was normalized separately, using its own min and max value? Does anyone have any ideas?


title: Kernel PCA with an SVD algo

text: Suppose that I have a great algo for calculating the SVD and I want to do Kernel PCA. It is possible to first apply the Kernel function to my data and then run the SVD algo on the transformed data?


title: Regression with variables in different scale (how to encode properly and what regression to use?)

text: I'm working here on my research proposal for my alma mattera. I'm in accountancy science field, however, I need to deal with some statistics on this occassion. Would be very grateful if you could help me. I have an idea to built this kind of regression model: Y= x1 + x2 + x3 The research question is: Do audit prices, non audit services provided for audit client and auditor's tenure has an impact on audit quality? Observation - financial reporting of particular company. Let's say, we have 40 companies and their audited financial reporting. Y (proxy for audit quality) - count of errors in financial reporting. I guess it would suitable to group the errors found in each observation. For instance, range from 1 to 5 (where 1 - from 0 to 2 errors, 2 - from 2 to 4... and so on.). Or it could be left in relative scale - just how many errors (mistakes in financial reporting) I will found. x1 - price of audit services of financial reporting. It could be left in model in a relative scale, or coded as 1 - price above the average price in the sample, 0 - price bellow the average price in the sample x2 - non-audit services provided for the audit client. 1 - there were other services provided, 0 - there where no other services provided. x3 - the tenure of auditors. Like: 1 - the same auditors are auditing the financial reporting 3 or more years in the row, 0 - the same auditors are not auditing the financial reporting more or equal 3 years. Or it could be put in the model in relative scale: 2, 3, 4, 5, ...8 years in a row. That's my assumptions how to encode the data. Could you express your opinion how the data should be encoded actually, what kind of regression should I use here and what statistical tests would be suitable for evaluating the results? Thank you very much


title: scikit-learn SkewedChi2Sampler - meaning of skewedness parameter

text: I am trying to understand the meaning of the "skewedness" parameter for scikit-learn's SkewedChi2Sampler and figure out how this value affects the output of the sampler. I have looked at the docs here and all it says is that this should be set using cross-validation but I'm hoping someone can provide some intuition on how this parameter should be set so that I can interpret whether the parameters returned by cross-validation makes sense. I have also checked the referenced paper "Random Fourier Approximations for Skewed Multiplicative Histogram Kernels" but I'm not sure which parameter in the paper this "skewedness" parameter corresponds to. I'm guessing it might be the parameter "c" in the chi2 kernel. Can anyone tell me whether I'm correct?


title: Distribution of $\bar{X^2} $ when $X\sim N \left( \theta, \sigma^2 \right) $

text: How can I derive the distribution of $$\bar{X^2}\quad \text{when}\quad X\sim N \left( \theta, \sigma^2 \right) $$ The context of this question is an exercise requiring me to show that $\bar{X^2}- \frac{\sigma^2}{n} $ is an unbiased estimator of $\theta^2$ and afterwards find its efficiency. Thus I must find its variance as well. Obviously $\bar{X} \sim N \left(\theta ,\frac{\sigma^2}{n} \right) $ but when it comes to its square I would not know. I suppose that my only option is to derive the distribution of $\bar{X^2}$ using the CDF technique and repeatedly use change of variables. I tried that but it leads to a mess. I cannot discern the resulting distribution either Is there perhaps a neater way to derive the distribution of $\bar{X^2}$? Thanks.


title: How does one generate the table mapping t-test values to p values?

text: In the dark ages, we would map the results of a Student's t-test to a null hypothesis probability _p_ by looking up _T_ and degrees of freedom in a table to get an approximate result. What is the mathematical algorithm that generates that table? _ie_ , how can I write a function to generate a precise _p_ given an arbitrary _T_ and _df_? The reason I ask is that I'm writing a piece of embedded software that continually monitors hundreds of populations with hundreds of samples each, and raises an alert if successive snapshots of a given population come to differ significantly. Currently it uses a crude _z_ -score comparison, but it would be nice to use a more valid test.


title: Calculating the transfer entropy in R

text: The transfer entropy, from information theory, is an effective way to measure the one-way information dependence between two variables. A nice high-level summary is here: http://lizier.me/joseph/presentations/20060503-Schreiber- MeasuringInfoTransfer.pdf I see that there is a package for entropy and _mutual_ information estimation (http://strimmerlab.org/software/entropy/), but not the one-way transfer metric. What is an efficient way to calculate this in R? Perhaps I can use a chart output or metric from the mutual information package as a startpoint.


title: Searching for geometric network dataset

text: I'm looking for a large dataset (~20k nodes,~100k edges) for a geometric network, like utility networks, etc (See) Ideally, at a minimum I would like a dataset that includes identified 'nodes' with geographic coordinates and a connectivity arrangement between the nodes. As a bonus, some timestamped quantities that are measured between the nodes would be great. I've been looking in all the regular places (data.gov, data4nr.net, d8taplex.com, data.gov.uk) but I can't seem to find something that could meet all this criteria. I know I am asking too much but would like to know if anyone is aware of an intermediate dataset. Thanks in advance.


title: What is the name for a "typical" range of values?

text: First off, let me say that I'm extremely poorly versed in statistics, and this is a question purely about terminology. I have a distribution for some quantity (height of person, say) and the most likely 95% of outcomes are within some range $h_0$ to $h_1$. I want to know what to call (ideally a 2 or 3 word phrase) this range. It is _not_ the case that the distribution is normal, though we can safely assume that it is unimodal. **Clarification:** Let's say I have a coin, which I think is biased for heads with $p_h=0.6$. Now I toss it $n$ times and I get $a$ heads. I want to write a sentence which basically states that the result is in some 95% "likelihood region", (something like "middle 95 percentile" seems a little clunky, if not ambiguous --- I'm not taking a region about some median!). Now, it might be that my $p_h$ is actually a prediction based on some data, and I have some uncertainties in it, but I _don't_ want to include those uncertainties! I want a term that refers to the fact that because I didn't toss the coin an infinite number of times, there is some variance in the ratio $a/n$.


title: Hierarchical modeling in R

text: Here i am again. From this script:               model1 <- zelig(Decision ~ as.factor(age) + as.factor(income) + as.factor(town.type) +                     tag(town.type|town), model="logit.mixed", data=Fish)      Does it mean that `town.type` is treated both as fixed effects and random effects and `town` treated as random effects nested in `town.type`? I am currently analyzing the effects of `age` and `income` to the decision of fishers to exit the fishery. In addition i added `town.type` (whether the respondent comes from developed (coded 2) or less developed (coded 1) towns ) as fixed effects. We surveyed many towns each of which we categorized into one of two town types. I'm sorry, i am new to R. Hope you could help me out. Thanks a lot.


title: Custom Neural network in matlab

text: I am trying to make a custom Neural network structure using 'network' command. I am a little confused.Can we change the connections between individual neurons?Like it is possible in weka?


title: Learning theory for search data?

text: Has learning theory ever been applied in practice for search log data? If so, what are some findings about generalization/learnability from this data? I'm interested in generalization about an individual user's search behavior from a time slice or generalization about a total population of users, given a sample. Paper references are a plus!


title: How to interpret coefficents from a polynomial model fit?

text: I'm trying to create a second order polynomial fit to some data I have. Let's say I plot this fit with `ggplot()`:               ggplot(data, aes(foo, bar)) + geom_point() +             geom_smooth(method = "lm", formula = y ~ poly(x, 2))      I get: ![](http://i.stack.imgur.com/t5uO4.png) So, a second order fit works quite well. I calculate it with R:               summary(lm(data$bar ~ poly(data$foo, 2)))      And I get:               lm(formula = data$bar ~ poly(data$foo, 2))          ...          Coefficients:                         Estimate Std. Error t value Pr(>|t|)         (Intercept)         3.268162   0.008282 394.623   <2e-16 ***     poly(data$foo, 2)1 -0.122391   0.096225  -1.272    0.206     poly(data$foo, 2)2  1.575391   0.096225  16.372   <2e-16 ***          ....      Now, I would assume the formula for my fit is: $$ \text{bar} = 3.268 - 0.122 \cdot \text{foo} + 1.575 \cdot \text{foo}^2 $$ But that just gives me the wrong values. For example, with $\text{foo}$ being 3 I would expect $\text{bar}$ to become something around 3.15. However, inserting into above formula I get: $$ \text{bar} = 3.268 - 0.122 \cdot 3 + 1.575 \cdot 3^2 = 17.077 $$ What gives? Am I incorrectly interpreting the coefficients of the model?


title: How can I smooth a set of discrete data points for the purpose of schedule planning?

text: Disclaimer: I do not have a background in statistics or the math behind filtering, save one long-time-ago college course. I have a well defined problem space. I am calculating hourly staffing requirements for retail stores. The raw values are based on a matrix of the full-time employees required for a given historical transaction volume. The calculation process generates a set of ordered pairs ( hour, people ) for a given store and day. The output of this process can produce single hour periods where the number of people required are much higher or much lower than the neighboring range, but the concern is only the lower values. So if the calculated data (hour, people) for a given store and day were to be:               8, 2   -> I need two people to work from 8a - 9a     9, 3   -> I need three people to work from 9a - 10a, etc.     10, 3     11, 2     12, 3     13, 3     14, 1     15, 3     16, 3     ...      I need a filtering approach to: turn (11,2) into (11,3) and turn (14,1) into (14,3) I would never want an outputted value to result in a lower value than the source data. **EDIT - Additional Information ** The practical need for this data and the smoothing process is for a application user interface. The store managers schedule their staff for a given week using this data as a guide. The raw data represents the minimum staff coverage for each hour. When store manager builds their schedule, the application flags time periods where insufficient coverage is detected, i.e. below the calculated minimum staff coverage. Schedules are built up based on shifts, which means it is impractical to have someone work 2 hours, take an hour off, then work another 2 hours. We would want that person to work the entire 5 hours. We want to modify the raw data in such a way as to avoid having these coverage spikes as shown in the example above.


title: Two questions about Bonferroni adjustment

text: I have two questions about Bonferroni adjustments: 1). Can one use the Bonferroni method to compare independent groups? The reason why I ask this is it seems that many examples I've encountered discuss the Bonferroni method in the context of comparing dependent groups - for example, multiple comparisons after repeated measures ANOVA. 2). I created a set of simulated data (see code below for the reproducible dataset).               set.seed(123)     data<-data.frame(x=rep(letters[1:4], each=5), y=sort(rlnorm(20)))      Then, I used `pairwise.t.test()` and set `p.adj="bonf"` (see below) to test pairwise comparisons.               pairwise.t.test(x=data$y, g=data$x, p.adj="bonf") #see results below:          #  data:  data$y and data$x      #    a       b       c           #  b 1.00000 -       -           #  c 0.38945 1.00000 -           #  d 8.3e-06 3.5e-05 0.00031          # P value adjustment method: bonferroni       However, these results are different from the results obtained by doing pairwise t-tests using `t.test()` and then adjusting for the p-values (see below)               t.test(y~x, data[data$x=="a" | data$x=="b",])$p.value*6     t.test(y~x, data[data$x=="a" | data$x=="c",])$p.value*6     t.test(y~x, data[data$x=="a" | data$x=="d",])$p.value*6     t.test(y~x, data[data$x=="b" | data$x=="c",])$p.value*6     t.test(y~x, data[data$x=="b" | data$x=="d",])$p.value*6     t.test(y~x, data[data$x=="c" | data$x=="d",])$p.value*6      The results are below:               # a vs. b = 0.0788128848                # a vs. c = 0.0001770066                # a vs. d = 0.0324680659     # b vs. c = 0.0137812904     # b vs. d = 0.0488036762     # c vs. d = 0.0970799045      These adjusted p-values are rather different from the ones obtained from individual t-tests. So I wonder why there are such big differences.


title: Why does asymptotic distribution of LRT not depend on specific functions used to express null set?

text: For a likelihood-ratio test (LRT) statistic, $\Lambda(x)$, the asymptotic distribution of the statistic $-2 \log \Lambda(x)$ is a chi-squared distribution, which occurs only under certain regularity conditions. Here, we have two hypotheses for the unknown vector $\theta$:   1. $H_0: \theta \in \omega$ versus   2. $H_1: \theta \in \Omega \cap \omega^C$. Here, $\Omega$ is the full parameter space, and $\omega$ is the restricted subspace for the null hypothesis. The subspace $\omega \subset \Omega$ is defined in terms of $q$ independent constraints of the form $g_1(\theta) = a_1, \ldots, g_q(\theta) = a_q$, where $0 < q \le p$ and $a_1, a_2, \ldots, a_q$ are constants. The functions $g_1, \ldots, g_q$ must be continuously differentiable. Why does the asymptotic distribution of $-2 \log \Lambda(x)$ not depend on what the functions $g_1, \ldots, g_q$ actually are as long as they are smooth and independent?


title: Checking assumption of dependence in a generalised linear model

text: I'm wondering if I've violated the assumption of independence in a generalised linear model. The data I've used to make to make the model can be downloaded with this code:               library(repmis)     goose_feeding <- repmis::source_DropboxData("goose_feeding.csv", "gz6v6hvjytqn14n", sep = ",", header = TRUE)          goose_feeding$datetime <- as.character(goose_feeding$datetime)     goose_feeding$datetime <- dmy_hm(goose_feeding$datetime)        * pecking = proportion of a goose flock that are feeding    * water height = water height in metres   * ssp = subspecific identity of goose flock   * days = days since 27th September 2013 This is my model:               glm_pecking <-  glm(pecking ~ ssp*poly(water.height, 2) + ssp*poly(days, 2), family = quasibinomial, dataScanSampling_sub)      To check for dependence across time, I simply plotted pecking against date:               library(ggplot2)     ggplot(goose_feeding, aes(datetime, pecking)) + geom_point()      But this plot doesn't take into account that groups of observations were collected at the same water height or for same subspecies. How can I check for dependence across time in `goose_feeding`?


title: Measuring and analyzing sample complexity

text: I recently stumbled upon the concept of **sample complexity**, and was wondering if there are any texts, papers or tutorials that provide:   1. An introduction to the concept (rigorous or informal)   2. An analysis of the sample complexity of established and popular classification methods or kernel methods.   3. Advice or information on how to measure it in practice. Any help with the topic would be greatly appreciated.


title: compare the values assumed by two trend lines in correspondence a certain x value

text: I'm comparing samples of remote sensing data points collected over two different areas, so basically I come up with a pair of data set and I run the same analysis on both of them. Comparing the linear models (t-test) I found that the slopes and intercepts are significantly different from each other. Nonetheless, since my independent variable ranges from 0 to 1 and since the two trend lines approaches each other as x increases, I would like to assess whether the values assumed by the trend lines for x=1 are significantly different.   * Do I have to use again the t-test, as I did to compare the two models? In this case are the degrees of freedom the same?   * Is it sufficient to consider the uncertainties derived by summing slopes and intercepts (and their relative standard errors) and check if the ranges overlap?   * Do I have to plot confidence intervals for the two distributions? Thanks in advance.


title: multiplication of 2 PDFs

text: If I multiply the two PDFs, does the variance of the result PDF becomes narrower than the two PDFs always? In other words, if I multiply likelihood and prior to get the posterior, is the variance of the posterior narrower than that of the likelihood and prior? Thanks


title: How to assess the similarity of two histograms?

text: Given two histograms, how do we assess whether they are similar or not? Is it sufficient to simply look at the two histograms? The simple one to one mapping has the problem that if a histogram is slightly different and slightly shifted then we'll not get the desired result. Any suggestions?


title: Comparing Correlations Between Two Groups with Replication

text: I would like to compare whether or not the correlation between 2 variables are the same in 2 experimental groups, each with 6 replicates. The two variables that I would like to compare are non-normally distributed. For each replicate, I have ~100-150 measurements of both variables and a total sample size of around 1,550. The popular test for this sort of comparison seems to be Fisher's r to z transformation, but this obviously would not take into account the structure of my data, given the multiple replicates per experimental group. Can anyone suggest an approach that I might use? I had considered the possibility of creating a generalized linear model that would look something like the following: Variable A ~ Variable B + Experimental Group + Experimental Group:Variable B + (1|Replicate ID) Would this be appropriate? My thinking is that this is probably not okay, given that I'm dealing with correlations here. Thanks very much for any advice.


title: How do I find the probability of a type II error?

text: I know that a Type II error is where H1 is true, but H0 is not rejected. ### Question How do I calculate the probability of a Type II error involving a normal distribution, where the standard deviation is known?


title: How to calculate OR and confidence interval for two binary explanatory variables with interaction

text: I have a logistic regression model with several categorical explanatory variables and one interaction term (between two binary variables, named A and B). I know how to calculate the odds ratios for the different levels of A and B (for A=1, e.g., I need to add the coeff for A and coeff for A*B, then exponentiate), but how do I get a confidence interval for this OR? I need to do this in R please, this is the only statistical package I have access to.


title: How auto.arima works?

text: I found `auto.arima()` function in `forecast` packages. As far I have understand how it works, It should find the best model for the data, my question is: Is it possible to understand if the variance is not constant during all the period(timeseries) using that function? I mean... If it fit a model it should means that the variance is homoscedastics because i set `stationary=TRUE` Let me know, thank you!


title: How to apply unsupervised classification to spatial data

text: I am trying to learn how to apply unsupervised classification to spatial data. In near infra red satelite pictures; the ocean is dark and the forest is white. Each pixel in such an image is an observation. I guess the variables are [x-position, y-position, grayscale]? Which methods could I use to classify the pixels into ocean and forest? I want to ensure that the points in the classes forms spatially connected regions (e.g. no single trees in the ocean). Thanks in advance for any answers!


title: Best algorithm for classifying time series motor data

text: I am working on a machine control project. We can measure the motor's current during operation. Sample data from two motors performing an operation successfully is below. The red trace shows the current from one motor, the blue trace the current from another. I'd like to try and come up with an algorithm for identifying problems with machine behavior. Problems could be excessively high motor current, near zero motor current, current increasing at the end of the operation, a shorter time series than normal, anything in general which doesn't look like a typical operation below. Can anyone suggest a good algorithm for achieving this? The only one I'm familiar with is a neural network. I have put an Excel file of actual data at motor currents ![Motor currents - good operation](http://i.stack.imgur.com/RQka3.png) ![Motor currents - jam at end of operation](http://i.stack.imgur.com/tOIpI.png)


title: How do I avoid overlapping labels in an R plot?

text: I'm trying to label a pretty simple scatterplot in R. This is what I use:               plot(SI, TI)     text(SI, TI, Name, pos=4, cex=0.7)      The result is mediocre, as you can see (click to enlarge): ![enter image description here](http://i.stack.imgur.com/stAuk.png) I tried to compensate for this using the `textxy` function, but it's not better. Making the image itself larger doesn't work for the dense clusters. Is there any function or easy way to compensate for this and let R plot labels that _don't overlap_? * * * Here is a small subset of the data I have:               Name;SI;TI     01_BAD_talking_head;6.944714;4.421208     01_GOOD_talking_head;5.680141;4.864035     01_GOOD_talking_head_subtitles;7.170114;4.664205


title: Density related to sparseness measure

text: Are there any multi-variate continuous distributions whose probability distribution functions give high values for sparse vectors and low values for dense vectors, i. e. indicating the sparseness of the vector? I found one sparseness measure in [1] which is as follows for $x \in \mathbb{R}^n$: $$s(x) = \frac{\sqrt{n}-\frac{\sum |x_i|}{\sqrt{\sum x_i^2}}}{\sqrt{n}-1}$$ It gives 1 if the vector has only one non-zero element and 0 if all the elements are equal as given below (the image is from the referenced paper). ![Sparseness values](http://i.stack.imgur.com/4a3mS.png) This seems to be a good candidate for a _potential_ instead of a _density_ and it helps in my case since I compute the likelihood on the given data. I wonder if there exists any known distributions which can be used in this case. And how can I generate random vectors belonging to that distribution? [1] P. O. Hoyer, _Non-negative matrix factorization with sparseness constraints_ , Journal of Machine Learning Research, vol. 5, pp. 1457–1469, 2004.


title: How best to reduce dimensionality of a dataset composed of events and trials?

text: I'm trying to reduce dataset dimensionality. PCA is a good metric but that gives me new dataset. My goal is to determine from number of events (e.g. 60) and number of trials (e.g. 6) which events are more relevant. For example:   * 1st, 3rd, 21st, 45th ... (N total) events are good enough to approximate behavior of dataset.   That will allow me to discard 60-N events, and to deal with only N. For now, I'm calculating covariance matrix, and take events for which correlation is smaller than some threshold.   Is there some official metric or math function for this? Thanks.


title: How to handle data imbalance in classification?

text: I am working on a text classification problem. My data is highly imbalanced. For example, one category has 700 documents while the other has 30. I have around 30 categories. I tried different classifiers and the performance is consistently poor. What is the best way to tackle this issue? Thanks


title: Getting rweibull to output n observations in (1,52) given specific shape, scale parameters

text: I was told the following will generate 52 observations in range (1,52), with a shape parameter of .5 and a scale parameter of 2. Yet, there is a 67 so that is perhaps wrong?   What I need is N observations in range (1,52) with the observations conforming to a Weibull distribution with specific parameters.               > print(rweibull(52,.5,2))           [1]  0.406647225  0.052710518  2.158733355  0.500174364  2.074537224          [6]  0.378695917  0.396309715  0.131268753  4.182770394  0.293154663         [11]  2.045214882  0.016280681  0.268003469 17.612960070  2.375416850         [16]  6.482297632  0.151320479  0.214247990  2.613026032  3.397453127         [21]  7.357103077  2.054468323  3.090147233  2.594981565 21.316501285         [26] 35.364512001  2.272848255  1.186967955  4.178822565  0.006659194         [31]  6.006123713 ***67.918717193***  0.051010700  0.706372359  2.419090546         [36]  1.858408184  9.545857635  2.182673324  0.892914986  0.244567333         [41]  0.065993717  0.225678795 27.655436513  2.965634299  1.294741473         [46]  0.069503896  1.266168448  3.949787345  3.391448817  1.663533337         [51]  0.604392557  1.232525921      Any help is greatly appreciated! **Edit** : I apparently have solved my own question. The following will give me the expected # of failures in the Tth month, given k as the shape, and 100 initial # of parts:               ceiling((pweibull(t,shape=k,scale=12)-pweibull(t-1,shape=k,scale=12))*100)


title: Sparse parameters when computing AIC, BIC, etc

text: I'm designing large-scale, regularized logistic regression models with lots of sparse, binarized features. e.g. isUS, isFR, etc. As a result, a lot of the weights in the model are zero. I'm wondering how I should compute the "number of parameters" in model- selection criteria like AIC, BIC, etc. Should I only count the number of non- zero weights or all of the weights? e.g. If there are 10 possible countries, but only 6 have non-zero weights, is the number of parameters 6 or 10?


title: Standard Errors in Winsteps: ERROR versus MODLSE

text: I am running anchored analyses in Winsteps on four data sets (a full data set as well as the data degraded by an additional 20%, 50%, and 70%). I noticed that the standard errors are labeled ERROR for the full data and the data that is degraded an additional 50%, but are labeled MODLSE for the data that is degraded an additional 20% and 70%. This standard error relabeling occurs in both the item file (IFILE) and person file (PFILE) outputs. Does anyone know if ERROR and MODLSE are the same? Why the difference? Thanks!


title: p-Value of mean of sample means

text: I suspect this question is related to: When combining p values, why not just averaging? Let me explain my problem:   * Suppose you have a set of sample means and associated p-values, standard errors and degrees of freedom. Let’s denote them as: $\mu_i$, $p_i$, $\text{SE}_i$, $df_i$ where $i=\{1,n\}$   * It’s straight forward to figure out the mean of sample means. But what is the resultant standard error and p-value of the mean of sample means? I suspect p-values multiply in some manner…


title: Lognormal with negative values

text: I have some logged increments from time series data and wanted to fit a lognormal distribution, but obviously some are negative. How can I do this?


title: How to split nodes in regression trees

text: I am looking for a comparison of different regression tree node splitting approaches within the random forest framework. I am looking at the trade-off between ensemble accuracy/reliability (holding forest size constant) and computational complexity of the split since I deal with large datasets. The standard approach is to minimise $\sum_{i \in Region_1} (y_i - mean(y_i)_{Region1})^2$ +$\sum_{i \in Region_2} (y_i - mean(y_i)_{Region2})^2$. But another approach I have seen in a PhD thesis which drastically reduces the computational complexity is to look at the average of $(\sum_{i \in Region_1} y_i)^2/(Cardinality of Region_1)$ versus the same for $Region_2$; this allows me to use a running sum of $y_i$.


title: Whether to log transform variable when untransformed variable has positive skew and transformed has negative skew with additional missing data?

text: I have performed a log transformation on my skewed data, however on my DV it went from positive skew to negative skew after the (log) transformation, further data was missing from my DV after the transformation. Please help


title: A strange coefficient in a negative binomial model

text: I already asked a first question here fortunately or unfortunately the answers support my understanding of the general problem. Too bad this doesn't solve the problem. So, here is the thing: I am fitting a random effects negative binomial model. The dependent variable counts monthly events in certain regions. From descriptive analysis of the data I know that one region sees way more events than the others (e.g. the total number of counts in this region is around 2000 while all other regions stay below 400). So I included a dummy variable for this region in my model. It is coded 1 for the region and 0 for the rest. Surprisingly (at least to me) what I get for the dummy variable is a coefficient of -0.3365 or an incidence-rate ratio of .7 . To my understanding this means that if the dummy variable equals 1 the expected number of counts should decrease. To me this is counter intuitive. Does anybody have a possible explanation for this observation (other than that there's something horribly wrong with my model)? Or an idea how I could further investigate this problem? So far I came up with this idea: Most regions only have a count of 0 to 4 each month. For this special region however, the highest observed count is 79. Still, more than 50% of all observed counts for this region are below 9 and there are at least 40 possible values between 9 and 79 that are not observed at all. Could it be that the model expects observations for all numbers between 0 and 79 and that because so many are not observed I get an incidence- rate ratio below 1?   * If I recode all counts higher than 8 as 8, I suddenly get a incidence-rate ratio of 5.2


title: Confusion related to derivation of conditional distribution

text: I was reading this paper http://stat.duke.edu/phd- program/gcc/resources/ResourcesDocuments/CodeExplanation.pdf related to Gibbs sampling. Suppose we have n iid samples $x_i$ from $N(\mu,\sigma^2)$ Prior probabilities for $\mu$ ~ $N(m,s^2)$ $\sigma^{-2}$ ~ $Gamma(a,b)$ Then ![enter image description here](http://i.stack.imgur.com/h7kbY.png) Can anyone help me with how the conditionals are derived. I tried for $p(\mu|\sigma^2,x)$ which I got something like $\frac{N(\frac{\sigma^{-2}\mu + s^{-2}m}{\sigma^{-2} + s^{-2}}, \frac{\sigma^2s^2}{\sigma^2+s^2})}{p(x|\sigma^2)} $. Then I didn't get how to proceed and got stuck. This is what I tried then $p(\mu|\sigma^2,x) = \frac{N(\frac{\sigma^{-2}\mu + s^{-2}m}{\sigma^{-2} + s^{-2}}, \frac{\sigma^2s^2}{\sigma^2+s^2})}{\int N(\frac{\sigma^{-2}\mu + s^{-2}m}{\sigma^{-2} + s^{-2}}, \frac{\sigma^2s^2}{\sigma^2+s^2})d\mu}$. I don't think this will give me the required expression. Suggestions?


title: pLSA using tempered EM

text: In an article by Hofmann pdf, he proposes: > initialize $β$ to one, run until convergence, then rescale $β$ by a factor > $η<1$, run again until convergence, and iterate this until changing $β$ no > longer improves the result. I used the TEM algorithm on a dataset for my project, but the problem is that the algorithm ends up traped in several local minima and doesn't give me full convergence. **My question is: what is the optimal condition under which I can assume that the algorithm converged?**


title: Model selection, when candidates have been fit to two separate data sets independently...anything beyond multiplying likelihoods?

text: I have a series of candidate spatial capture-recapture models--but 2 or 3 of these models end up generating data vectors that R can't allocate when considering the entire dataset. I've split observations into two groups, and fit each candidate model twice, and have been able to piece together estimates and variance using the delta method...but I'm not sure how to go about model selection. After multiplying each model's log-likelihoods, differences in candidate AIC values that were approximately 1-2 when considering the entire data set are >1,000. Is it defensible to average each model's log-likelihoods or AIC weight for comparison? Thanks!


title: Does doing predictive regression need a smaller sample size than exploring risk factors with logistic regression?

text: Cardiologists have a tool called EUROScore used to adjust the risk associated with performing heart surgery. It comes into play when, for example, one surgeon is recognised as being more expert and so takes on the riskier patients: it adjusts his/her performance measures to take this into account. It also helps the surgeon when advising patients/family of the risk involved. We are developing a similar tool for use in Percutaneous Intervention (for example: use of stents). We started by using logistic regression but it was clear that the sample size required was greater than the sample size that we had available. In any case, I came round to thinking that propensity scoring might be the way forward. I would welcome any comments on that. More recently I have come round to thinking that predictive regression is the way forward. Again, I would welcome any comments. Assuming that predictive regression IS the way forward, am I correct in thinking that the sample size required to do this is much less than the sample size required for logistic regression? If so, is there a rule of thumb for determining the sample size required. Thanks in advance, Martin


title: Uniform Distribution

text: Suppose we want to bid on a piece of land when one other bidder is interested. The seller announced that the highest bid in excess of \$10,000 will be accepted. Assume that the competitor's bid $X$ is a random variable that is uniformly distributed between \$10,000 and \$15,000.   1. Suppose you bid $12,000. What is the probability yours will be accepted?   2. What amount should you bid to maximize the probability that you get the property. My approach to this problem is as follows: Opponent's mean bid = \$12,500 Opponent's Std Dev of bid = (\$15,000-\$10,000) / 6 = \$833.33 12000-12500 = -500 -500/833.33 = -0.6 Zscore = 0.27425 27.425% that the 12000 bid will be accepted?


title: How to model a multi-dimensional feature set for classification

text: I am new to statistical modelling and so please pardon if the question appears trivial. I have a set of multi-dimensional data ($T$) where each dimension represents features ($f_i$) obtained from a mammogram. For example $T_1=(f_{11}, f_{21} f_{31},\ldots, f_{n1}); T_2= (f_{12}, f_{22}, f_{32},\ldots, f_{n2});\ldots; T_M$ . Referring to the paper > Ge, Srinivasan & Krishnan (2002)   >  _Cardiac arrhythmia classification using autoregressive modeling_ ,   >  BioMedical Engineering OnLine, 1:5 (link), I am interested to know how to fit a polynomial or an AR model to begin with so that a general model may be obtained for classifying the feature sets into two classes; how many parameters are required, etc? I fail to understand how and why to use the coefficients as features when we have the feature set? How to employ the features for deriving an AR model? Say if I begin with AR model, how would I formulate so as to reach a general AR(p) model and so on? I would be glad for pointers such as research papers which have explained model-based fitting techniques and other pointers would be really helpful.


title: Justifying the distribution for the maximum likelihood estimator in a linear regression example

text: Data $(x_1, y_1), \dots, (x_n, y_n)$ is modelled with $x_i$ being non random and $y_i$ being observed values of $$Y_i = \alpha + \beta (x_i - \bar x) + \sigma \epsilon_i$$ with $\epsilon_i \sim N(0,1)$. I have calculated the MLE's of the parameters $\alpha, \beta$ and $\sigma^2$ with $\hat \alpha = \bar y$ and $$\hat \beta = \frac{\sum y_i(x_i - \bar x)}{\sum (x_i - \bar x)^2}$$ we now have a corresponding estimator for $\beta$ given by $$B = \frac{\sum Y_i(x_i - \bar x)}{\sum (x_i - \bar x)^2}$$ which apparently (according to my notes) is normally distributed with a mean of $\beta$ and a variance of $$\sigma^2 \over \sum(x_i - \bar x)^2$$ and I can't justify to myself why this is the case! Could someone please help explain this? Thanks!


title: Defining norm of a matrix of MFCC coefficients

text: **Scenario Construction:** I have a MFCC generator block which gets the speech samples from the user and generates a rectangular matrix say $A$ of the order $m \times n$, whose elements are the Cesptral Coefficients(MFCC). Now, suppose I maintain a database which are previously stored containing the user speech signal. Through an LPC filter, I generate the speech sample and then direct it to the MFCC generator block, with a constraint that I don't give the entire set of samples from the database for the filter to generate the speech signal. Rather, I give a part of the speech sample. Now this predicted speech signal will be now directed towards the MFCC generator block to generate the predicted signal's Cepstral coefficients which again turns out to be another rectangular matrix say $B$ of the same order $m \times n$. Then I use a matrix norm along with a heuristically chosen threshold to compare two matrices(find the error) and authenticate the user. If it fails, the input speech samples for the prediction is linearly increased and again the constraint is checked. **An insight into matrix A,B as defined before.** The rows of the matrix represent the number of coefficients to be generated per speech frame. The columns are the concatenation of coefficients of all the frames for the entire speech sample. A and B have the same setting. During MFCC generation we make use of a window of fixed size, operate on the samples under the window which yields coefficients for MFCC for that frame and then slide the window such that the slide steps are less than the window size, i.e.,: every successive window overlaps. **Question:** I have seen this http://stackoverflow.com/questions/6932096/matching-two- series-of-mfcc-coefficients link. I found it somewhat useful. Yet, I have a few concerns to raise in regards to the problem that I just defined. Even when an authenticated user speaks ( utters the exact word that is stored in the database) it is not necessary that the MFCC (positions of each element in the matrix) should be in exact match with the one that is generated during prediction.If both the rectangular matrix is converted into a vector, there may be a time delay between the samples. If so, the norms defined in the link I mentioned needn't work out even for an authenticated user. How do I fix it? Is there any other form of solving the above problem.


title: Kolmogorv Smirnov Test in R

text: I want to proof the "Relative Age Effect" of a football team. I have a list of birth dates of the team members (about 20 numbers between 1 and 365, the day of the year). I now want to use the KS-Test to test against a uniform distribution. this is how I did it, but it looks very wrong to me: > ks.test(dates,"punif") > > D = 1, p-value = 3.331e-16 alternative hypothesis: two-sided almost all players are born in the first half of the year, yet the p-value seems far too low in my opinion.


title: Finding degree of polynomial in regression analysis

text: I am working on a machine learning project where I am trying to fit a curve on data. Unfortunately the date has somewhat high feature vector. So, I can't really plot them on a 2D or 3D space to guess how the shape of the data looks like. So, other than hit and trial, does there exist a mathematical way to find degree of polynomial that might best fit my data. I mean I do know that I can look into a least square error for each degree and then choose the one with the minimum error, but then it will be a double optimization problem as the first optimization loop would consist in finding a set of weights for the curve that fits the data, while the second loop would be used to check the degree. Any suggestions?


title: Limit of c.d.f. of Poisson($n$) when $n$ goes to infinity

text: I'm trying to prove that $\lim_{n\rightarrow\infty}F_{X}(n)=1/2$ when $X\sim \text{Poisson}(n)$ without success. Could someone help me ?


title: Negative binomial GLM, the most complex model always has lowest AIC (all interaction terms)

text: I apologize for all those questions on modelling. It is the very first time that I try GLM and I am really lost even after reading a lot of papers. I have divided my covariates according to their theme e.g. topographic variables together. I was trying to first get the most parsimonious model for each set of variables. I am using the `step` function and I am running GLM **negative binomial** models. My continuous variables were scaled. Now, my problem is that whatever model I run the most complicated one is always the one with the smallest AIC. For example if I run : `glm(responsevar~Cov1*Cov2*Cov3)` the best model will be : `responsevar~Cov1+Cov2+Cov3+Cov1:Cov2+Cov1:Cov3+Cov2:Cov3+Cov1:Cov2:Cov3` My sample size is 2780. Would it be better to use BIC to select the models? Should I limit the interactions to what I think is coherent?


title: Generalized estimating equation and data format

text: I first summarize the analysis of some papers that I want to imitate. The regression looks like: `Y2(subjects) = Y2(friends) + Y1(subjects) + Y1(friends) + Z(subjects)` Where `subjects` refers to focal people and `friends` refers to those subjects connect to. Some one may be a subject in a relationship and a friend in another. For this question, just skip the validity of this model. The idea is to use the outcome variable of friends at time 2, `Y2(friends)`, to explain outcome variable of subjects at that time (`Y2(subjects)`), controlling for outcomes of both subjects and friends at time 1. Other covariates of subjects may be included, but it's not the main consideration. Following this model, I think the data should be in the **wide format** , which looks like: `SubjID FrID Y1Subj Y1Fr Y2Subj Y2Fr Z` The authors apply the **generalized estimating equation** in their analysis. However, I don't know how to perform GEE with the wide data format (let's say `PROC GENMOD` in `SAS`), if that's the case. Can you tell if I'm missing some thing (e.g., reconstruction of data in some way so that GEE can be used)?


title: Minimum Sample Size Required to Estimate the Probability $P(X \le c)$ for a Constant $c$ (Given a Confidence Level & Confidence Interval)

text: I have a large population of size $n$ from an unknown continuous random variable $X$, and I do not know the underlying distribution of $X$. Given a constant number $c$, I want to determine the minimum sample size I need to estimate the probability $P(X \le c)$ given a confidence level, $p_c$, and confidence interval, $I_c$ (I am not sure if we need them! ). How can I find the minimum sample size to estimate this probability? I have found the following discussion in Wikipedia which is independent of the number of population. I am not sure if it is a good way to determine sample size! ![enter image description here](http://i.stack.imgur.com/t09Kn.jpg) I have also found some methods to determine sample size for data to be analyzed by nonparametric tests.you don't have to make any assumption about the distribution of the values. That is why it is called nonparametric. Now I am confused if these nonparametric methods can be used to solve my problem or the method I found in Wikipedia is the correct way to solve my problem, or there exists a better solution. Thanks for your help.


title: Constrained Least Squares/Fixing Data

text: So I have maybe an unusual question. I have some simulation data that I use as input into a theoretical model. The problem is that the simulated data is noisy and this causes divergences when I try to input it into my theory. Well I figured out what conditions my simulated data must satisfy to ensure that I don't get the divergences. While they come close they don't satisfy the conditions to a high enough level of precision. I have 3 sets of data these three sets are related through three conditional equations. Now what I want to do is correct the noise of the simulations ever so slightly so as to ensure the three equations are satisfied. Working with my adviser we devised a plan to use a least squares fitting with known constraints used to sort of tweak the data points to ensure the constraints are met. So the idea is to fit the simulated data to itself with a vector of "tweaking" parameters. The idea is a data point from a respective simulation $$sim_x(r_i)$$ has its own tweaking point $$c_x(r_i)$$The idea would be $$\epsilon =\sum\limits_{i=1}^n((sim_x(r_i)+c_x(r_i))-sim_x(r_i))^2 $$So my least squares looks like this $$ min ||Ex - f|| $$ where $$ x = \{c_1(r_1),c_1(r_2),...c_1(r_n),c_2(r_1),c_2(r_2),...c_2(r_n),c_3(r_1),c_3(r_2),...c_3(r_n)\}_{3n} $$ is the tweaking vector to solve for and $$E=I_{3n}$$ and $$f=\{0,0,...0\}$$ So obviously just solving this is trivial as the solution that will minimize this is $$x = \{0,0,0...\}$$ So the thought is by adding the conditional statements it will make this non-trivial. $$Cx=d$$ Where $$C \in \mathbb{R}^{3n,3}$$ and $$d \in \mathbb{R}^3$$ I should note that all three constraint equations are linearly independent. Just to test things out I tried solving my constraint equations simply using Mathematica's Linear solve.(I think it just does solves it using the Pseudoinverse of C?) And it gives me my "tweaking" vector x with c values that do work and give me my desired result. However, I would much prefer to have a unique solution hence why I started with this whole least squares thing to begin with. I have tried various methods for doing least squares with constraints. Of note I tried the method outlined in "On the Method of Weighting for Equality- Constrained Least-Squares Problems" by Van Loan as well as Methods from "Solving Least Squares Problems" by Lawson and Hanson and just using Lagrange Multipliers. In all of the methods I end up back the exact same answer from just solving the constraint equations using Mathematica. This leads me to believe that likely I am not actually doing what I think I am doing which is solving the least squares at all. This method is a bit of a stretch so it is possible it won't work. Linear Algebra isn't my forte and though I am learning a great deal the reason it isn't working isn't obvious to me. Also, if anyone has any suggestions of other things I might try doing it would be greatly appreciated


title: Convert linear SVM answers to class probability?

text: How to do this? I use SVMlight that returns me some scores (which say how sure SVM is that something belongs to a class?) The questions is - can I do something to convert it to a % probability? Any formula, method that I could code for myself? Or is this impossible?


title: Suggestions for training in data visualization?

text: I'm a UI Designer/Developer working at a start-up with a LOT of data. I'm looking to do some training in Data Viz techniques. Anyone have suggestions for classes in NYC? Or online, if there are some good sites.


title: Coefficient of determination ($R^2$) and sample size

text: Is there any relationship between $R^2$ and sample size - does the $R^2$ increase with sample size? And does the adjusted $R^2$?


title: Is it possible to do data analysis in Open Office Calc?

text: I have recently come across the mircosoft data analysis tool add-on, and I was wondering if there was something comparable that can be used with OpenOffice? However, if there is nothing for OpenOffice, is there any other software that can be used in Ubuntu-Linux that can do similar things?


title: How to standardize a data-set

text: I have a data-set consisting of N p-dimensional observations (all quantitative variables). I want to apply a hierarchical clustering algorithm to those data. As explained on page 505 in Elements of Statistical Learning, when using weighted average to combine the distances of the individual variables, it is often desirable (I have found no clues to the contrary in my scenario), to set the weights for each variable such that all variables have equal influence on the distance of the observations (formulas can be found in the book). The problem is, that I use scipy.spatial.distance for calculating the distances, and this does no let me specify any weights. My question is, if I standardize my observations (multiply each dimension with 1/average of that dimension), does that solve the problem?


title: What distribution is suited to modelling the difference between bus schedule and arrival times?

text: I am looking for a one-dimensional distribution that can be fit to bus-delays. The data are real numbers (positive and negative), and exhibit both skewness and kurtosis. Which distribution should I use to fit these data? I would like to compute MAP or ML estimates, and I am using Matlab.


title: How to transfer a table into a chart?

text: I have the following table and I would like to make a representative chart.                                          Start      End        Duration     Implementation of project  01/05/2012 30/04/2014 24     XRD and Polymerisation     01/05/2012 01/06/2012  1     DFT                        01/05/2012 01/11/2012  6     Supervision                01/09/2012 01/03/2014 18     Fund raising               01/02/2013 01/03/2014 13     Proposal writing           01/02/2013 01/03/2014 13     Conferences                01/12/2013 30/04/2014  5     Future post                01/11/2013 30/04/2014  6      Can you help me to make a nice chart?


title: How 'good' are Holt-Winters forecasts with unusual alpha, beta and gamma values?

text: I'm using this python script for Holt-Winters forecasting (https://gist.github.com/andrequeiroz/5888967) that I believe chooses values of alpha, gamma and beta via RMSE optimisation. Sometimes the output gives strange values, 1.0 and 0.0 are very common. However, I know from some literature that values should typically fall between 0.02 and 0.2. How should I approach this? Should I force values between the 0.02 and 0.2 range?


title: Post-hoc test after 2-factor repeated measures ANOVA in R?

text: I have problems finding a solution regarding how to run a post-hoc test (Tukey HSD) after a 2-factor (both within-subjects) repeated-measures ANOVA in R. For the ANOVA, I've used the aov -function:               summary(aov(dv ~ x1 * x2 + Error(subject/(x1*x2)), data=df1))      After reading answers to other questions, I gathered that I would first have to re-run the ANOVA using some other fuction (e.g., lme). This is what I came up with.               Lme.mod <- lme(dv ~ x1*x2, random=list(subject=pdBlocked(list(~1, pdIdent(~x1-1), pdIdent(~x2-1)))), data=df1)     anova(Lme.mod)      Both main effects were significant, but there were no interaction effects. Then, I used these functions for the post-hoc comparisons:               summary(glht(Lme.mod, linfct=mcp(x1="Tukey")))     summary(glht(Lme.mod, linfct=mcp(x2="Tukey")))      However, there were some problems: First of all, the R Help file says that "The mcp function must be used with care when defining parameters of interest in two-way ANOVA or ANCOVA models (...) multcomp version 1.0-0 and higher generates comparisons for the main effects only, ignoring covariates and interactions (older versions automatically averaged over interaction terms). A warning is given." And sure enough, I received the following warning message:               Warning message:     In mcp2matrix(model, linfct = linfct) :     covariate interactions found -- default contrast might be inappropriate      Another puzzling thing was that although both main effects were significant, there were no significant differences in the post-hoc comparisons for one of the factors (x1). I've never encountered this before. Are the scripts/analyses correct/appropriate, or is there something that I'm missing? Any help would be most appreciated!


title: High SRMR despite good fit based on other indices in SEM (latent growth curve)

text: **Context** : Latent Growth Curve Modeling for continuous variable with 15 time points. One `intercept`, two `slopes` (for first and second half of all time points). N=146, >90% data coverage, `MLR` estimator in MPlus, no covariates. **Problem** : High SRMR=0.20 despite good fit based on other indices: chi- square goodness of fit p=.08, RMSEA=0.04, CFI/TLI=0.97. **My beginner's attempt** : There is one time point where standardized residual z-score for the mean is significant (freeing that time point didn't help resolving this problem), and another time point that had large and significant z-score for standardized correlation residuals with other time points, the combination of the two might have driven SRMR high?? **Question** : I can't think of any theoretical reason why this is happening... would this affect the overall interpretation of the model? What can I do to fix it?


title: What kind of experiment might Hiawatha have designed?

text: I am trying to understand the poem Hiawatha Designs an Experiment because it looks like the kind of in-joke that it would be nice to get, as a statistician. Here is what I would like to know: **Unhelpful answer:** "The poem is about the bias-variance tradeoff." **Helpful answer:** The poem is critical of the following kind of experimental design for the following reasons: .... **Great answer:** Here is a detailed set of Cliff-like notes explaining all the terms used in the poem (Anderson, Bancroft, truncated normals, etc.) Here are some links to books or papers espousing the kind of experimental designs that are being criticised. **Ideal answer:** Further to the info from "Great answer", here is an actual example of a study which might be criticised by the writer of the poem. I'd also be really interested in some historical context. Is the poem commenting on some sort of trend in statistics which was later reversed? Or is it still relevant?


title: Using multiple imputation for Cox proportional hazards, then validating with rms package?

text: I've been researching the mice package, and I haven't yet discovered a way to use the multiple imputations to make a Cox model, then validate that model with the rms package's `validate()` function. Here is some sample code of what I have so far, using the data set `veteran`:               library(rms)     library(survival)     library(mice)          remove(veteran)     data(veteran)     veteran$trt=factor(veteran$trt,levels=c(1,2))     veteran$prior=factor(veteran$prior,levels=c(0,10))          #Set random data to NA      veteran[sample(137,4),1]=NA     veteran[sample(137,4),2]=NA     veteran[sample(137,4),7]=NA          impvet=mice(veteran)     survmod=with(veteran,Surv(time,status))          #make a CPH for each imputation     for(i in seq(5)){         assign(paste("mod_",i,sep=""),cph(survmod~trt+celltype+karno+age+prior,             data=complete(impvet,i),x=T,y=T))     }          #Now there is a CPH model for mod_1, mod_2, mod_3, mod_4, and mod_5.      Now, if I were just working with one CPH model, I would do this:               validate(mod_1,B=20)      The problem I'm having is how to take the 5 CPH models (1 for each imputation), and be able to create a pooled model that I can then use with `rms`. I know that the `mice` package has some built-in pooling functions but I don't believe they work with the `cph` object in `rms`. The key here is being able to still use `rms` after pooling. I looked into using Harrell's `aregImpute()` function but I'm having some trouble following the examples and documentation; `mice` seems simpler to use.


title: NP-Hard optimisation problems that require approximate methods

text: What are some examples of NP-hard optimisation problem that requires approximate methods (such as Monte Carlo? I have done a lot of research but I can't find a suitable problem to implement apart from the Travelling Salesman Problem


title: Line graph has too many lines, is there a better solution?

text: I'm trying to graph the number of actions by users (in this case, "likes") over time. So I have "Number of actions" as my y-axis, my x-axis is time (weeks), and each line represents one user. My problem is that I want to look at this data for a set of about 100 users. A line graph quickly becomes a jumbled mess with 100 lines. Is there a better type of graph I can use to display this information? Or should I look at being able to toggle individual lines on/off? I'd like to see all the data at once, but being able to discern the number of actions with high precision isn't terribly important. # Why I'm doing this For a subset of my users (top users), I want to find out which ones may not have liked a new version of the application that was rolled out on a certain date. I'm looking for significant drops in the number of actions by individual users.


title: Which classifiers work well with unbalanced data?

text: I have a binary classification problem which is very unbalanced - it can have 98% of data from one class. Which classifiers work well with this sort of data? I have an unlimited supply of training data, since I produce it using a pseudo random number generator. However, I found that to get a neural network to produce decent results, I had to generate balanced (50:50) data. This is the equivalent of over-sampling. The problem with this approach is that the training data is then not representative of real life.


title: Given the location and scale parameters of a Gumbel distribution for variable X, how does one calculate the mean and variance of X^2?

text: I am working with predictive models for wind speeds, which have been given as Gumbel distributions. I need to convert the wind speeds to wind pressures using the formula: $Pressure = Density * Velocity^2$ So my question is how do I determine the distribution parameters for the wind pressure given the distribution of the velocity? Also, can I assume that $X^2$ also has a Gumbel distribution? The problem is fairly easy if I can just assume that it has a normal distribution and that: $E(X^2) = \mu^2 + \sigma^2$ $var(X^2) = E(X^4)-(E(X^2))^2$


title: Technical term for a survey scale with opposite values at both ends?

text: I am having trouble finding the technical term for a scale used in survey research with opposite values at both ends of the scale. For instance, how can I refer to a scale of the following kind: In your opinion, who is responsible for health care? Government (1) --- (2) --- (3) --- (4) --- (5) Citizens?


title: Total Probability and Events Negation

text: I would like to ask if I have an event when flipping a coin for $n$ times until a **first** head is observed. Event $E$: the first head on even numbered flip Event $H$: head on first flip Let $P(H) = p$ I am required to find $P(E)$ I know $P(E) = P(E|H)P(H) + P(E|H^c)P(H^c)$ where $H^c$ be the negation or complement of $H$; and $P(E|H)=0, P(H) = p, P(H^c)=1-p$. But I am stuck at $P(E|H^c)$. The solution stated that $P(E|H^c)=P(E^c)=1-P(E)$ but I do not understand. What is the negation or complement of $E$, let's denote is as $E^c$, in plain English? My assumption of $E^c$ is that it refers to first head on odd numbered flip. However I don't see why the probability of having a first head on even-flipped while having a tail on first flip equals to the probability of having the first head on odd numbered flip. Can anyone enlighten me?


title: required: good and straightforward method to detect change points in dependent univariate time series using r

text: I know that there are many related threads, packages and papers. Currently I`m reading through many of them. However, I don't plan to dig too deeply into this topic. I need a sound method that works with strongly dependent univariate time series (in large samples). My champion so far is the Kokoszka & Leipus test as evaluated and discussed in this paper ANDREOU (2002) DETECTING MULTIPLE BREAKS IN FINANCIAL MARKET VOLATILITY DYNAMICS but there seems to be no implementation of this test in R? Maybe someone knows better? I've spotted the Inclán and Tiao (1994) test (which is also discussed in this paper) in the `changepoint` package. But this method is too sensitive to outliers. It is originally developed for independent data, so additional transformations might be required first. The following paper also develops an interesting method and offers R code: Fryzlewicz (2014) Multiple-change-point detection for auto-regressive conditional heteroscedastic processes Currently I'm also checking the `bcp`package and the `strucchange` package (while the latter seems to be made for roughly independent data). I`m a bit lost between all these alternatives. And due to deadline constraints I have not the time to evaluate all these different possibilities (and many more that I have not listed here). Is there an R user who can recommend a certain method based on his own experience? The application will be to financial returns for a project at uni.


title: Fisher's exact test, contingency tables

text: Are there any other methods for an $m\times n$ contingency table with $m$ or $n$ greater than 2 for use with small samples ($np<5$) other than Fisher's exact test?


title: What is the minimum sample size for discrete choice experiment?

text: I was wondering if there is a minimum sample size for conducting discrete choice experiment. From what I know, if choosing the number of sample size is a problem, one can resort to using the magic number of 400+. Although it would be nice to have such sample size, but then this kind of experiment is expensive, so 400+ or more may be impractical. I have read several journal articles about DCE and I was surprised that their sample sizes did not even reached 400.


title: Interpreting standardized coefficients on a natural log response variable in OLS multiple regression

text: I am working on a housing problem in which I use dichotomous and ratio data to predict housing production (units constructed in a year-ratio) in a 17 year time period. At this time, I am using OLS and as I get better at stats, I shall attempt this problem using time-series analysis. That said, I have used R to standardize all of my ratio predicting data and left the dichotomous data raw. And I have also transformed the response variable to a Natural log to normalize the distribution (i.e. many, many zeros>>yes, I know Poisson or Zero-populated counts in the future). I have read the post on "interpret coefficients from a quantile regression on standardized data" and also the "convert my unstandardized independent variables to standardized." Based on those, I think that can do the following interpretation based on the following output. The variable `region_id` is dichotomous, `supply` is standardized.               Coefficients:                           Estimate Std. Error t value Pr(>|t|)         (Intercept)          2.687e+00  2.171e-01  12.379  < 2e-16 ***          region_id            1.805e+00  1.383e-01  13.049  < 2e-16 ***          supply              -2.205e+01  2.204e+00 -10.005  < 2e-16 ***      Region Interpretation:   For every on city that is located in the Houston region, you can expect that annual housing production will increase by 1.8%. Supply Interpretation:   For every one-unit increase in the standard deviation of housing supply, you can expect that annual housing production will decrease by -22.05%. Nota bene.   I am not a stats or math person at all, but I have been using R for the past three years and I am quite familiar with OLS, but if you throw up an equation it will look "appropriately" Greek to me. :)


title: simple sampling method for a Kernel Density Estimator

text: I have developed a simple _Kernel Density Estimator_ in Java, based on a few dozen points (maybe up to one hundred or so) and a Gaussian kernel function. The implementation gives me the PDF and CDF of my probability distribution at any point. I would now like to implement a simple _sampling method_ for this KDE. An obvious choice would of course be to draw from the very set of points making up the KDE, but I would like to be able to retrieve points that are slightly different from the ones in the KDE. I haven't found so far a sampling technique that I could easily implement to solve this problem (without depending on external libraries for numerical integration or complex computations). Any advices? I don't have specially strong requirements when it comes to precision or efficiency, my main concern is to have a sampling function that works and can be easily implemented. Thanks!


title: Is it valid to take a mean of p-values during cross-validation, when comparing the predicted output of a model to the actual output?

text: I am doing a cross-validation study, training a model on an input to predict a target. During training, my model generates an output vector that is guaranteed to be the same size as the corresponding training target vector. I can use metrics such as $R^2$ or RMS error to quantify this. During testing, my model can produce an output vector that is not the same size as the input (but they are the same order of magnitude). I'm wondering if there are any ways to quantify the similarity between the model output and the testing set targets. What I've come up with so far is to compare the _distributions_ under the null hypothesis that the model output distribution is the same as the test target distribution. I'm using things like the Kolmogorov–Smirnov test, Ansari- Bradley test, or a permutation test. For each cross-validation fold, there is 1 p-value. Is it valid to report a mean of p-values to summarize this? Or are there better ways to do this?


Recovering values from an estimated Pareto distribution A couple of weeks ago I asked a question about the Pareto distribution, wanting to understand how the units of measurement are related in order to make the usual statements about the invariance of the relationship between the percentage of people and the percentage of wealth. I got a really excellent, clear, thoughtful, and detailed answer to that question, here: Getting the units right for the Pareto distribution of wealth: x = people, dollars, dollars per person? Or rather, to the last half of my question (which _is_ what I emphasized). I now realize that the question was incoherent. It is silly to ask about units on a ratio of percentages, which is of course unitless. And I am still confused about the units. So I am reframing my question, as the answer to my previous question helped me to understand that the question I asked was not the one I really wanted to ask. (I always feel particularly indebted to those who help me to clarify my questions. (Especially relative to those who make fun of them, dismiss them as trivial, etc.)). Here is my second try. Suppose I have fit a Pareto distribution to a data set of income from individual households with income above a specified limit – call it the Pareto floor. It is the constant by which x is divided in the usual parameterization. Suppose, moreover, that incomes above a certain predetermined level (higher than the floor) have been censored. Call this level the Pareto ceiling. So this is not a true random sample, though it may be conditionally random given that there is censoring. (Actually they may censor a specified share of the population and select the ceiling level to fit. But let’s ignore that,) Please assume that I am fitting the parameters in a reasonable and intelligent manner, given the nature of the distribution and the censorship. Now, from the parameters of the Pareto that we have found, I want to determine both the number of people and the amount of income above and below any given income level. In addition to the Pareto parameters, you may assume that I have the following information: first, I know the Pareto base and ceiling; second, I know the total number of people above the Pareto base; and third, I know the total income belonging to all persons with incomes between the Pareto base and the Pareto ceiling (or at least I have a good estimate of it from a large sample). Now, taking it as given that the original assumption that the data are well- modeled by the Pareto is correct, assume we have people (conceptually) lined up by wealth or income from poorest to richest. I would like, for any income level, to be able to extract from the parameters of the Pareto and the other information given above, the number of people with incomes above and below that amount, and the number of dollars received by all the people above and below that amount. And I want to understand how to set up the problem so it comes out in people when I want people, and in dollars when I want dollars. I keep ending up with both of them proportional to the area under the CDF, and I know that’s not right. Warmest regard, andrewH


title: q-value and p-value

text: I have heard about q-values, in lines of p-values. What is the main benefit of using q-values over p-values? Why is it that p-values have to be adjusted? Wikipedia only has two lines of mention as follows: "The q-value is defined to be the FDR analogue of the p-value. The q-value of an individual hypothesis test is the minimum FDR at which the test may be called significant. One approach is to directly estimate q-values rather than fixing a level at which to control the FDR."


title: Quantifying homogeneity

text: I am hydrologist and I am looking for a solution for the following problem: Suppose that I have $k$ regions each yielding a sample of the same variable (e.g., annual peak flow) drawn from different sites within the given region (each sample $i=1,..K$ has a record length $n_i$). The homogeneity test (e.g., $K$ sample Anderson-Darling test) answers the question: is the region homogeneous? As a null hypothesis test, there are only two possibilities of answer: i.e., either the region is homogeneous ($F_1=F_2=...=F_k=F$) or it is heterogeneous. This binary answer is not flexible and do not allow me to compare homogeneity among several regions to point out the most homogeneous of them. So I am looking for a measure that overcomes this limitation. By analogy, this measure would be the same as the Akaike Information Criterion (AIC) for linear models. Indeed, the AIC compares different models in order to determine which one is the better to fit to data. Thus, this measure would compare different region and measure of how much sampling sites are alike within regions.


title: Which number is of value in linear regression: R Squared or P?

text: I am attempting to compare a two pitching stats (W/L% and ERA) to determine to what extent the latter can predict the former. After entering my data and performing the appropriate linear regression analysis I am left with two seemingly contradictory results. 1) I have an R-sq value of .257 which would lead me to state that the relationship is fairly weak 2) I have a P-value indistinguishable from 0 which would lead me to state that ERA is a near perfect predictor of W/L% What do I not understand about the meaning of these two numbers? What do they tell me about the correlation?


title: What is the expected partial value function really called?

text: If f is a pdf, the integral of x*f(x) over the entire range where f(x) > 0 gives, of course, the expected value. Suppose that integrate the same function, x*f(x) from negative infinity up to t, giving a new function, G(t). So G(t) = E(x | x < t), and is measured in the same units as x, as the expected value is (while probability measures like the CDF are unitless). I believe that G must be a well-known function with a name I should recognize and, on seeing it, clonk myself on the head and say “Of course!”. But I can not for the life of me think what it is. Am I right that this _is_a function with an established name? And if so, what is that name? I have some problems to solve in this area and without a name it’s hard to search the literature. Help me out here, please. Warmest regards, andrewH