title: Treating categorical variables in logistic regression in SAS

text: I am scoring my dataset using a logistic model. For getting the betas (coefficients) I used `proc logistic y = x1 x2 x3`. To tell the model which are class variables I am using`class` statement. Now, lets say if I have a categorical variable (with name ppsc), which has 4 categories, the betas are generated for top 3 categories (ppsc1, ppsc2, ppsc3 ) and I guess the fourth category is taken as reference. Now to score, I am using `proc score` to score, therefore I have to generate 3 binary variables ppsc1 ppsc2 and ppsc3 such that               if ppsc =1 then ppsc1= 1 else ppsc1 = 0.      ... ans same for ppsc2 and 3. Now my questions are:   1. If my category 4 is most important (number vise) and the logistic is generating betas for ppsc1, ppsc2 and ppsc3, what do I do? So for 4th category my logit(p) will be zero for that variable if the value of ppsc is 4? How do I handle this problem? Would this hold for other categorical variables too?   2. I dont want to make those ppsc1, ppsc2 and ppsc3 binary variables. Can't I specify at the time of scoring that my following variables are categorical variable by some class statement just like we do in `proc logistic`? I hope my problems are clear... I am BTW more concerned about 1st problem as that is something fundamental to understanding how the score is getting generated.


title: Two correlated DVs in a multiple regression

text: I have a multiple regression with an interaction term and I would like to run it on two separate DVs. The two DVs are quite highly intercorrelated (but conceptually different). How can I control for that intercorrelation?


title: Can Spearman's correlation be run on z scores?

text: I understand the logic of standardising raw data that is based on different scales into z scores so that they can compared, for example comparing a score of 75 out of 100 in one test versus 65 out of 120 in another test.   * Can I run the Spearman's correlation test on the z scores then? (i.e. forget about the raw data).    * How do I then relate the result back to the raw data?


title: Using information on both sides of a 'gap' in time series data for imputation

text: As with my previous question, I'm looking at ways to impute missing data in a hierarchical time series data. With al my other procedures, including the experimentation of imputation packages (`Amelia`, `HoltWinters` from `Forecast` and `MICE` imputation) I've only been able to use the time series data prior to the missing gap.                    Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec     2001 220 194 238 190 217 244 242 225 242 259 267 244     2002 212 246 250 236 261 286 265 269 226 267 234 246     2003 202 199 297 272 236 266 235 226 260 183 226 265     2004 211 215 219 213 240 236 273 266 262 244 241 235     2005 212 198 233 251 259 282 305 267 241 264 222 269     2006 182 220 250 287 279 281 286 332 300 272 221 233     2007  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA     2008 193 215 235 242 246 315 326 280 279 239 236 258     2009 246 189 257 241 268 223 260 288 234 260 216 195      I'm trying to do simple imputation procedure that uses forecasting and backcasting estimates from the time series model. Forecasting using prior data to predict the future and backcasting using the later data to “predict” the past. I would then like to combine the forecast and backcast value to use as imputation. After which I will look at the fit etc. How do I go about this in coding? For example, I'm able to determine what SARIMA model exist for the first period 2001-end2006. But not the full period (because my basic functions I know from R does not support the NA values.) This is only for the period 2001-end2006:               ARIMA(2,0,2)(1,0,1)[12] with non-zero mean           Call: auto.arima(x = ts.datt)           Coefficients:              ar1      ar2      ma1     ma2    sar1     sma1  intercept           1.3610  -0.8258  -1.2407  0.9191  0.8982  -0.7560   244.8374     s.e.  0.0884   0.0960   0.0878  0.1127  0.2190   0.3335     6.1894          sigma^2 estimated as 605.9:  log likelihood = -335.01     AIC = 686.02   AICc = 688.3   BIC = 704.23      Should I just model the first period, predict by `forecast`; model then the last period separately and then backcast? How will I do this backcasting (ie. 'predicting' the past)? **EDIT:** What I'm asking: 1) How do I use the data from years 2008 & 2009 to BACKCAST? I already know how to use 2001-2006 to forecast. 2) How do I determine the SARIMA model for the whole period? (2001-2009) ie.


title: Expected number of uniques in a non-uniformly distributed population

text: I am wondering if there is any reasonably simple way of calculating the following problem: Drawing, with replacement, $n$ balls from a bin of $N$ different colored balls, with a known probability of drawing each color of ball, what is the expected number of "unique" balls, _i.e.,_ balls with no other ball of the same color? _e.g._   $P(red) = 0.25$   $P(blue) = 0.3$   $P(green) = 0.2$   $P(yellow) = 0.25$ Some example outcomes with 5 balls: $\{red, red, green, blue, yellow\}$ - 3 unique balls   $\{red, red, green, green, blue\}$ - 1 unique ball   $\{blue, blue, blue, yellow, yellow\}$ - 0 unique balls Or, with 3 balls: $\{red, green, blue\}$ - 3 unique   $\{red, red, blue\}$ - 1 unique   $\{red, red, red\}$ - 0 unique For 1 ball, it's trivially 1; for 2 balls, it's 1 - the probability of the outcomes where the two balls are the same color * 2 balls, after that it starts getting more complicated.


title: Get the number of parameters of a linear model

text: Is there a way to get the number of parameters of a linear model like that?               model <- lm(Y~X1+X2)      I would like to get the number 3 somehow (intercept + X1 + X2). I looked for something like this in the structures that `lm`, `summary(model)` and `anova(model)` return, but I didn't figure it out. In case I don't get an answer, I'll stick on               dim(model.matrix(model))[2]      Thank you


title: How to find probability distribution of a multi-attribute datapoint in a dataset

text: I have a dataset with certain number of multi-attribute tuples. Each of the attribute values is a continuous random variable. I want to model each tuple (rather each attribute of the tuple) by a probability distribution function. The ultimate aim is to cluster the tuples. Can anyone suggest to me the best method to represent probability distribution of each tuple of the dataset.


title: Regression clustering

text: I am looking for references about classical methods in regression clustering. My problem is the following: I have a cloud of points that are assumed to have been generated by inverse functions with different coefficients. So that I have a family of regression functions. I only know the form of the model and i don't know the number of clusters. I would like to implement an algorithm for best fit using n regression functions of same form the inverse function $y = \frac{c}{x}$ with $c$ some constant. Do you know most popular methodologies? Do you know methods especially suited for this type of regression functions family?


title: Testing normality and independence of time series residuals

text: The simplest form of a white noise process is where its observations are uncorrelated. We can check this by applying e.g. a portmanteau test such as Lung - Box or Box - Pierce. The series might be Gaussian white noise where the observations are uncorrelated and also normally distributed and hence independent. We can test this with a normality test and a portmanteau test. As far as I know there is a third case where the observations are uncorrelated and independent without being normally distributed. In that case how can we test whether the observations are independent? Is there a statistical test for this?


title: statsmodels: error in kde on a list of repeated values

text: I am using Python's statsmodels module to plot a **violin/bean plot** of some data. I get the error _"LinAlgError: singular matrix"_ in the **KDE** calculation whenever a single violin plot is drawn from a list of **repeated instances of the same number**. Regular boxplots work just fine. Is this a bug or is there a reason behind this?               /home/me/myscript.py in plot_data_in_violinplot()        3435         plot_opts={'violin_fc':(0.8, 0.8, 0.8), 'cutoff':True,                               'bean_color':'#FF6F00', 'bean_mean_color':'#009D91',                              'bean_median_color':'b', 'bean_median_marker':'+'}     -> 3436         sm.graphics.beanplot(data_to_boxplot, ax=ax, jitter=jitter, plot_opts=plot_opts)               /usr/local/lib/python2.7/dist-packages/statsmodels/graphics/boxplots.py in beanplot(data, ax, labels, positions, side, jitter, plot_opts)         333     for pos_data, pos in zip(data, positions):         334         # Draw violins.     --> 335         xvals, violin = _single_violin(ax, pos, pos_data, width, side, plot_opts)         336          337         if jitter:          /usr/local/lib/python2.7/dist-packages/statsmodels/graphics/boxplots.py in _single_violin(ax, pos, pos_data, width, side, plot_opts)         170     pos_data = np.asarray(pos_data)         171     # Kernel density estimate for data at this position.     --> 172     kde = gaussian_kde(pos_data)         173          174     # Create violin for pos, scaled to the available space.          /usr/local/lib/python2.7/dist-packages/scipy/stats/kde.pyc in __init__(self, dataset, bw_method)         186          187         self.d, self.n = self.dataset.shape     --> 188         self.set_bandwidth(bw_method=bw_method)         189          190     def evaluate(self, points):          /usr/local/lib/python2.7/dist-packages/scipy/stats/kde.pyc in set_bandwidth(self, bw_method)         496             raise ValueError(msg)         497      --> 498         self._compute_covariance()         499          500     def _compute_covariance(self):          /usr/local/lib/python2.7/dist-packages/scipy/stats/kde.pyc in _compute_covariance(self)         507             self._data_covariance = atleast_2d(np.cov(self.dataset, rowvar=1,         508                                                bias=False))     --> 509             self._data_inv_cov = linalg.inv(self._data_covariance)         510          511         self.covariance = self._data_covariance * self.factor**2          /usr/local/lib/python2.7/dist-packages/scipy/linalg/basic.pyc in inv(a, overwrite_a, check_finite)         381         inv_a, info = getri(lu, piv, lwork=lwork, overwrite_lu=1)         382     if info > 0:     --> 383         raise LinAlgError("singular matrix")         384     if info < 0:         385         raise ValueError('illegal value in %d-th argument of internal '     LinAlgError: singular matrix


title: Best way to combine binary and continuous response

text: I am trying to come up with the best way to predict payment amount for a collections agency. The dependent variable is only non-zero when a payment has been made. Understandably, there are an overwhelming number of zeros because most people cannot be reached or cannot pay back the debt. There is also a very strong negative correlation between the debt amount and the probability of making a payment. Typically, I would create a logistic model to predict the probability of pay / non-pay but this has the unfortunate consequence of finding the people with the lowest balances. Is there a way to combine a logistic pay/non-pay model with a separate model predicting the payment amount?


title: I need to convert a yearly data into a quarterly and monthly data?

text: Can anybody please help me convert yearly data into monthly and quarterly data?


title: Use of a negative binomial model for fitting alternative splicing event

text: I am working on RNA-Seq data (on alternative splicing). Let's say I am looking at a particular type of alternative splicing event - _exon skipping_. For each intron (or junction), I look if it is _normally spliced_ or if there is an _exon skipping_ event happening. I have three biological replicates ( for each junction, I have 3 set of values). To sum it up, my data set looks like this:               Junction 1:               Rep1    Rep2    Rep3     exonSkip    8      0       0      normal      12     6       8              Junction 2:               Rep1    Rep2    Rep3     exonSkip    5      9       8      normal      58     60      44         ....     ....       My objective is, for each junction, from these replicates, to find out if that particular exon-skipping event is _statistically significant_. Initially, I summed up all the values for _exonSkip_ and _normal_ separately (in the first case, 8 and 26) and then concluded there are at least 2 exonSkip events. However, I came to know its not the best and that there are better ways. 1) From the literature (on gene expression), I came to know that the biological replicates should be used to obtain an _estimate_. Since these are reads = count data, and they happen to have a high degree of dispersion among replicates, a _negative binomial distribution_ is suggested. So, I used a `glm.nb` model from the `MASS` library as follows:               # R-code      # Junction 1        require(MASS)       dat1 <- data.frame(y=as.numeric(c(8,0,0,12,6,8)),                         exonSkip=as.factor(c("yes","yes","yes","no","no","no")))       out1 <- glm.nb( y ~ exonSkip, data=dat1)       summary(out1)     # Junction 2        dat2 <- data.frame(y=as.numeric(c(5,9,8,58,60,44)),                         exonSkip=as.factor(c("yes","yes","yes","no","no","no")))       out2 <- glm.nb( y ~ exonSkip, data=dat2)       summary(out2)      For Junction 1: I got `p=0.143` for exonSkipyes. **Question:** Does this mean that the nb model fit can not be trusted? For Junction 2: I got `p<2e-16`. However, the test `glm.nb` gave a               *warning: In theta.ml ... Iteration limit reached*.       **Question:** Is this okay? Now, from the model, I then used the `predict` function to estimate the values for _exonSkip_ and _normal_ from these 3 replicates with the model fit.               exp(predict(out1, dat1))       # Result: 2.6667 and 8.6667 = 3 and 9       exp(predict(out2, dat2))       # Result: 7.333 and 54.000 = 7 and 54        **Questions:**   1. Is this method of estimating, assuming negative binomial glm model right? Particularly, the use of `predict` function.    2. Can I still use `predict` (as I have used above) in case of a non-significant fit? If not, then how else can I get an estimate?  2) Even if manage to get the estimate, I have again 2 numbers: 1 for _exonSkip_ and other for _Normal_. From here, I would like to obtain a measure (or p-value) of how significant it is. How I can go about this? I think the p-value from the `glm.nb` is how significant the model fits the data...   One way I thought of is: If there were totally _X+Y_ ( _X_ = total _exonSkip_ and _Y_ = total _Normal_ ) events, then from a sample of _b_ events, if I get _a_ exonSkip, then it would follow a _hypergeometric distribution_ and I could obtain the p-value as,               sum(dhyper(a:b, X, Y, b))      Is this correct? I'd appreciate any feedback.


title: Multiple determination in terms of simple determination

text: I'm working through _Applied Linear Regression Models_ (Kutner) and having a hard time figuring out this problem (it's actually problem 6.26): > Given a First-Order Model with Two Predictors: $Y = \beta_0 + \beta_1X_1 + > \beta_2X_2 + \epsilon$ , show that the coefficient of simple determination > between $Y$ and $\hat{Y}$ equals the coefficient of multiple determination > of $R^2$. I've also seen this problem in other texts as well - it seems to be a common problem, yet I have not been able to find any kind of proof online. Does anyone have a link to or can someone show me how to show they are equal?


title: Selecting regression type for Dickey-Fuller test

text: I have a time series and I want to check whether it has a unit root or not. I want to use the Dickey-Fuller test. How should I choose the regression equation from the three options: regression without constant and trend, with constant (drift) and with constant and trend? Is there a procedure that I should follow to select the regression? On what criteria is the choce of regression is based? If the errors of the chosen regression are correlated then should I run the augmentd DF test with the same regression chosen in the first step?


title: Which AIC value to use from R's sarima() function for model comparison

text: I'm using R's 'astsa' package and I get the following output from sarima. Which AIC value would I use to compare this model (let's call it A) against others? When trying another model (B), model A's fit$AIC (858.19) is greater than model B's, but model A's AIC (12.38841) is less than model B's, so I'm not sure which model to choose. What's the difference between the two AIC's, AICc's, and BIC's? I've checked the sites below, among others, but haven't been able to figure it out. Any help is much appreciated. http://stat.ethz.ch/R-manual/R-patched/library/stats/html/arima.html http://www.inside-r.org/packages/cran/astsa/docs/sarima               $fit     Series: xdata      ARIMA(0,1,1)(1,1,1)[12]                              Coefficients:              ma1    sar1     sma1          -0.3282  0.5529  -0.8835     s.e.   0.3290  0.4751   0.8635          sigma^2 estimated as 82513:  log likelihood=-425.1     AIC=858.19   AICc=858.93   BIC=866.5          $AIC     [1] 12.38841          $AICc     [1] 12.42448          $BIC     [1] 11.48327


title: Bayesian vs Maximum entropy

text: Suppose that the quantity which we want to infer is a probability distribution. All we know is that the distribution comes from a set $E$ determined, say, by some of its moments and we have a prior $Q$. The maximum entropy principle(MEP) says that the $P^{\star}\in E$ which has least relative entropy from $Q$ (i.e., $P^{\star}=\displaystyle \text{argmin}_{P\in E}D(P\|Q)$) is the best one to select. Whereas the Bayesian rule of selection has a process of selecting the posterior given the prior which is supported by Bayes' theorem. My question is whether there is any connection between these two inference methods (i. e., whether the two methods apply to the same problem and have something in common)? Or whether in Bayesian inference the setting is completely different from the above mentioned setting? Or am I not making sense?!


title: Standardised or unstandardised?

text: I've been trying to teach myself some statistics (all social science-related) using SPSS. It's all going well so far, I'm amazed at how much I'm managing to understand given how terrible I am at maths. However, there is one thing I don't fully understand. Say I run a multiple OLS regression. Then SPSS will give me a thing called a 'standardized coefficient' and another called an 'unstandardized coefficient'. I understand that, essentially, the standardised coefficient is useful if you want to know which of the different independent variables is more important. The unstandardised coefficient, on the other hand, is the one you use to find out what effect 1 unit change in x will have on y. The practical question is this: if I'm writing a paper using a multiple regression, I will have to complete a table with all the results. The problem is that, in this paper, I would want to comment on the standardised coefficient (to compare effects of variables) **as well as** the unstandardised one (to calculate unit changes). So essentially, which coefficient am I supposed to report in a paper in which I make reference to both? Should I include one in the results table and the other one somewhere else? Ta!


title: Calculating joint probabilities from conditional probabilities

text: I need to calculate the probability of $P(L,D)$, where $L$ and $D$ are not independent. I have estimated $P(L)$ and $P(D)$ with two distinct models and I also know $P(L|D)$ and $P(D|L)$. As far as I can see it I have two different estimates of $P(L,D)$: $P(L,D) =P(L) * P(D|L)$ $P(L,D) =P(D) * P(L|D)$ Am I best off just averaging these two estimates for the joint $P(L,D)$?


title: Probability of negative return & target return

text: I need to implement the below two things: Probability of negative return formula link http://www.styleadvisor.com/content/probability-negative-return & Probability of Target Return formula link http://www.styleadvisor.com/content/probability- target-return. If I have just 12 or 24 monthly portfolio return then inside CCM or CCV calculation they are taking E[R]; which ideally means expected return but then how to calculate that. One possible way is if I receive 12 or 24 same period benchmark return and then calculate Portfolio Beta using variance-covariance method and put it in CAPM equation to calculate E[R]. Later on why/how Standard normal cumulative distribution have to be applied on them. As how is concerned, I think I can apply NORMSDIST function of excel to show the demo result, later on use R/SAS for implementation. Can somebody tell my thought process/approach is correct. My question is more on why because standard normal cumulative distribution mu and sigma is fixed, i.e. 0 and 1 so how the expectation/probability factor getting captured. Any clarification would be great.


title: What is the statistical efficiency of L-moments?

text: In particular I am interested in the scale estimator. Hopefully it is much better than that of IQR.


title: Trouble understanding two scale graph

text: I am trying to plot a graph which has two y-axes and one x-axis. Consider three variables A,B,C. My first data set G1, has values of A and B. While getting this dataset in the script I kept C constant. In second dataset G2, I have B and C and was generated for some fixed value of A. Now I am plotting this with `matplotlib.pyplot`. My question is if I keep Y1=A, X=B, Y2=C and plot this data. It comes out to be weird looking which I think is because for G1, C has only one value and for G2 A has only one value. I get two different lines in the final graph and I do not understand what this graph is showing? To me it feels that if I plot graph between A,B,C variables then it should show me the relationship between them rather than plotting two independent different lines. ![http://imgur.com/OahwFir](http://i.stack.imgur.com/JqK6Q.png)


title: Standard error of the mean vs. standard error of mean difference

text: When plotting within-subject data for condition A vs. condition B (significance tested through a t-test), should the error bars reflect standard error of the mean for each condition independently? Or should they be SE of the mean difference between A and B? When there are more conditions for an ANOVA -- A vs. B for condition 1 and A vs. B for condition 2 -- should the error bars correspond to the SE of mean difference of A-B at 1 and then the SE of mean difference at 2? **Edit:** This paper gives a way to calculate SE/CI for within-subject designs. The gist is that the subjects are normalized to reduce the between-subject contribution to the error bars and better reflect the results of a repeated measures ANOVA. Cousineau, D. (2005). Confidence intervals in within-subject designs: A simpler solution to Loftus and Masson’s method. Tutorial in Quantitative Methods for Psychology, 1(1), 42–45. PDF


title: Generating random matrices with sum and maximality constraints

text: I'd like to generate a random square matrix such that the rows are normalized to one and the diagonal elements are the maximum of their column. If there an efficient way to sample these matrices uniformly? $2 \times 2$ matrices are straightforward to create. The first column is generated by sampling two independent uniforms from $[0,1]$ and moving the maximum to the first entry. The second column is then the complement of the first. My attempts to follow an analogous procedure in higher dimensions feel very ad-hoc, and most troublesomely, give a very different distribution for the final complemented column than the others. Is some type of rejection sampling going to be my best bet? The motivation is to generate likelihood matrices, with the rows being states of the world and columns being observed signals. Each signal is most likely to occur in its corresponding state, but isn't necessarily the most likely signal in that state. **Edit:** The exact distribution doesn't matter, other than it seem natural. I would like to avoid the diagonal also being maximal for the rows though. Here is my own previous ah-hoc attempt.               gen.likelihood.matrix <- function(T){         # T: Number of states              mat <- matrix(runif(T^2), nrow=T) # Intialize with uniform variates         diag(mat) <- 1 # To guarantee diag is max in column              candidate.mat <- matrix(nrow=T, ncol=T)         # Loop to check for constraint satisfaction         for(k in 1:100){             weights <- runif(T-1)             weights <- weights / sum(weights)                  for(i in 1:T){                 # Rescale columns                 candidate.mat[i,-T] <- mat[i,-T] * weights                 # Replace last column with complement to meet row constraint                 candidate.mat[i,T] <- 1 - sum(candidate.mat[i,-T])             }                  # Does the last column meet the maximality constraint?             if(candidate.mat[T,T] > max(candidate.mat[-T,T])){                 # Shuffle to minimize the distortion from last column                 reorder <- sample(1:T,T)                 return(candidate.mat[reorder,reorder])             }         }              # Recurse to get new initial values if necessary         return(genLikeliMatrix(T))     }      This works in terms of giving me matrices efficiently, but I get a weird bimodal distribution on the diagonal due to how I handle the last column. ![Distribution of diagonal elements](http://i.stack.imgur.com/5QPuS.png)


title: Plot a subset of categories on the x-axis in ggplot

text: everyone! I'm trying to plot pvalues of a test by categories (in this case, genetic loci). So, the x-axis contains gene names. I'm getting the plot to work, but the x-axis was very busy, so I am thinking of splitting the screen into several plots, each containing only a few loci (i.e., a few x values). When I attempted to plot a subset of my data with the script below, I got a plot that still contained far too many x-values on the x-axis. In other words, the bars that were plotted were all crowded onto the left-hand side of the plot, and many x-values had nothing. Any idea what I'm doing wrong? It's hard to read the loci names on the x-axis, so I can't tell whether it's repeating the same locus names (i.e., maybe a problem having to do with melt()) or whether it's plotting all of the Loci values from data instead of from split1_data. Alternatively, any different recommendations for how to make a plot with hundreds of bars easy to see? Thank you so much!               library(reshape)     library(ggplot2)     require(ggplot2)          setwd("/Users/markfisher/Desktop")          sink(file="/Users/markfisher/Desktop/Pvalue_HWE_output.txt")     data=read.csv("Pvalues_of_all.csv", header=TRUE)     attach(data)     print(data$Loci[1:10])           split1_data<-subset(data,data$Loci %in% data$Loci[1:10])          split1_datam<-melt(split1_data,id="Loci")     print("split1_data$Loci")         print(split1_data$Loci)      sink()          pdf('/Users/markfisher/Desktop/pvalue_sensitivity.pdf', bg = "white")     p <- ggplot(split1_datam, aes(x =Loci, y = value, color = variable, width=.15))     p + geom_bar(position="dodge") + ylab("P-value")+ geom_hline(yintercept=0.05)     dev.off()      Update: I added made the suggested change of using droplevels(), and got an error about an unexpected numeric constant. This is the output from my R console:               > library(reshape)     > library(ggplot2)     > require(ggplot2)     > setwd("/Users/markfisher/Desktop")     > sink(file="/Users/markfisher/Desktop/Pvalue_HWE_output.txt")     > data=read.csv("Pvalues_of_all.csv", header=TRUE)     > attach(data)     > print(data$Loci[1:10])          > split1_data<-droplevels(subset(data,data$Loci %in% data$Loci[1:10]))         > split2_data<-subset(data,data$Loci %in% data$Loci[11:20])         > split3_data<-subset(data,data$Loci %in% data$Loci[21:30])         > split4_data<-subset(data,data$Loci %in% data$Loci[31:40])         > split5_data<-subset(data,data$Loci %in% data$Loci[41:51])         > split6_data<-subset(data,data$Loci %in% data$Loci[52:62])         >          > split1_datam<-melt(split1_data,id="Loci")         > print("split1_data$Loci")     > print(split1_data$Loci)        >      >      > sink()     > pdf('/Users/markfisher/Desktop/pvalue_sensitivity.pdf', bg = "white")     >      > p <- ggplot(split1_datam, aes(x =Loci, y = value, color = variable, width=.15))1     Error: unexpected numeric constant in "p <- ggplot(split1_datam, aes(x =Loci, y = value, color = variable, width=.15))1"     > p + geom_bar(position="dodge") + ylab("P-value")+ geom_hline(yintercept=0.05)     Error: object 'p' not found     > dev.off()     null device      1      >      If anyone needs to see my output file, it looks like this:               The following object(s) are masked from 'data (position 3)':         All, Loci, X1_only, X1_removed, X2_only, X2_removed, X3_only, X3_removed,         X4_only, X4_removed, X5_only, X5_removed, X7_only, X7_removed, X78_only,         X8_only, X8_removed     The following object(s) are masked from 'data (position 4)':         All, Loci, X1_only, X1_removed, X2_only, X2_removed, X3_only, X3_removed,         X4_only, X4_removed, X5_only, X5_removed, X7_only, X7_removed, X78_only,         X8_only, X8_removed     The following object(s) are masked from 'data (position 5)':         All, Loci, X1_only, X1_removed, X2_only, X2_removed, X3_only, X3_removed,         X4_only, X4_removed, X5_only, X5_removed, X7_only, X7_removed, X78_only,         X8_only, X8_removed     The following object(s) are masked from 'data (position 6)':         All, Loci, X1_only, X1_removed, X2_only, X2_removed, X3_only, X3_removed,         X4_only, X4_removed, X5_only, X5_removed, X7_only, X7_removed, X78_only,         X8_only, X8_removed     The following object(s) are masked from 'data (position 7)':         All, Loci, X1_only, X1_removed, X2_only, X2_removed, X3_only, X3_removed,         X4_only, X4_removed, X5_only, X5_removed, X7_only, X7_removed, X78_only,         X8_only, X8_removed     The following object(s) are masked from 'data (position 8)':         All, Loci, X1_only, X1_removed, X2_only, X2_removed, X3_only, X3_removed,         X4_only, X4_removed, X5_only, X5_removed, X7_only, X7_removed, X78_only,         X8_only, X8_removed     The following object(s) are masked from 'data (position 9)':         All, Loci, X1_only, X1_removed, X2_only, X2_removed, X3_only, X3_removed,         X4_only, X4_removed, X5_only, X5_removed, X7_only, X7_removed, X78_only,         X8_only, X8_removed     The following object(s) are masked from 'data (position 10)':         All, Loci, X1_only, X1_removed, X2_only, X2_removed, X3_only, X3_removed,         X4_only, X4_removed, X5_only, X5_removed, X7_only, X7_removed, X78_only,         X8_only, X8_removed     The following object(s) are masked from 'data (position 11)':         All, Loci, X1_only, X1_removed, X2_only, X2_removed, X3_only, X3_removed,         X4_only, X4_removed, X5_only, X5_removed, X7_only, X7_removed, X78_only,         X8_only, X8_removed     The following object(s) are masked from 'data (position 12)':         All, Loci, X1_only, X1_removed, X2_only, X2_removed, X3_only, X3_removed,         X4_only, X4_removed, X5_only, X5_removed, X7_only, X7_removed, X78_only,         X8_only, X8_removed     The following object(s) are masked from 'data (position 13)':         All, Loci, X1_only, X1_removed, X2_only, X2_removed, X3_only, X3_removed,         X4_only, X4_removed, X5_only, X5_removed, X7_only, X7_removed, X78_only,         X8_only, X8_removed     The following object(s) are masked from 'data (position 14)':         All, Loci, X1_only, X1_removed, X2_only, X2_removed, X3_only, X3_removed,         X4_only, X4_removed, X5_only, X5_removed, X7_only, X7_removed, X78_only,         X8_only, X8_removed      [1] Baez       Blue       C147       C204       C21        C278_PT    C294       C316            [9] C485       C487_PigTa     62 Levels: Baez Blue C147 C204 C21 C278_PT C294 C316 C485 C487_PigTa C536 Carey ... Yellow     [1] "split1_data$Loci"      [1] Baez       Blue       C147       C204       C21        C278_PT    C294       C316            [9] C485       C487_PigTa     Levels: Baez Blue C147 C204 C21 C278_PT C294 C316 C485 C487_PigTa      I would also show what my plot looks like, but alas I do not yet have enough cool points on CrossValidated to do so. Thanks again for all of your help with this. I hope my updates clarify things a little bit...


title: Multi-class AUC in Matlab

text: I would like to compute the area under the ROC-courve (AUC) metric for a classifier with multiple classes. Do you know (reliable) functions for Matlab that implement methods for that, like e.g. in [1]? Thank you! [1] Fawcett, T. (2006). An introduction to ROC analysis. Pattern recognition letters, 27(8), 861-874.


title: Can I use a paired-samples t-test when the 2 groups that I compare have different numbers of data points?

text: I want to compare the effect of a drug (immunomodulator) on antibody levels in a group of patients $(n=18)$. Blood samples were taken pre- and post-treatment in each patient and the means of antibody levels compared for each study subject. At this point I guess that the best and simplest option is a "paired t test". However, in some cases, I have a different number of blood samples in the pre- vs. the post-treatment group. For example, for one patient I have 16 samples pre-treatment and 5 post-treatment.   * Can I still use the paired t test?    * Should I use the independent t test?


title: How to combine values based on standard errors?

text: Let's assume we have two means with corresponding standard errors.               mean_1 = 3.75, se_1 = 0.64     mean_2 = 2.90, se_2 = 0.94      I would like to average the means and account for the standard error. I would like `mean_1` to contribute more (since the standard error is lower) to the average than `mean_2`. A weighted mean if you will. One way would be to specify weights as the relative magnitude of the standard error like so (I use R's function `weighted.mean` found in the base `stats` package):               > weighted.mean(c(3.75, 2.0), w = c((1 -0.64/(0.64+0.94)), (1- 0.94/(0.64+0.94))))     [1] 3.041139      Does that make sense? How would one go about this problem? I would be happy to hear your critiques and suggestions.


title: How to calculate the probability of winning an election?

text: I'm analyzing a survey data that contains the following 3 variables (questions):   * If the election was today who would you vote for in the following list? (A, B, C, D, or E)? Let's imagine the respondent choose A in the previous item. Then, she is asked:   * What if A is not participating in the contest, who would you vote for? (B, C, D, or E)? Let's imagine the respondent choose C in the previous item. Then, she is asked:   * What if A and C are not participating in the contest, who would you vote for? (B, D, or E)? I have the answer for the three questions and I want to create an algorithm to help me calculate the chances of victory of A against all other candidates, the chances of B against all the others, and so on. My goal is to report who are the most dangerous opponents of each candidate. Anyone has an idea/suggestion?


title: Conditional Inference trees on correlated observations

text: I have a dataset containing a binary response variable and a few numeric predictor variables. I would like to use Conditional Inference Trees algorithm in R using the ctree function in party package. However, the outcomes in the dataset are correlated because a single individual is contributing to more than one record in this data. Is there an underlying assumption of independence in ctree that may prevent the use of this technique on my dataset? If so, are there any standard techniques by which I can get around this? Please note that I am using GEE (PROC GENMOD in SAS with the REPEATED option) as one of the modeling options. I would like to also provide an alternative solution using a tree-based approach and I would definitely prefer to use Conditional Inference Trees over CART of CHAID, if possible. Thanks.   * KC


title: What's the best way to visualize the effects of categories & their prevalence in logistic regression?

text: I need to present information about the main predictors of a candidate's votes using a public opinion survey data. I have run a logistic regression using all the variables that I care about, but I can't find a good way to present this information. My client doesn't care about the size of the effect only, but about the interaction between the size of the effect and the size of the population with such attribute. How can I deal with that in a graph? Any suggestions? Here is an example: > The $\beta$ of variable SEX (Male=1) when the dependent variable is Vote/Not > in a candidate is 2.3, which is a big number after having been exponentiated > and treated as odds ratio or probability. However, the society in which this > survey was run only had 30% men. Therefore, although man supported this > candidate quite a lot, their numbers are insignificant for a candidate > trying to win a majoritarian election.


title: How to measure the performance of a regressor?

text: I have a set of data with 6 feature and an attached score for each data point within the set. I have build a linear regressor for my data. The thing is that the output values are pretty huge, approximately 10000 and I obtain values between 9600 and 10000. How can I measure the error of my system such that it falls between 0 and 1?


title: Comparing Pareto fitting methods

text: I am comparing different method of fitting Pareto distributed random numbers. What seems very strange to me is that fitting a straight line in log-log scale has to be the worst numerical method and indeed it is, but when I superimpose the log-log histogram with that straight line they superimpose very good. The other fitting method superimpose very bad. How is it possible this paradox? I am doing some mistakes? Here is my code in R for check:                pareto.MLE <- function(X)      {         n <- length(X)         m <- min(X)         a <- n/sum(log(X)-log(m))         curve(-(a+1)*x+log(a)+a*log(m),add=T,col='black')         return( c(a) )       }           pareto.Alpha <- function (X)      {           m <- min(X)         mu <- mean(X)         a <- mu/(mu-m)         curve(-(a+1)*x+log(a)+a*log(m),add=T,col='green')         return( c(a) )       }                pareto.median <- function (X)      {           m <- min(X)         med <- median(X)         a <- log(2.)/log(med/m)         curve(-(a+1)*x+log(a)+a*log(m),add=T,col='blue')         return( c(a) )       }            pareto.fit <- function (X)            {         dd <-hist(X,plot=F)              for(i in length(dd$counts):1)           {      if(dd$density[[i]]==0)       {       dd$density[[i]]=dd$density[[i+1]]         }     }             cc <-lm(log(dd$density)~log(dd$mids))        curve(x*(coefficients(cc)[2])+(coefficients(cc)[1]),add=T,col='red')        return( c(-1*(coefficients(cc)[2]+1)))               }             library(PtProcess)        alpha=1.        hh1<-(matrix(rpareto(100,alpha,0.1),ncol=1))             pp<-hist(hh1,plot=F)        plot(log(pp$mids),(log(pp$density)))             parameter.1<-pareto.MLE(hh1)        parameter.a <- pareto.Alpha(hh1)        parameter.m <- pareto.median(hh1)        parameter.f <- pareto.fit(hh1)                  cat('Alpha: ',alpha,'\n')        cat('MLE: ',parameter.1,'\n')        cat('Mean: ',parameter.a,'\n')        cat('Median: ',parameter.m,'\n')        cat('Fit: ',parameter.f,'\n')


title: How can I get stepwise selection for multinomial regression in SPSS to process?

text: I am using SPSS to analyze my data. It is multinomial regression I am using. I have 14 predictor variables of which 9 are categorical and 5 continuous. The response variable is categorical with three categories. Then after pressing the appropriate buttons to reach at the multinomial regreesion the result displays but the error message I am facing is stepwise regression is not processed. Where do you think is the problem? And what is numerical problem and how can I solve it?. Can I reduce the number of variables step by step considering the VIF till it approaches to 1?


title: Simple random walk and sudden drift, how to detect the change?

text: I have got the following question relating to random walks. I would like to determine the moment when a random walk changes from being a simple random walk to start to drift at certain time. This sort of scenario could happen when temperatures start to rise or sea levels droppping. The idea is to detect that change as soon as possible. Clearly if the system is left to drift for a long time, it is fairly obvious it has changed from a given initial level (given by the simple random walk at the start), but I would like to know if it is possible to detect that change very quickly after it happens. What sort of analysis would be needed? Any suggestions greatly appreciated. Many thanks!! Michael and whuber, I really appreciate your comments. My initial question was in fact related to Quality Control, as I am trying to model a CUSUM as a Random Walk: When the system is in-control, I see it as a simple random walk, when it starts to drift and eventually goes out-of-control, I see it as a drifting random walk. Thus my question as to how (and critically, how quickly) to detect this change from stationary to drifting within random walks. I can see it is not easy to detect the drift, specially when the drift is small due to the stochastic nature of the system. But I think the ARIMA idea will help me, I had not looked at this theory before and thank you for this advice. Mili.


title: Time Series Decomposition : Box Cox for Additive Decomp

text: Coming from basically no time series back ground, this is likely a simple question, but what is the relationship between "being able to" use an additive decomposition of a series into seasonal, trend and remainder and a Box Cox transformation? From Professor Hyndman's blog: > Because not all data could be decomposed additively, we first needed to > apply an automated Box-Cox transformation. I was wondering: 1) What makes an additive decomposition attractive relative to a multiplicative one (which I understand is basically the other choice). 2) What is the requirement for an additive decomp and what does Box Cox do to make this possible? I think of Box Cox for ANOVA and reducing heteroskedasticity. Is there a tie in with decomposition of a series?


title: Missing values replacement?

text: I have a dataset of 500 people and am trying to fit a prediction model using both quantitative and categorical variables. I have adjusted the dataset as much as possible, but still have one variable (which is important in the analysis) with 19 missing cases. What would you say is the best approach/program to use to replace the missing values/instead of deletion? Would the EM estimator from SPSS provide good enough replacement? Seeing that I have important categorical predictors that will be used in the analysis, I understand that I can't use them in Amos. What other programs could be suited for this type of model? Thank you!


title: Estimating mass proportions from fractional samples

text: Say you are on a beach that is composed of a mixture of different colored sands (for simplicity, black and white) and you wish to estimate how much there is of each type. You have a bucket that you can use to take samples of the sand and a technique to separate the different types and weigh them. Say also that the beach is not too big, so you know the fractional amount of sand in your bucket samples compared to the total sand on the beach. Let the true fractional portion of white sand on the beach be called $w^*$, and therefore (in this two-type case) the true amount of black sand be $(1 - w^*)$. Our bucket samples give weights $\{w_i\} = w_1, w_2, …, w_n$. I'd like to estimate the distribution of $w^*$ given my measurements $\{w_i\}$. Assume the sand is iid uniformly at random across the beach and that the bucket samples are taken without replacement. Can anyone give me some direction on how to model this problem? I would ultimately like to generalize to more than two types of sand, but the two-type case, even with just a single sample, should at least get me started. I began trying to consider the sampling as a Bernoulli process, but I ran into trouble with understanding how to deal with the non-integer nature of the samples: i.e., the measurement values are continuous-valued masses, rather than some sort of count.


title: Slope Derivation for the variance of a least square problem via Matrix notation

text: I have a question to solve the following matrix problem: $$E[( (X'X)^{-1}X \epsilon )^T ((X'X)^{-1}X \epsilon )]$$ into the solution $$= \Sigma^{2} (X'X)^{-1}.$$ Where $\Sigma$ is the covariance matrix of $\epsilon$. My question is how to solve this problem using the properties of matrix algebra. This equation occurs when using least squares to solve a linear regression problem. (as can be seen in: http://en.wikipedia.org/wiki/Proofs_involving_ordinary_least_squares , section Unbiasedness and Variance of $\hat\beta$). I saw a solution for a specific case in this link: Expected Value and Variance of Estimation of Slope Parameter β1 in Simple Linear Regression. However, I was interested in the generalized solution described on the equation. That is what I obtained so far: Using properties of transpose, we obtain: $$E[( (X'X)^{-1}X \epsilon )^T ((X'X)^{-1}X \epsilon )] = E[( (\epsilon' X' ((X'X)^{-1})' )((X'X)^{-1}X \epsilon )].$$ $(X'X)^{-1}$ is a symmetric matrix, thus: $$( (X'X)^{-1} )' = (X'X)^{-1}.$$ And I obtain $$E[( (X'X)^{-1}X \epsilon )^T ((X'X)^{-1}X \epsilon )] = E[( (\epsilon' X' (X'X)^{-1})^{2} X \epsilon )].$$ I could not solve any further than this. I also tried the fact a symmetric matrix, $S$, can be transformed in $S = QDQ'$, for some orthogonal matrix $Q$ and diagonal matrix $D$. Can anyone help me with this solution?


title: I can't test a contrast in my repeated-measures ANOVA in Stata 12. Why?

text: The anova command is               anova score method/subject|method prediction method#prediction / prediction#subject|method infdir method#infdir / infdir#subject|method prediction#infdir method#prediction#infdir if(agegroup==0 & block==1), repeated(prediction infdir)      where `method` is a between-subjects factor, and each level of `infdir` is nested in `prediction` (repeated within subjects: level 1 and 2 of `infdir` are prompted for each level of `prediction` within subjects, where each of those is between subjects in `method`). I get the lovely ANOVA output I want that I've replicated in SPSS. However, I cannot use the `contrast` function at all. It, `margins`, and similar functions all return `(not testable)` whenever I attempt to examine simple effects. Does anyone know why? Edit: For some reason I can't answer my own question in such a way that it counts as answered so other people can reference it later. I'll fix it, but until then: I figured it out. I needed to add               emptycells(reweight)      After my `contrast` command. Empty cells in the interaction between the factors were throwing off the functions. I just had to RTFM more.


title: the approach for checking whether a process is stable?

text: I have two scenarios for time series data. 1) I have a uni-variate variable spanning across the time axis, are there any approaches or statistic to check whether this process is stable? 2) I have a multivariate variable spanning across the time axis, are there any approaches or statistics to check this multivariate variable is stable?


title: Determining the variables when training a/classifying with a random forest

text: Based on a response to the question "Statistical classification of text", I've decided that I want to use a random forest to perform various classifications of text. Specifically, I am getting YouTube video data from their XML feeds. Here is an example:               <entry gd:etag="W/"DUQBR347eCp7I2A9WhdXFEQ."">         <id>tag:youtube.com,2008:video:ltkMAcv2cas</id>         <published>2010-12-23T18:11:28.000Z</published>         <updated>2011-08-28T02:09:16.000Z</updated>         <category scheme="http://schemas.google.com/g/2005#kind"             term="http://gdata.youtube.com/schemas/2007#video"/>         <category scheme="http://gdata.youtube.com/schemas/2007/keywords.cat"              term="SSF4AE"/>         <!-- More categories/keywords here -->         <title>SSF4AE Mago 2D God VS Kazunoko aka Inoue</title>         <!-- More video information here, description, etc -->     </entry>      With that, I have very specific data that I believe I can use as the basis for variables in a random forest:               Variable    Value     --------    -----     Title       SSF4AE Mago 2D God VS Kazunoko aka Inoue     Category    SSF4AE      Looking at the video and based on the title/categories, I know that this is a video of the game Super Street Fighter 4: Arcade Edition (SSF4AE). Looking at the alglib Decision Forest, I also read the general principals of data analysis methods, which speaks about encoding specifically for the algorithms on alglib; while I won't use it specifically, I believe it expanded my understanding of how data should be represented for a random forest. Basically, using the above data, I'd need to transform each data point into a number of some sort. For example, I could transform this into the following variables:               Variable                 Value     --------                 -----     Has "SSF4AE" in Title        1     Equals "SSF4AE"              1      Now, assuming I know that this video is for SSF4AE, and that I know that SSF4AE is a common acronym used to identify the game, would using the above as variables work for a random forest for training and classification? Of course, that's a small set, and it's a very rough set of variables; my understanding is that a random forest basically will create nodes where there is a threshold for the variable, values over this fall on one side of the tree, values under on the other side. In the above example, there are only two options, true or false. However, I imagine I can do something like get the term-frequency/inverse document frequency and use that (using Lucene.NET or elasticsearch), like so (numbers are not accurate)               Variable                               Value     --------                               -----     Term Query for "SSF4AE" in Title        0.36     Term Query for "SSF4AE" in Category     1.00      And then train the random forest on that. I assume that if the variables were done as above, then I could also place things such as boosts on the variables (in the tf/idf calculation) where I see that they are more prominent or I feel I want them to have more weight. And finally, I can add more variables, but they have to be very specific. For example, I _know_ that if there is a character Yun that only appears in SSF4AE, (not previous versions of Street Fighter 4), as well as in Street Fighter III. To that end, I can add two more variables:               Variable                               Value     --------                               -----     Term Query for "SSF4AE" in Title        0.36     Term Query for "SSF4AE" in Category     1.00     One character is Yun                       1     Term Query for "3S" in Title            0.00     Term Query for "3SO" in Title           0.00      ("3SO" and "3S" being common acronyms for Street Fighter III: Third Strike) However, for another video, which I know is for 3S, I the variables might look like this:               Variable                               Value     --------                               -----     Term Query for "SSF4AE" in Title        0.00     Term Query for "SSF4AE" in Category     0.00     One character is Yun                       1     Term Query for "3S" in Title            0.54     Term Query for "3SO" in Title           0.00      Where the training set would indicate that video is for 3S. That said, are these feasible ways of going about determining the variables of a random forest? Also, if I understand correctly, the random forest will also indicate which variables contributed the most to the classification, so I can discard them during classification later (if I want) or add different ones to get better classifications from my training data.


title: Statistical Test against randomness

text: I have the following problem: vector = {1,0,0,1...}, only consists of 1 or 0, and size of vector is N. I am calculating proportion of 1 in such vector. Say, I got prop_1 = 0.96 Precisely, my null hypothesis is: H0: 0.96 is random HA: 0.96 is not random And, I would like to reject H0 to show that 0.96 is not random. What could be a procedure? Bootstrapping? Or some other ways?


title: Predicting from a simple linear model with lags in R

text: I have a dataset that I want to fit a simple linear model to, but I want to include the lag of the dependent variable as one of the regressors. Then I want to predict future values of this time series using forecasts I already have for the independent variables. The catch is: how do I incorporate the lag into my forecast? Here's an example:               #A function to calculate lags     lagmatrix <- function(x,max.lag){embed(c(rep(NA,max.lag),x),max.lag)}     lag <- function(x,lag) {      out<-lagmatrix(x,lag+1)[,lag]      return(out[1:length(out)-1])     }          y<-arima.sim(model=list(ar=c(.9)),n=1000) #Create AR(1) dependant variable     A<-rnorm(1000) #Create independant variables     B<-rnorm(1000)     C<-rnorm(1000)     Error<-rnorm(1000)     y<-y+.5*A+.2*B-.3*C+.1*Error #Add relationship to independant variables           #Fit linear model     lag1<-lag(y,1)     model<-lm(y~A+B+C+lag1)     summary(model)          #Forecast linear model     A<-rnorm(50) #Assume we know 50 future values of A, B, C     B<-rnorm(50)     C<-rnorm(50)     lag1<-  #################This is where I'm stuck##################          newdata<-as.data.frame(cbind(A,B,C,lag1))     predict.lm(model,newdata=newdata)


title: Predicting a maximum value with little data

text: My problem is i'm trying to figure out how many servers might be required to handle a theoretical maximal load of data requests. To do that I need to know what the maximum number of requests in a second might be. I have a fairly limited set of data that 156,000 transactions occurred in one hour. This is the maximum recorded transactions for a rolling 60 minute period over 4 months. 4500 devices have access to the server and each client-server transaction takes and average of 45 seconds. I was wanting to predict what the maximum number of client-server transactions in a single second might be with 95% certainty. Each server can service a maximum of 15 concurrent transactions before they start to experience lag. It's a long time since I've done any stats so if someone could point me in the right direction then that'd be brilliant. I figured I might need to know the max/min/average transaction time to be able to get to the variance but i'm probably just babbling now. Thanks a million.


title: What should the appropriate model for my data look like?

text: My research objective is that I would like to test if some sounds in languages are similar to each other in terms of duration. In my situation, I have 4 languages with 5 sounds. The first three languages have all 5 sounds but only one language has 3 sounds. My items are divided into 6 vowels and all languages have all these 6 vowels. Then I have speakers of 4 languages but 2 languages are produced by the same group of speakers - so actually I have 3 groups of speakers instead of 4 (each group has both males and females). As duration is continuous variable, I am using Linear mixed model (LMM) in R with lme4. My fixed effects are language + vowel + sex + sound. My random effects are speaker + item. However, I am not sure what the proper model should look like in my case. I have tried this formula: `lmer4 <\- lmer(Duration ~ (1|item) + (1+sound*vowel|speaker) + vowel*sex*sound*Language, data=data.frame,REML=FALSE)` but the model cannot be run because some sounds do not exist in some languages. Does anyone have any idea about the appropriate model for my case? Your help would be exactly grateful.


title: Kalman filter vs. smoothing splines

text: Q: For which data is it appropriate to use state-space modeling and Kalman filtering instead of smoothing splines and vice versa? Is there some equivalence relationship between the two? I'm trying to get some high-level understanding of how these methods fit together. I browsed through Johnstone's new Gaussian Estimation: Sequence and Multiresolution Models. It surprised be that there was not one mention of state-space models and Kalman filtering. Why wouldn't that be in there? Isn't that the most standard tool for these sort of problems? The focus, instead, was on smoothing splines and wavelet thresholding. I'm now very confused.


title: How do calculate Present Value in Google Spreadsheet?

text: Let say I have this kind of data:               1 01/1/1980     2 01/2/1999     3 03/12/2000     -1 03/6/2005     -5 07/07/2007      how can I calculate the Present Value (PV) to them in respect to current date, let say with 5% interest rate, in spreadsheet? `Current date` means today. **[Update]** I am probably misunderstanding terms `PV` and `FV`. I am trying to find some ready function similar to this one. Whatever method you use it must work with the above data. Please, stop spam.


title: Evaluate improvement in questionnaire questions

text: I asked two different groups of subjects to answer a questionnaire of 2 yes/no questions. The "correct" way of answering the questionnaire was to answer yes to both questions. I want to prove or disprove that group2 is answering better than group1. My idea is to check if a get a "low" p-value in a Chi-square independence test on a 2x2 contingency table (1 degrees of freedom) where rows are group1 and group2 and columns are the number of subjects that replied yes-yes and the number of subjects that replied differently. However, I am not sure this is the right way to proceed. Should I use a 4 columns contingency table (3 degree of freedom) dividing subjects in 4 categories (those that replied yes-yes, no-yes, no-no, and yes- no to the two questions) ? Thanks


title: Balls in urns of different sizes

text: I have 7 balls to be randomly placed in four urns of size {5,3,2,2}, what is the probability of exactly three urns containing at least one ball? (Phrased differently, twelve equally likely spots for the seven balls, but the first five is considered the first urn, next three the second urn, etc.) Through simulation got an answer of about 32.4%, but am wondering how to solve it in the general case? (n balls in urns of size {m_1,m_2,...,m_k}, probability of exactly r urns containing at least one ball?)


title: Calculating probability question from a student

text: I have a simple question for you guys but I am studying for exams and I honestly forgot how this was done. Here is the problem: > "A vase contains 4 white marbles, 3 red and 2 green. You draw marbles 2 > without putting the first one back after the second. > > What are the chances that you draw a red marble after the first one was > green? > > What are the chances that you draw 2 white marbles? > > What are the chances that you draw a white and a red marble?" I see I gotta calculate P(A ∩ B) / P(B) but how do I calculate P(A ∩ B) then? Also the third question from my problem is completely puzzling me even more than the others I don't have any answers and I have come up with like 100 different answers and I don't even know if any of them are right since I can't compare them to the actual answers and I don't know how to calculate it at all


title: Likelihood at MLE and transformations, the multivariate normal case

text: Given a univariate sample $\vec X = X_1, ..., X_n$ with standard deviation 1 and a strictly monotone transformation $t: R \to R$ with the property that the standard deviation of $t(\vec X)$ is also 1 (where $t(\vec X)$ is $t$ applied to each $X_i$). If one know fits a normal distribution to $\vec X$ and $t(\vec X)$, one observes the following fact: the likelihood at the respective MLEs of $\vec X$ and $t(\vec X)$ is the same. The reason for that is, that the standard deviation is the MLE for the $\sigma$ parameter of the normal distribution, so just the mean changes, which does not change the likelihood. In some computations with _multivariate_ data I observed the same fact, using transformations such that $\det(cov(t(\vec X)))=1)$. But I did not see a straightforward proof for that! Note that in the multivariate case, the MLEs for the variance parameters do change, even if the above determinant is always 1. I am sure I am not the first to notice, but my (somewhat limited) literature on the multivariate normal distribution does not give me a clue how to prove the above statement, i.e. that the multivariate normal likelihood does only depend on the determinant of the $cov(t(\vec X))$, regardless of choice of $t$. A proof boils down to looking at the terms of the form $x'\Sigma x$ for $\det \Sigma = 1$. Do you have more clue? Thanks, Philipp


title: Randomization Testing help

text: I am scratching my head over this one...any help would be greatly appreciated. I want to know if the average travel time between Guell Park and the beach in Barcelona using the bus or the metro differ and which one is faster. Since I do not know the distribution of travel time, I can't use a t-test; therefore, I was suggested to use a randomization test (so I don't make any assumptions about the data). Here are the data: Metro (in min): 34, 39, 28, 42, 32, 37, 41, 40, 47, 36, 38, 48, 31 Bus (in min): 40, 44, 49, 50, 48, 45, 47, 46, 39, 55 How would I go about proving there is or is not a different between the forms of transportation (randomization tests)?


title: Two dice problem

text: Two players are throwing each one die. The one that has higher value, receives a number of points equal to the difference of the values on both dice. How do you estimate probability for a winner to collect more than 100 points in 100 throws?


title: How to re-cluster new instance in centroid base clustering?

text: I have applied clustering algorithms like k-mean, k-medoid and DBSCAN on my patients dataset. For each algorithm RapidMiner generated a clustered model (centroid table and graphs etc) and a clustered set (shows which examples are part of which cluster). Now I want some way that when a new patient come i want to assign him a cluster based on previous trained model. I am confused of the way to do this… Is it something like that (I may be wrong)?   1. For each attribute value of new patient - that attribute value from centroid table summing all the differences of attributes of patient and taking average.   2. Then assign him cluster whose average is minimum with respect to that patient. If this the right way then how will I re-cluster, i.e when a new patient comes our algorithm is assigning him cluster, that's mean. centroid moves and then I have to re-cluster with each record insertion. How to handle this in my scenario?


title: Is a Gaussian AR process with white noise independent?

text: I was just wandering if, given the AR process \begin{equation} X_t = \alpha X_{t-1} + \varepsilon_t, \quad \varepsilon_t \overset{iid}{\sim} N(0,1), \end{equation} the $X_t$ values are independent due to the independent white noise?


title: normalization in max-sum algorithm (loopy belief propagation)

text: I was implementing the max sum algorithm for a general graph (i.e., the ones with a cycle). I updated the messages as indicated in http://www.cedar.buffalo.edu/~srihari/CSE574/Chap8/Ch8-GraphicalModelInference/Ch8.3.3-Max- SumAlg.pdf. What I found is that, however, in the case of a factor graph with a cycle, if the factor value is between 0 and 1, the messages will keep decreasing as the iteration goes. It is because the log of the factor value is negative and the looping in the cycle will keep adding the (negative) factor value. This can be easily verified with a factor graph with two random variable nodes and two factor nodes connecting the two random variable nodes. It seems to me that there should be a kind of normalization after each message update but the slides I am refering to do not mention it. Does anybody know the solution for this phenomenon?


title: General questions regarding text-classification

text: I'm new to Topic Models, Classification, etc… now I'm already a while doing a project and read a lot of research papers. My dataset consists out of short messages that are human-labeled. This is what I have come up with so far:   * Since my data is short, I read about Latent Dirichlet Allocation (and all it's variants) that is useful to detect latent words in a document.   * Based on this I found a Java implementation of JGibbLDA but since my data is labeled, there is an improvement of this called JGibbLabeledLDA   * In most of the research papers, I read good reviews about Weka so I messed around with this on my dataset   * However, again, my dataset is labeled and therefore I found an extension of Weka called Meka that had implementations for Multi-labeled data   * Reading about multi-labeled data, I know most used approaches such as one-vs-all and chain classifiers... Now the reason me being here is because I hope to have an answer to following questions:   1. Is LDA a good approach for my problem?   2. Should LDA be used together with a classifier (NB, SVM, Binary Relevance, Logistic Regression, …) or is LDA 'enough' to function as a classifier/estimator for new, unseen data?   3. How do I need to interpret the output coming from JGibbLDA / JGibbLabeledLDA. How do I get from these files to something which tells me what words/labels are assigned to the WHOLE message (not just to each word)   4. How can I use Weka/Meka do get to what I want in previous question (in case LDA is not what I'm looking for) I hope someone, or more than one person, can help me figure out how I need to do this. The general idea of all is not the issue here, I just don't know how to go from literature to practice. Most of the papers don't give enough description of how they perform their experiments OR are too technical for my background about the topics.


title: What is meant by the standard error of a maximum likelihood estimate?

text: I'm a mathematician self-studying statistics and struggling especially with the language. In the book I'm using, there is the following problem: A random variable $X$ is given as $\text{Pareto}(\alpha,60)$-distributed with $\alpha>0$. (Of course, you could take any distribution depending on one parameter for the sake of this question.) Then a sample of five values $14$, $21$, $6$, $32$, $2$ is given. First part: "Using the method of maximum likelihood, find an estimate $\hat{\alpha}$ of $\alpha$ based on [the sample]." This was no problem. The answer is $\hat{\alpha}\approx 4.6931$. But then: "Give an estimate for the standard error of $\hat{\alpha}$." What is meant by this? Since $\hat{\alpha}$ is just a fixed real number, I don't see in what way it could have a standard error. Am I to determine the standard deviation of $\text{Pareto}(\hat{\alpha},60)$? If you think the question is not clear, this information would help me as well.


title: Can univariate linear regression be used to identify useful variables for a subsequent multiple logistic regression?

text: Does the $R^2$ (or some other statistic) from a univariate linear regression tell me anything about how it would work in a logistic model? What if I normalized the data to mean zero? I'm doing multiple logistic regression, and I'm wondering if I created a ton of univariate linear regressions, if I can use the $R^2$ as an indicator of variables that are worth trying or not. Perhaps even adding a brute force search of interesting interactions or adding polynomials, etc.


title: Variable is significant through stepwise regression but not in final model's summary; which should I report?

text: I used generalized linear mixed models (with the `glmmADMB` package) to identify environmental factors related to parasite abundance in rodents. I used stepwise backward elimination to sequentially simplify the full model until only significant factors and interaction terms remained (Crawley, 2007). However, when I checked the final model's summary, one of the factors is no longer significant (`p = 0.087`). However, in the ANOVA tests that led to the final model, this variable (temperature) is significant: **Final model's summary:**               Coefficients:                                Estimate Std. Error z value Pr(>|z|)                     (Intercept)     3.131658   0.213630   15.33  < 2e-16                 temperature    -0.004019  0.0033143   -1.71  0.08673 .      **anova**               Model 1: parasites ~ age + vegetation + synchrony            Model 2: parasites ~ temperature + age + vegetation + synchrony            NoPar  LogLik Df Deviance Pr(>Chi)             1     7 -550.75                                  2     8 -546.81  1    7.898 **0.004949**          My questions are which significance value should I report? Does this mean that temperature is not important? I've seen manuscripts that report values from final's model summary but also some that report significance values from model comparisons, which one is correct?


title: Correlations circles: Interpretation

text: I'm trying to interpretate some **circles of correlations** in the frame of a **Canonical Analysis of Correlations**. I have been told that only points in the halo could be interpretated with those **following rules** :   1. If their angle is close to $0$ they are positevely correlated.   2. If their angle is close to $\pi$, their are negatively correlated.   3. If their angle is close to $\frac{\pi}{2}$ they are independent.   4. If they are close to the extern edge of the circle, they are higly correlated with the basis you choose. But is it possible to say something regarding the **points near the center** , in the small circle or not? I just know that they are not correlated to the basis, so you could change the dimension to choose. But are there other ways to interpretate them or are other rules? Here is an example of what I have: ![enter image description here](http://i.stack.imgur.com/qKIDG.jpg) Thanks for helping.


title: Confidence interval of metric from empirical distribution

text: I'm familiar with using `boot` in R to estimate a confidence interval for VaR. For example, I have a set of monthly returns for an asset, and I want a range on the 95% VaR. Now, however, I'm given an empirical PDF, and I want to estimate the CI for VaR. I have no information on the distribution other than this PDF. It is an estimate of losses due to hurricanes – it generally has positive skew, $\ln$ or beta in shape, but certainly not analytic. A brute-force idea would be to generate a large number of random numbers in [0,1], invert from the empirical CDF, and get the 95% of the inverted values. Then repeat 1000 times to get a distribution of the 95%ile value. `boot` doesn't seem to apply here. Is there a more elegant solution? Thanks.


title: What test should I use to see correlations between overlapping groups and a score?

text: To put it simply, I have the courses students have taken and scores on an exam. The students come from different course backgrounds. (Some have taken only course A some only D some A and D some A and B some B C and D etc etc) What can test can I use to account for this?


title: How to do feature selection for learning from positive and unlabeled examples?

text: I have a binary classification task for German webpages for which I only have positive examples. That is why I use learning from positive and unlabeled examples as described on this page, also known as _partially supervised learning._ At the moment, I'm just excluding very rare features which occur only once and very frequent features, i.e., stop words and those which occur in more than 50% of all positive examples. As for classifiers, I want to test Naive Bayes (the positive example webpages are quite short, so I favor the multivariate version with boolean features) and Support Vector Machines (SVM). I've read that feature selection is not so important for SVM as it does not effect the classification results very much. Is that true? As I'm not so familiar with feature selection algorithms, can you recommend an algorithm that works especially well with features from positive examples only and which generally yields better results than just cutting off very rare and very frequent features? If it is not possible to give a general answer to this question and if it highly depends on my data set, then please say so as well. Thanks a lot!


title: Quantile regression with dummy variable that's equal to 0 over most $t$

text: I have the model $Y_t = a + b*X_t + c*D_t + e_t$, where $t \in T = \{1,...,3000\}$ and $D_t$ is a binary variable equal to $0$ over $T \backslash \{20,21,...,30\}$, and equal to $1$ over $\{20,21,...,30\}$. $c$ is highly significant in OLS. However if I run quantile regression $Q_{Y_t}(\tau | X_t,D_t) = a(\tau) + b(\tau)*X_t + c(\tau)*D_t + e_t$ then $c(\tau)$ is insignificant for any $\tau \in [0,0.4]\cup[0.6,1]$. Now what I want to know is if QR could give false tail results for a dummy variable that's equal to unity over a tiny portion of $T$? As in, could the population relationship be strong, but QR isn't "strong" enough to pick up any tail relationship when $D_t = 0$ over almost all $t \in T$?


title: How to derive MA(1) from Wold decomposition

text: I want to derive the moving average process 1 from the Wold decomposition theorem. I don't know how to define the coefficient of MA(1), though it seems that it should be easy.


title: $E(x^k)$ under a Gaussian

text: What would be the expectation of $|x|^k$, where $x\sim\mathcal{N}(0,1)$, $k>0$ and $k$ is not an integer?


title: Granger Causality Testing With Panel Data

text: I'm trying to apply a _Granger Causality test_ to panel data. I've found enough literature to understand that topic. However, I've been unable to find and **R** package to carry out that analysis. I'm wondering if anybody know whether there is around any package to deal with that. Thanks! ### _I'm adding a potential solution, but new questions arose._ The solution that I found is apply a _Granger Non- Causality test_ and using _Generalized Method of Moments_ ( _GMM_ ). In the Erdil & Yetkiner’s (2004) paper you can find a description of _Granger non-causality test_ with panel data. To perform a _GMM_ I used the `plm` package for **R**. If you have a look at its tutorial (Croissant & Millo, 2008), you’ll see that the built-in function `pgmm` (page 17) removes the individual effect by the first difference and time dummies are included. The function’s summary also provides some tests to assess the model. For instance, to check serial autocorrelation in the residuals, Wald tests for coefficients and for time dummies and the Sargan test to evaluate if there is correlation between the instrumental variable and the residuals. Then I performed a Wald test (the first one in Erdil & Yetkiner, 2004) with the sum of squared residuals of an unrestricted model ( _SSRu_ ) and of a restricted model ( _SSRr_ ). Now, my questions for the audience are: 1) Do the time dummies remove the time effect? I think so. 1.1) What if the time dummies aren't significant? 2) Therefore, if I got rid of the individuals and time effects, is the Wald test ( _SSRr-SSRu_ ) be as a Wald test applied to an `OLS` model? I think so. If so, I’m not sure about the freedom degrees. Let’s see first the test suggested by Erdil & Yetkiner (2004): $$W=\frac{(SSRr-SSRu)/Np}{SSRu/[NT-N(1+p)-p]}$$ where _N_ = number of individuals, _T_ =years and _p_ =number of lags. Note that they didn’t get rid of individuals and time effects (at least that's what I understood). Now, if I got rid of individuals and time effects the Wald test as applied to OLS models would be: $$W=\frac{(SSRr-SSRu)/m}{SSRu/ (n-k)}$$ where _m_ = number of restrictions (number of coefficients that were removed from the unrestricted model to turn it restricted), _k_ = total number of coefficients in the unrestricted model and _n_ = number of observation. More questions: 3) What is number of observation? 3.1) Is it the number of year or number of years*number of individuals? If it is number of years it seems reasonable, but if it is the product between years and individual it doesn’t. For example, in my case I have 328 individual and 13 years, so it is 4264; therefore, the numerator in the Wald test will be very, very small and I’ll be rejecting everything. Finally, 4) Am I right doing as I did? Again, any help will be much appreciated


title: Fitting data based on an unobserved variable

text: I have pairs (x, y) which I would like to regress (x independent, y dependent). Plotting them, I see distinct bands which can be attributed to a third variable. Unfortunately this variable cannot be directly measured. Is there a good way to both fit the data to an equation and determine the value of the binary variable (for each x)? The equation I was considering is $$y=ax^b(1+cz)$$ where $z$ is the unobserved binary variable. (Of course there are two ways to fit the data depending on which group is chosen as 0 and which as 1, but restricting $c>0$ should make the choice essentially unique.) Unless there's reason to prefer something else I would look for a least- squares fit.


title: Error with caret, using "out-of-bag" re-sampling

text: When I run this code (R 2.13.2, caret 4.98 on windows 7):               library(caret)     data(iris)     TrainData <- iris[,1:4]     TrainClasses <- iris[,5]     rfFit <- train(TrainData, TrainClasses,                      method = "rf",                      tuneLength = 3,                      trControl = trainControl(method = "oob"))      I get the following error:               Fitting: mtry=2      Fitting: mtry=3      Fitting: mtry=4      Aggregating results     Selecting tuning parameters     Fitting model on full training set     Error in structure(list(method = method, modelType = modelType, results = performance,  :        object 'resampledCM' not found      What am I doing wrong? /edit: I updated to caret 4.99 using R forge, and the problem goes away. It seems to have been unique to caret 4.98 and R 2.13.2.


title: R DFA eigenvalues

text: I'm doing some followup to MANOVA analysis and I was wondering how in R you get eigenvalues for DFA model? There is a attribute `$svd` in the return value of `lda` function, but there is a conflicting information online regarding whether those are actually eigenvalues or not. I'm reading _"Discivering statistics with R"_ by Field et al. and there is written that _"proportion of trace are eigenvalues expressed as proportion"_. I've looked how proportion of trace is calculated and this is the part in `MASS:::print.lda` where they are calculated               if (length(svd) > 1L) {         cat("\nProportion of trace:\n")         print(round(svd^2/sum(svd^2), 4L), ...)     }      From this seems to me that eigenvalues are calculated as `svd^2`. However online I've found several examples (e.g., example1 and example2) where people are using values of `svd` as eigenvalues. Help page for `lda` function it says that `svd` attribute is a list of > the singular values, which give the ratio of the between- and within-group > standard deviations on the linear discriminant variables. Their squares are > the canonical F-statistics. However, I'm not sure what is correct, and I would like to report those in some analysis I'm conducting and don't want to put wrong values.


title: Finding number of gaussians in a finite mixture with Wilks' theorem?

text: Assume I have a set of independent, identically distributed univariate observations $x$ and two hypotheses about how $x$ was generated: $H_0$: $x$ is drawn from a single Gaussian distribution with unknown mean and variance. $H_A$: $x$ is drawn from a mixture of two Gaussians with unknown mean, variance and mixing coefficient. If I understand correctly, these are nested models since the model that $H_0$ represents can be described in terms of $H_A$ if you constrain the parameters of the two Gaussians to be identical or constrain the mixing coefficient to be zero for one of the two Gaussians. Therefore, it seems like you should be able to use the E-M algorithm to estimate the parameters of $H_A$ and then use Wilks' Theorem to determine whether the likelihood of the data under $H_A$ is significantly greater than that under $H_0$. There's a small leap of faith in the assumption that the E-M algorithm will converge to the maximum likelihood here, but it's one I'm willing to make. I tried this in a monte carlo simulation, assuming that $H_A$ has 3 more degrees of freedom than $H_0$ (the mean and variance for the second Gaussian and the mixing parameter). When I simulated data from $H_0$, I got a P-value distribution that was substantially non-uniform and enriched for small P-values. (If E-M weren't converging to the true maximum likelihood, the exact opposite would be expected.) What's wrong with my application of Wilks' theorem that's creating this bias?


title: Analysis of interaction tables

text: I have two tables of interactions. One represents the number of times that a clownfish of a given species is found in an anemone of a given species. One represents the number of times that a clownfish of a given species is found in the same anemonae than a clownfish of another species.               Table_anemone_fish                Fish_a    Fish_b   Fish_c     An_A    23        56       12     An_B    12        5        5     An_C    23        10       68      Note: Summing this table gives the total number of anemones in my sample.               Table_fish_fish                  Fish_a    Fish_b   Fish_c     Fish_a   NA         7        3     Fish_b    7        NA        1     Fish_c    3         1        NA      Note: Most of the times the fishes live with no other species and they are therefore not represented in this table. The frequency of each fish species can be obtain with `Table_anemone_fish` The questions I want to answer are:   * Is there differential preference of fish species in choosing their anemones?   * Is there different preference of fish species in choosing with which other species to live with? To answer this post one can simply gives the kind of test to use and their philosophy or even better, one might give an example of how to analyze these data with R.


title: Poisson process - calls arriving

text: Already posted on MSE. Had no answer, so will post here. Assume the number of calls per hour arriving at an answering service follows a Poisson process with $\lambda = 4$. Question: If it is know that $8$ calls came in the first two hours. What is the probability that exactly five arrived in the first hour? Attempt: Isn't this just a combinatorial question? So the answer is ${8 \choose 5}/2^8$


title: How to to calculate the topic distribution of a document

text: I have a simple (may be stupid) question. I want to calculate Kullback–Leibler divergence on two documents. It requires probability distribution of each document. I do not know how to calculate probability for each document. Any simple answer with layman example would be much appreciated. Let's say we have follow two documents:               1 - cross validated answers are good      2 - simply validated answers are nice      (wording of the documents is just bla bla to give you an example) How do we calculate probabilities for these documents? Let's say we add one more document:               3 - simply cross is not good answer      If we add another document then how would it impact probability distribution? Thanks


title: Systematic way of choosing a threshold?

text: I am performing some analysis on understanding various characteristics of a certain population. To narrow down my selection (as I will be going through them manually), I am using different thresholds on certain measurable properties (which tell me how relevant someone is) to draw the following frequency graph. For instance, the line marked `Threshold 1` tells me that when I fix my threshold at `t1`, there are approximately 360 people with at least 1 friend, 290 people with at least 2 friends etc. I am aware that I will be having a hard time explaining what I am doing but I will try my best. I am trying to fix two thresholds: one that is used to obtain the line and the other that I use to filter noise. The '0' on the graph marks the thresholds. In this particular case, I am using `Threshold 3` and filtering my population to contain only those people that have `>= 2 friends`. While this is working well for me, I am trying to determine if there is systematic procedure to explain why these thresholds are better than others. For instance, `Threshold 2` is not such a good filter because it behaves like `Threshold 1`. Utilizing `Threshold 4` is not so useful because the returns in going from `Threshold 2` to `Threshold 3` is much higher than that going from `Threshold 3` to `Threshold 4`. Any suggestions on where I should start? Basically, I want to do something similar to cost-benefit analysis (that's the closest I could think off). ![enter image description here](http://i.stack.imgur.com/ehzXe.png)


title: What are the measure for accuracy of multilabel data?

text: Consider a scenario where you are provided with KnownLabel Matrix and PredictedLabel matrix. I would like to measure the goodness of the PredictedLabel matrix against the KnownLabel Matrix. But the challenge here is that KnownLabel Matrix have few rows only one 1 and other few rows have many 1's(those instance are multi labeled). An example of KnownLabel Matrix is given below.               A =[1 0 0 0         0 1 0 0         0 1 1 0         0 0 1 1         0 1 1 1]      In the above matrix, data instance 1 and 2 are single label data, data instance 3 and 4 are two label data and data instance 5 is the three label data. Now i have PredictedLabel Matrix of data instance using an algorithm. I would like to know various measure which can be used to measure the goodness of the PredictedLabel Matrix against KnownLabel Matrix. I can think of frobeinus norm difference between them as one of the measure. But im looking for the measure such as $accuracy (= \frac{Correctly\_predicted\_instance}{total\_instance})$ Here how do can we define the $Correctly\_predicted$ for multiple data instance?


title: Performance testing: detecting changes in distributions over time

text: I have an interesting application I’d like some advice on. A common task when working with stochastic systems in an engineering context is testing for regressions and improvements in performance. First, it’s important to have a confidence one is collecting enough data points to have an accurate assessment of the distribution of performance. And second, it’s important to have a good strategy for comparing those distributions to identify any discernible performance changes. However, I’m not sure there’s a widely accepted solution for how to tackle these challenges. Or if there is, I don’t know about it yet. I’ll describe the problem I’m facing in more detail. I have a process whose objective is to reach a certain goal in as little time as possible. The process is stochastic, so the runtime has some random variations. In fact, the variance of the runtime could be fairly large. Every day I collect one “sample” by running the process to completion for N repetitions and recording the elapsed time per repetition, thus collecting N data points per sample. My recording precision is 0.1 seconds. After I’ve collected a sample I’d like to compare it with yesterday’s sample to measure if there’s been a change in performance. What do I mean when I say "a change in performance?" I'm interested in significant changes to the average runtime (the mean) and to the variation in the runtimes (variance). Part of the issue is I'm not quite sure how to define "significance" in this context. My definition for now is any change in the mean greater than four times the recording precision, and any change in the variance greater than a few percent. There's an additional complication: I'd like to find a method where I don't need to make assumptions about special features of the distribution, like the specific scale or shape. The reason for this is I'm performing the process I described above for several stochastic processes which each have distinct properties. For example, one produces a bimodal distribution; another produces a one-tailed distribution. I had tried setting up problem-specific thresholds for the mean and variance, but tuning them proved quite difficult, especially since each has high variance. Given the problem I’ve just described, I have two questions: **1) How do I know when I’ve collected enough data points per sample to make an accurate comparison of two sample distributions?** I would love to collect millions of datapoints per sample but it takes several hours to collect 1000 datapoints, so I’m dealing with a classic tradeoff. **2) What are some traditionally used methods for comparing two samples in order to detect any significant change in the mean and variance of the stochastic process?** This is assuming I’ve collected enough data points to have fairly accurate views of the distributions of the two samples. I would include the work I’ve done so far to tackle this problem, but I don’t want this post to be any longer. My knowledge of statistics is fairly basic, which is why I’m appealing to the good folks at crossvalidated for advice or pointers to related materials. I've found a few related posts but nothing which described my problem exactly. Thank you for reading and I appreciate anything you can give me.


title: What distribution function might describe the following statistical experiments?

text: I have made several experiments, the result graphs are below, I need to figure out what is the distribution function of the following random variables... Where do I start with statistical analysis of my experiment? How to perform these analysis with spss? ![Sequence 1](http://i.stack.imgur.com/arAJ4.png) ![Sequence 2](http://i.stack.imgur.com/CYwAD.png) **EDIT** Explanations: I'm testing algorithm efficiency for 2 algorithms, the graphs describe efficiency distribution.   1. I made 500000 tests for 2 algorithms, which produced 2 sets of 500000 floating point values.   2. I've computed MIN and MAX value of the set (The smaller the values the better the result is, the bigger the values the worser the result is)   3. I wanted to see the distribution of values according to 500 groups, so I've computed step = (MAX-MIN) / 500   4. I've computed the amount of values that fall into each of 500 segments [MIN + 0 * Step..MIN + 1 * Step], [MIN + 1 * Step, MIN + 2 * Step], ... [MIN + 498 * Step..MIN + 499 * Step] The graph you see below represents on x-axis an index of the segment, on y-axis the amount values fallen into the segment divided by the total number of values, and also divided by the `Step = (Max - Min) / 500`. As a result this should represent the distribution function of my algorithm efficiency... **EDIT** Here is DOC containing SPSS > Analyze > Descriptive Statistics > Explore output. Can anybody help reading that?


title: Gaussian Regression With Multiple Inputs?

text: Is it possible to use a Gaussian Process to relate multiple independent input variables (X1, X2, X3) to an output variable (Y)? More specifically, I would like to produce a regression graph like the example shown below where confidence interval reduces around clusters of data (i.e. variance is high at x = 1 where there is no data, but x = 0.3 the regression is tight due to the clustering of input variables) and Instead of having one input variable on the x-axis, where would be multiple inputs. ![enter image description here](http://i.stack.imgur.com/4YEAb.png) For example, is it possible to develop a regression relationship that relates the price of houses (HPrice = [125000, 63000, 500000]) to Floor Area (FArea = [856,497,1300]) and the Number of Bedrooms (BedR = [2,2,4])? Ideally I would like to do this in R, and wonder if there are any recommendations/example available? Thanks!


title: How to deal with negative coefficients in logistic regression

text: I have 5 single factor regression models. Each has intuitive (positive) coefficients which is good. If I run multiple logistic regression one of the coefficients turns negative. If I run stepwise selection, the negative coefficient stays in.. it's not redundant.. How should I deal with this? Is there a way to enforce positive coefficients? Thanks


title: Multilevel regression using lmer function in R and Stata

text: I have a large dataset and have performed a multilevel regression in Stata, the model is the following: `xtmixed dependent independen1 independent2 independent3 independent4 || independendt5:` So there is one grouping factor: `independent5` In R I did the following:               lmer(dependent ~ independent1 + indepdendent2 + independent3 + independent4 + 1 | independent5, REML=TRUE)       A few questions: is this identical? Stata output gives me number of groups, 100, while R gives number of groups as 99. Furthermore, the variances and standard deviations are not the same. Also I would like to know how to obtain p values and coefficients from the R ouput. I have done `fixef(model)` and `ranef(model)` but when I do `coef(model)` it says: `Error in coef(model) : unable to align random and fixed effects` Also Stata only gives me 1 coefficient for each predictor while the R `fixef(model)` gives one for each group. So someone familiar with both Stata and R could help me with this.


title: What statistical methods should I use for such case?

text: I have a set of 20 chemical indicators measured in soil samples collected from one (the same) place over some time. What i would like to find out is: - To find which indicators are somehow correlated with others. I don't want to assume dependent / independent variables, just to see if some parameters seem to be correlated somehow. What methods should I use? I thought of cluster analysis - is that right choice? What other options I do have.


title: ANOVA terminology: "repeated measures" vs . "within/between subjects"

text: What does the term "Repeated" mean in the context of ANOVA? I totally understand the distinction between within-subjects test and in-between subjects test, but am not familiar with the wording "repeated" in the following context: "ANOVA:Repeated Measures - between factors" and "ANOVA: Repeated Measures - within factors" (in G-Power 3).


title: Conditions for positive dependence

text: Consider two random variables $X$ and $Y$ with joint distribution $F_{X,Y}$ and strictly positive density function $f_{X,Y}$. Additionally, let $x^*$ be the value of $x$ that solves: $$ \Pr[Y<y\mid X=x^*]=1-k $$ for some constant $k \in [0,1]$. I would like to know the minimum conditions I should impose on $F$ in order to ensure that $$\frac{\partial x^*}{\partial k}>0$$ It seems quite intuitive that $COV[X,Y]>0$ would do the job, but this is not the case. I found in the literature that Affiliation* or Decreasing Inverse Hazard Rate** are sufficient, but I would like to know if there are weaker conditions. *The pdf $f$ is afﬁliated if $x \leq x$ and $y \leq y'$ imply $f(x, y')f(x',y)\leq f(x', y')f(x,y)$ ** $Y$ is Inverse Hazard Rate Decreasing in $X$ if $\frac{F(y|x)}{f(y|x)}$ is non-increasing in $x$ for all $y$ (where $f(y|x)$ denotes the pdf of $Y$ conditional on $X=x$). NOTE: This question has also been posted on the Stackexchange Mathematics site (this is the link)


title: Prioritizing data collection

text: In toxicology machine learning methods are used to estimate the probability of compound toxicity. Unfortunately most toxicology datasets contain on the order of 100 compounds and sometimes contain many more variables than compounds. When organizing the creation of new data one might ask how we can use a generative model to determine which chemicals would be valuable to test next, and which tests to perform on those chemicals. The goal of such data collection would be to improve the underlying model as quickly as possible. After doing some searching I have not found any papers approaching this question. Is anybody aware of methods used to organize data collection to improve the predictivity of a machine learning model as quickly as possible?


title: Choosing elements of an array based on given probabilities for each element

text: I'm not well-versed in statistics at all. I'm currently working on an implementing an optimization algorithm, and am having a bit of trouble on this part. I have an array of $n$ elements $A = \left\{a_{i}\right\}$ where $a_{i}\in\mathbb{R}$ and $1 \le i \le n$. Furthermore, I have an array of probabilities given by $P = \left\{p_{i}\right\}$ with $p_{i}\in[0,1]$. (Also, $\sum_{i} p_{i} = 1$ of course) If it matters, the probability distribution is triangular with $i=1$ having the highest probability and $i=n$ having the smallest probability. Now I want to choose $q < n$ distinct elements from $A$ using the probabilities of selection from $P$. I'm currently coding the algorithm in MATLAB, but I'm not sure how to go about choosing the $q$ elements.