title: Binomial GLM - Predicting the time of buying of a product

text: I'm trying to predict the frequency of buying a product based on data I have. I have two things to consider: 1- Predicting the success of buying this product based on frequencies (0: Failed ; >0: success) 2- Predicting the time of buying this product I could build a model to predict the success of buying this product, but I need to calculate also the time (how successful is the Hour/Day/Month/Year) ? If I can calculate also the time with the same model, I would have two variables (Predicting_the_success,Predicting_the_time) I am wondering now, I can't apply this formula as I am not sure if the both variables are not independent. > Total_Proba_Success = Predicting_the_success * Predicting_the_time A professor of mine said that I need to think about using regression here, What do you suggest as a model to solve this problem ? **UPDATE 3** I have predicted the success of buying this product based on frequencies (0: Failed ; >0: success). Based on comment of Gavin: I should turn this predicted values to discrete prediction say >= 0.5 == Buy, < 0.5 = Let's say that I have proba P : calculate discrete prediction of buying product success. Now I need to calculate based on that P, the probability of buying an item based on hour/day/month. Gavin suggested to use Binomial GLM: "binomial GLM, where the response is a vector of 0s and 1s (indicating product not bought and product bought respectively) and then vectors of the hour of day, day of week, and month of year, for each observation you made. The estimated coefficients from such a model will tell you the effect in terms of log odds (if a logit link is used) for each of those three predictors" I am confused! - I found this challenging, in my dataset, time is mentioned only in items where product bought. So, I ended-up using only time of products that have been bought, while in response, he mentioned that I should use a vector of 0s and 1s (indicating product not bought and product bought respectively).   * Should I calculate all the hours between two dates ? and include it in my model.  let's say for eg: 01-01-2014 -> 01-02-2014.               -+--------+-----------+--------+-------------------------+        hour       day       Month     number products bought       -+--------+-----------+--------+-------------------------+        12:00     01           01             2      In this case: should I transform my model to include all the hours that haven't been mentioned with frequencies 0. I mean, if we only use time of items that have been bought, we should consider another response vector ? Last thing, why do we use GLM Binomial instead of using Time-series prediction models? is it because we have "discrete times" and not "continuous times".               -+--------+-----------+--------+-------------------------+        hour       day       Month     number products bought       -+--------+-----------+--------+-------------------------+        00:00     01           01             00        01:00     01           01             00        02:00     01           01             00        03:00     01           01             00        04:00     01           01             00        05:00     01           01             00        06:00     01           01             00        07:00     01           01             00        08:00     01           01             00        09:00     01           01             00        10:00     01           01             00        11:00     01           01             00        12:00     01           01           **02**        13:00     01           01             00        14:00     01           01             00        15:00     01           01             00        16:00     01           01             00        17:00     01           01             00        18:00     01           01             00        19:00     01           01             00        20:00     01           01             00        21:00     01           01             00        22:00     01           01             00        23:00     01           01             00      **UPDATE 4** Based on Gavin's comment: "model the number of products bought as the response with hour, day, and month as predictors in a Poisson GLM" Just to confirm:   * Should we model the problem separately ?  For hours:                  Predictor      response      -+------------+-------------------------+        hour          number products bought       -+------------+-------------------------+        00:00               01        01:00               00      For days:                 Predictor      response      -+------------+-------------------------+        day         number products bought       -+------------+-------------------------+        00               03        01               01      For months:                 Predictor      response      -+------------+-------------------------+        month         number products bought       -+------------+-------------------------+        00               20        01               12      or should I do it this way                   Predictors                 +  response     -+--------+-----------+--------+-------------------------+        hour       day       Month     number products bought       -+--------+-----------+--------+-------------------------+        00:00     01           01             00        01:00     02           01             02      I have data of two years; if I do it this way hour-day-month as predictor, I will end-up having 2 values for everyone. which is not enough to predict -based on what I've learnt from other comments, the more data we have, the most accurate-. Another comment : Should I consider holidays. if I model it by Hour-day --> response Thursday 14:00 ==> 12 : 'cause it's Thursday 1st may (most workers are in holidays) Thursday 14:00 ==> 01 : 'cause it's Thursday 8th may (most of ppl working at this time ) After thinking about it; I am seeing that modeling the problem can make sense this way : Option1 - GLM Poisson :                   Predictors                 +  response     -+--------+-----------+--------+------------+-------------------------+        hour       day       month    isHolidays     number products bought       -+--------+-----------+--------+------------+-------------------------+        00:00     01           00       00          00        01:00     02           01       01          02      isHolidays: based on national holidays.   * Strength : It gives us all the details   * Weaknesses : we have a lot of cases, and for ( hour-day-month) we have 3 cases for everyone as we have data of three years. Option2 - GLM Poisson : Assuming that months doesn't make difference                   Predictors                 +  response     -+--------+----------+-------------------------+        hour       day        number products bought       -+--------+----------+-------------------------+        00:00     01           01              20:00     06           21              * Strength : we have many samples of data for every case. We don't need info about months as product is not seasonly (people need it during the whole year).   * Weaknesses : less details, Assuming that months doesn't make difference is little bit wrong, as we have people who are traveling during summer.  Option3 - GLM Binomial :               -+--------+-----------+--------+-------------------------+        hour       day       Month     number products bought       -+--------+-----------+--------+-------------------------+        00:00     01           01             00        01:00     01           01             00        02:00     01           01             00        03:00     01           01             00        04:00     01           01             00      I should add all the -possible- times whether product have been bough or no. then select time where product has been bough by using 1, 0 for no. Now, what is the most logical model here. They all have weaknesses, what do you recommend ?


title: What approach could be used for modelling virus infections?

text: I'm a physics graduate who ended up doing infosec so most of the statistics I ever learned is useful for thermodynamics. I'm currently trying to think of a model for working out how many of a population of computers are infected with viruses, though I assume the maths works out the same way for real-world diseases so references in or answers relevant to that field would be welcome too. Here's what I've come up with so far:   * assume I know the total population of computers, N.   * I know the fraction D of computers that have virus-detection software (i.e. the amount of the population that is being screened)   * I know the fraction I of computers that have detection software _that has reported an infection_   * I don't know, but can find out or estimate, the probability of Type I and II errors in the detection software.   * I don't (yet) care about the time evolution of the population. So where do I go from here? Would you model infection as a binomial distribution with probability like (I given D), or as a Poisson? Or is the distribution different?


title: What is the probability that a person will die on their birthday?

text: I am curious about what the probability is that a person will die on their birthday? I am sure there are a number of ways to approach this, plus I have heard that actual numbers point to a higher rate on birthdays, hence why I am asking it here.


title: Programs and procedure for Blackwood test?

text: I am looking for programs and procedure for the Bradley-Blackwood test, which is normally used to assess test-retest reproducibility.


title: How Can I predict the frequency of buying of a product

text: I'm trying to predict the frequency of buying a product based on their previous purchasing history. The approach I took as of now is,   1. Create a vector by finding the duration distance of each purchase (Xnt - X(n-1)t   2. Sort the resulted Duration distance   3. Remove Top 10% and bottom 10% results. Say X is 10% count then X = Number of elements / 10   4. Find the mean and standard deviation of remaining vector   5. Approximate duration: D = Mean + or - standard Deviation   6. Next purchase is: Last purchase date + Approximate Duration Is there any other way to model the prediction of buying? How can I enhance or remodel?


title: Cascade Combination of Kernel Functions

text: I have a question regarding machine learning and specifically kernel functions. Suppose we have a Kernel function, say $K(x)$, and also another distinct one, say $K'(x)$. I want to know is $K(K'(x))$ a kernel function as well? That is, if one feeds the output of a kernel function to another kernel, what does it mean? does it make sense or not? Another question is about the expected behavior of linear combination of well- known kernels,such as RBF, polynomial and MLP. Suppose the MLP kernel yields 60% of accuracy in a classification task and RBF yields 85%. Does necessarily the RBF+MLP yield a better accuracy compared to the one resulted by MLP?


title: Volume contribution decomposition in Log-log models

text: I am currently working on pricing analysis : the effect of competitor SKU pricing on the number of units sold of my SKU. The model was built on the log(units) sold. I want to measure the contribution of each of the independent variables in the model. I have a log-log model and want my contributions in actual units unlike a linear multiple regression model were my Total Contribution (or my predicted) = intercept + beta1*Price1 + beta2.Promo1 etc where beta1*Price1 measures my price contribution to predicted units sold, Here Im guessing in the log-log model there is a different way to go about the process. Any help would be appreciated. Thanks!


title: Factorial ANOVA

text: I want to compare means of two groups that have been controlled to have similar verbal characteristics: The proportion accuracy and mean reaction time (RT) to correct replies as two dependent variables. Priming, attention and Secondary Task Type are independent variables. Should I submit the data to two separate 2×2×2 mixed factorial design with Priming Condition (Repeated, Non-repeated) manipulated within subjects and Attention (Divided Attention, Full Attention) and Secondary Task Type (Perceptual, Conceptual) manipulated between subjects?


title: Binomial count data - use glmer, lmer or just average it all?

text: I'm new to both mixed-effects models and R, so please excuse me if this is a bit of a silly question! I'm struggling with choosing the best method to analyse my data for a speech perception study, and would appreciate some help! Some info about the data: In the study, participants listen to sentences in different accents and embedded in different noise levels, and then repeat what they understand. So for each sentence there is the number of words correctly repeated, and the total number of words in each sentence. I'm trying to figure out how to analyse this data in a mixed-effect model. I thought a binomial glmer would be best, as it's count data, but I had no idea if I was doing this right (particularly specifying the proportion of words correctly repeated), so I tried the same model (without interactions) in glmer and lmer. glmer:               > mixed.glmer<-glmer(formula=cbind(Correct, Max-Correct)~Accent+Noise_Factor+Accent+(1|Participant), data=noise.data, family="binomial")      lmer:               > all.lmer<-lmer(Correct/Max~Accent+Noise_Factor+(1|Participant), data=noise.data, REML=FALSE)      Then I compared the models:               >               Df     AIC     BIC  logLik deviance Chisq Chi Df Pr(>Chisq)         > all.glmer      6 13441.2 13480.1 -6714.6  13429.2                                 > all.lmer       7  2564.6  2610.0 -1275.3   2550.6 10879      1     <2e-16 ***      Which seems to suggest that the lmer model is a much better fit, but surely the data isn't being modelled well by linear methods? Then, a colleague suggested that I average across all noise levels for each accent (so I can treat the data using linear methods), giving only three data points per participant. So I tried this model, which seems to be a much better fit than either of the above models with all of the data:               > average.m1<-lmer(Average.score~Accent+(1|Participant), data=mixed.model, REML=FALSE)     > AIC       BIC    logLik  deviance      > -102.1759  -93.1426   56.0880 -112.1759       However, this seems to be throwing away a lot of data (particularly now that the noise factor has disappeared), even if the model does fit better - I thought the whole point of mixed models was to avoid doing the averaging step, like I'd need to for an ANOVA. This seems like very much like the easy way out of this problem! So my questions are:   1. Should I bother trying to fit a model to the whole dataset, or go with the easy way and average it?   2. Am I right to think that for the whole data set, the lmer model is a better fit than the glmer model, or is this like comparing apples and oranges?    3. If I do use the gmler model, am I specifying the DV correctly? It should be the proportion of correct responses (i.e.: Correct/Max). I'm finding it very, very hard to find instructions that I can actually understand!  Any help would be highly appreciated!


title: multiple regression and bootstrapping

text: I performed multiple regression and some of the assumptions are violated like normality of the outcome variable. Can I use bootstrap to correct for violated assumptions? Or doing it would be pointless!


title: Pricing Model for desktop computers

text: I am trying to implement a pricing model to automatic price desktops based on some attributes(say memory, disk and brand for now). I have collected some sold history data from the internet and they are more like in a "melted" format. Say, there are three transactions have been done in the history.               TransactionID    Brand    Memory    Disk    SoldPrice     1                Lenovo   1G        60G     $100     2                IBM      NA        60G     $150     3                Lenovo   4G        NA      $200     ...      And the data in I have looks like this:               SoldPrice    Attribute    Value     $100         Brand        Lenovo     $100         Memory       1G     $100         Disk         60G     $150         Brand        IBM     $150         Disk         60G     $200         Brand        Lenovo     $200         Memory       4G      I know I can use reshape(2) package in R to cast the data back to the wide format which might be a perfect machine learning problem using whatever methods(Linear regression, or omnipotent random forest). However, my gut feeling tells me will it actually be a perfect fit for Bayesian algorithm if I don't change anything? So you can easily get the probability distribution for each conditions. So it will be easier to calculate.               P(Price|Brand*Memory*Disk) = P(Price|Brand) * P(Price|Memory) * P(Price|Disk)      Then we can have a plot of the probability for different prices vs the given computer spec. will be very easy to calculate and actually will have a better accuracy than other models? I am more on the programming side so please provide more detail if you are mentioning graduate level math. I wish I can get some mathematical suggestion or guidance about my thought here instead of Stackoverflow.


title: What is the most suitable statistical techniques to be used to analyse realationship between Three(03) IVs and Two(02) DVs?

text: Research Topic: Effect of leadership leadership styles(Transformational, transactional and Laissez-faire) on Employee Competencies and Opportunities for their Development. Data Type : Ordinal.(0=never,1=rarely,2=sometimes,3=Most of the time,4=always) (For Every 5 variables) Independent Variables : Transformational Leadership, Transactional Leadership and Laissez-faire leadership Dependent Variables : Employee Competencies, Opportunities for Employee Development Research Method : Case Study Sample Size : 38 leaders(Unit of analysis) number of total participants to the research : 206 ( subordinates, supervisors and leaders themselves) Data analysis to be done using SPSS.


title: Random effects model with PLM: "System is computationally singular"-Error?

text: I am currently trying to estimate some panel data models in R using PLM package. This includes the estimation of basic pooled, fixed effects and random effects models. Therefore I make use of this code:               # read in data     mydata<- read.csv2("Panel.csv")     attach(mydata)          # define dependant variable     standarddeviation <- cbind(sd)     # define independant variable     x <- cbind(ratio1, ratio2, ratio3, ratio4, mean)          # Set data as panel data     pdata <- plm.data(mydata, index=c("id","t"))          # Pooled OLS estimator     pooling <- plm(standarddeviation ~ x, data=pdata, model= "pooling")     summary(pooling)          # Between estimator     between <- plm(standarddeviation ~ x, data=pdata, model= "between")     summary(between)          # First differences estimator     firstdiff <- plm(standarddeviation ~ x, data=pdata, model= "fd")     summary(firstdiff)          # Fixed effects or within estimator     fixed <- plm(standarddeviation ~ x data=pdata, model= "within")     summary(fixed)          # Random effects estimator     random <- plm(standarddeviation ~ x, data=pdata, model= "random")     summary(random)      Now here's the problem: I can without any problem estimate all models except for the random effects model. After entering the "random"-formula, R produces the following error:               "Error in solve.default(crossprod(X.m)) :        system is computationally singular: reciprocal condition number = 9.57127e-023"      First guesses:   * Linear combinations in x? A first guess would be that there are exact linear dependencies of the exogenous variables in x. The data is balance sheet data and I would like to explain the standard deviation (y) of a specific balance sheet position by other balance sheet positions (or the ratio of the position and the balance sheet sum). Of course, the variables in x are related to each other. For example some of the ratios are calculated by dividing by the mean which is also a separate independant variable. And the dependant variable, which is the standard deviation, is also calculated by using this mean. But again: There should be no EXACT correlation. But: If I exclude some of my exogenous variables, the problem disappears, but I have to include them actually.   * Unbalanced panel / NAs? The data is unbalanced and there are NAs. Fixed effects output says: n=16, T=18-40, N=455. Probably the unbalanced data or the NAs are the reason for the error? Traceback-Code: > random <\- plm(standarddeviaution ~ x, data=pdata, model= "random") Error in > solve.default(crossprod(X.m)) : system is computationally singular: > reciprocal condition number = 9.57127e-023 traceback() 8: > solve.default(crossprod(X.m)) 7: solve(crossprod(X.m)) 6: > diag(solve(crossprod(X.m)) %*% crossprod(X.sum)) 5: swar(object, data, > effect) 4: ercomp.formula(formula, data, effect, method = random.method) 3: > ercomp(formula, data, effect, method = random.method) 2: plm.fit(formula, > data, model, effect, random.method, inst.method) 1: plm(standarddeviation ~ > x, data = pdata, model = "random") Is there anybody who can give me a hint what this error does actually mean and especially: how to solve the problem? How do I have to correct the code in order to get results? Thanks a lot!


title: Smoother lines for ggplot2

text: This question probably has a simple solution, still the thing is I've written a code to plot mortality in 2 different groups and that is, death in obese patients vs not obese. Now their are 2 groups t1 and t2 (obese vs normal BMI). The graphs works just fine with this code, I'ts just that I wanted the lines to look smoother and not so jagged. I've tried stat_smooth but just cant get it to work. Im guessing that the code is longer then necessary.               allmortality<- read.table(header=TRUE, sep=";", text =                        "ar_year;t1_all_estimate;t2_all_estimate                       ar2001-2002;208.0960242;170.924898                       ar2003-2004;217.48718;151.6087781                       ar2005-2006;205.9649097;141.4196023                       ar2007-2008;212.2112923;135.2787361                       ar2009-2010;205.628018;136.4582058                       ar2011-2012;166.5654204;146.9776943")          require(ggplot2)            ggplot(allmortality, aes(x=ar_year)) +        geom_point(aes(y = t2_all_estimate, size=4)) +       geom_line(aes(y = t2_all_estimate, group=1)) +        geom_point(aes(y = t1_all_estimate, size=4)) +       geom_line(aes(y = t1_all_estimate, group=1)) +       theme_bw() +         theme( legend.position="none",              axis.text.x = element_text(angle = 45, hjust = 1),              axis.text=element_text(size=12),              plot.background = element_blank()              ,panel.grid.major = element_blank()              ,panel.grid.minor = element_blank()              ,panel.border = element_blank()) +       theme(axis.line = element_line(color = 'black')) +       ylab("Mortality") +       xlab("Year") +       ggtitle("mortality rates in obese vs not obese patients")       ![enter image description here](http://i.stack.imgur.com/i895e.jpg)


title: Bootstrap USL/LSL for non-normal data

text: This is new statistics space for me, so please excuse my ignorance. I have some data (N=180) that isn't normally distributed (verified in Minitab, P<0.005) and I want to use the data to get a +/- 4 sigma specification limit. Basically, I want to say: "We can expect 99.993666% of parts to have a measured dimension between DATA_LSL and DATA_USL." None of the distributions in Minitab fit with good confidence, and so I found that sometimes you can do bootstrapping, which is something like re-sampling your data to get an idea of the uncertainty in some statistic based on 1000+ random re-samples. I think I got that right, but I may be mistaken. Does anyone know how to bootstrap the standard deviation in Minitab? I found a macro that does it for the mean, and I can put in whatever confidence interval that I want. I think that only tells me about the mean, though, and not the standard deviation. There are some papers that talk about bootstrapping for LSL/USL, but nothing that I can find that's an example/macro in Minitab. Any help or advice would be greatly appreciated :) Thanks, Adam Hartman Mechanical Engineer - Flextronics stat_stackexchange.15.arrakis@spamgourmet.com


title: Causality, omitted variable bias

text: This might be a basic question, but I want to be sure that what I'm doing is right. I have a model that suggests that variable X causes both Y and Z. When I regress Y on X, or Z on X, I get positive and significant coefficients as expected. Now, when I regress Z on Y, I still get a positive significant coefficient. **Question 1:** is this an omitted variable bias? **Question 2:** is it legitimate to regress Z on Y **and** X to test whether the relationship between Z and Y is spurious? **Question 3:** if it is legitimate and if I get positive significant coefficients on **both** Y and X what does that mean? Does it mean "X causes Y and Z, but Y still has marginal explanatory power on Z"? Many thanks, Dave


title: liblinear one vs rest learn parameters

text: Liblinear (http://www.csie.ntu.edu.tw/~cjlin/liblinear/) does not support for probability estimates. Say I have three classes C1, C2 and C3. I want to learn the model paramters for each 'one vs rest' cases:                      C1 vs C2&C3,             c2 vs C1&C3 and            C3 vs C1&C2      How can I do it?


title: Statistical Model from distributions

text: I am aware of GLM. There the dependent variable $Y$ is assumed to be generated according to a certain distribution and for the mean it holds that $$ E[Y] = g^{-1}(X \beta). $$ Is there a common theory for models of the form $$ E[Y] = g^{-1}(X \beta) $$ and we assume distributions for the components of $X$ (e.g. $X_1$ is Guassian, $X_2$ is Gamma)? I am thinking about an MLE approach in this context. I guess this is related to GLM but maybe I am missing something. Thank you!


title: how to calculate future state probability if present probability is given by Markov chain

text: There is 50 elements transmitted from transmitter to receiver in 1st time where the probability of each element is not received by receiver is 0.1 and it will be received by receiver with probability 0.9. Again transmitter transmit same 50 elements 2nd time but this time i do not know what is the probability of each element is not received by receiver?Like this, transmitter in total 9 times transmit this 50 elements. I know that the second time the probability of each element is not received by receiver is 0.0494 and 3rd time it is 0.02284 but i do not how it is come. Anyone help is very much appreciated. Thanks in advance


title: Regression model F test on different response variables?

text: I am comparing two different models (multiple linear model). Basically the explanatory variables remain the same while response variable are different. $Model\ 1 - Y_1 = X_1+X_2+X_3$ $Model\ 2 - Y_2 = X_1+X_2+X_3$ I have results from R for both models. I am just wondering if two of these results eg $F$-statistics, $p$-value, $R^2$, RSS are comparable. So can I say model 1 is better because p-value is smaller. etc ? From what I learned from class, I only compared models which had the same response variable and different explanatory variables.


title: R: How to interpret the QQplot's outlier numbers?

text: How to interpret the labels with the outlier numbers when you plot the following in R (QQplot)               set.seed(1)     y <- rnorm(100)     x <- rnorm(100)     plot(lm(y ~ x), which=2)   # which = 2 gives the plot      It gives a number 61 on the top. What is it? I figured it might be the index of the outlier couple. It appears to be connected to a score of around `y = 3` and `x = 3`. But when:               cbind(y,x)[61,]          >  y         x      2.4016178 0.4251004       How to read these numbers in R's QQplot?


title: What test is this for endogenous variables?

text: Can somebody tell me whether the following R code (for econometrics endogenous variables) is for a Hausman test, a Nakamura test, or some other test?               etudes1 <- lm(EDUC  ~      EXPER+EXPERSQ+SMSA+SOUTH+FATHEDUC+MOTHEDUC                   )          Test    <- lm(LWAGE ~ EDUC+EXPER+EXPERSQ+SMSA+SOUTH                  +residuals(etudes1))     summary(test)


title: Comparison between groups with unequal sample sizes

text: I have calculated modulation indices (scalar) of neuronal responses for two categories of neurons $x, y$ recorded from two hemispheres $a, b$, with $n_{xa}=120, n_{xb}=80, n_{ya}=20, n_{yb}=20 $. I would like to know whether each category of neurons has significantly different modulation indices between the two hemispheres. Since I am not sure whether modulation indices are normally distributed I took the conservative approach and chose a non- parametric Wilcoxon rank sum test. I found that ${xa}$ is significantly different to ${xb}$ but ${ya}$ is not different to ${yb}$. However, given the small sample size of ${ya}$ and ${yb}$ I worry that the outcome might be a result of a Type II error. To account for the unequal sample sizes I came up with the idea to draw 20 samples from ${xa}$ and 20 from ${xb}$ (with replacement) for e.g. 1000 bootstraps, use for each sample the Wilcoxon statistic and come up with a distribution of p-values. Then calculate whether this distribution is significantly different to p=0.05. Do you find this approach valid? Any other ideas or suggestions would be appreciated.


title: Should I use Bayes' formula to generate simple data about business competition and the market place?

text: I run a small business, and I'm keen on understanding my market a little more. Market surveying and similar tools are not available to me, mainly due to time and budget. I don't require rock-solid data, just something that allows me to add a "pinch of salt" to my reasoning once in a while. Is it appropriate to use Bayes' theorem to determine belief that my customers will increase as the market (population) increases, even as my competition increases? This assumes homogenous geographic distribution of customer populations and competition, which is inaccurate, and assumes that all customers are drawn equally to all businesses, which is again (probably) inaccurate. But at least it gives me a starting point beyond my ill-informed hunches. I am not a statistician, which no doubt will be very evident, and have encountered several questions while trying to make sense of this approach. Apologies in advance for the number of questions here. In my limited understanding, Bayes' theorem might help me modulate my beliefs regarding my competition, the market place, and so on. I have a number of questions:   1. is this an appropriate use of Bayes' theorem? **Step 1: Constructing some formulae** Let's say: > $ P(C) = \text{probability customer numbers will increase} \\\$ $ P(M) = > \text{probability market population will increase} \\\$ $ P(K) = > \text{probability competition will increase} $ ...then I presume I should be able to construct something like the following: > $ P(C|M) = \frac{P(M|C)P(C)}{P(M)} $ and $ P(C|K) = \frac{P(K|C)P(C)}{P(K)} > $   2. Am I doing it right? :) **Step 2: Using appropriate data** The numbers I have are as follows:     * Professionals in my field, in the UK: approx 4,500.     * Growth in my field: approx 1,000 per year.     * Shrinkage in my field (due to retirement etc): approx 500 per year.     * Number of people in my city, and the rate at which it's growing.     * Percentage of a typical population that uses my profession's services: 10%.     * Month-on-month customer data in my business. I can probably get more precise figures regarding the number of professionals in my field. Regarding customer data, I have 18 months of past figures in a spreadsheet. In 9 of 17 completed months, customer numbers are higher than the previous month. This suggests to me that: > $ P(C) = \frac{9}{17} = 0.53 $ Regarding the market (population), let's say local government projects growth from 1,000,000 to 1,100,000 in the next year. This suggests to me that: > $ P(M) = \frac{1100000}{1000000} = 1.1 $ Regarding the competition, I'd assume: > $ P(K) = \frac{4500 + 1000 - 500}{4500} = 1.\dot{3} $   3. The figures may be accurate, but am I using them appropriately? **Step 3: The results** > $ P(C|M) = \frac{P(M|C) \cdot 0.53}{1.1} $ This is where I run out of steam. How do I determine $P(M|C)$, the probability of the market increasing given that my customers also increase? Since my business has negligible effect on market size should I simply consider that to be 100%, i.e. 1? My business will also have minimal effect on competition (or so I believe, although maybe new local businesses are deterred by its presence) so I'd assume that similarly $ P(K|C) = 1 $.   4. Is that a sensible approach? Somehow I suspect not. If that approach _is_ right, the answers become: > $ P(C|M) = \frac{1 \cdot 0.53}{1.1} = \frac{0.53}{1.1} = 0.4\dot{8}\dot{1} $ > > $ P(C|K) = \frac{1 \cdot 1.\dot{3}}{1.1} = \frac{1.\dot{3}}{1.1} = > 1.\dot{2}\dot{1} $   5. How should I then combine these two beliefs that customer numbers will go up as the market increases but competition also increases? Since they're probabilities, can I simply multiply them? That would suggest: > $ 0.4\dot{8}\dot{1} \cdot 1.\dot{2}\dot{1} = 0.584 $ Many thanks for any help you can offer, and apologies again for asking a lot of questions. I'm perfectly happy to be pointed toward useful reading materials, although my statistics and mathematics abilities are (clearly) limited, so novice-friendly info would be appreciated.


title: A naive question about the Kolmogorov Smirnov test

text: I have some times and I would like to test if they are distributed according to an exponential distribution. I can apply the Kolmogorov-Smirnov in python using `scipy.stats.kstest(data, 'expon')` . However, I assume I have to normalize my data somehow first otherwise it is comparing with an exponential distribution with unknown rate lambda. Is this right and how should I do that?


title: Relative variable importance values vs. magnitude of effect

text: I have ran a series of models to see which best fit the response variable and I got the following (for the model average of all models with a $\Delta AIC < 2$). I am currently learning models so this output is one of my first:                        Estimate Std. Error z value Pr(>|z|)         (Intercept)   102.7190     5.5300  18.575  < 2e-16 ***     HDr            -1.5495     0.3451   4.490 7.11e-06 ***     MF.vs.OF2      -7.6780     3.7507   2.047  0.04065 *       NHDp           -0.5145     0.2909   1.769  0.07695 .       NHDr           -1.4164     0.4663   3.037  0.00239 **      Site2           6.1477     2.7400   2.244  0.02485 *       tide.h.l2      -7.2546     2.6914   2.695  0.00703 **      tide.inc.out2  -5.8486     2.6187   2.233  0.02553 *       HDp            -0.3773     0.2732   1.381  0.16731         mean.for.rate  -0.3966     0.3220   1.232  0.21807         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Full model-averaged coefficients (with shrinkage):       (Intercept)        HDr  MF.vs.OF2       NHDp       NHDr      Site2  tide.h.l2 tide.inc.out2        HDp       102.718962  -1.549499  -5.734171  -0.239550  -1.416373   5.336532  -7.254627     -5.848553  -0.044795      mean.for.rate          -0.081734          Relative variable importance:       (Intercept)           Age           HDp           HDr mean.for.rate      MF.vs.OF          NHDp          NHDr               1.00          0.00          0.12          1.00          0.21          0.75          0.47          1.00               Site      tide.h.l  tide.inc.out               0.87          1.00          1.00       I have seen the estimates described as the effect of the variable and this is discussed in results sections as an important value to report (in regards to the size of them and their direction (+ve/-ve)). (the paper I was reading was stating that those with the bigger or smaller numbers had the greatest effect, even quoting that one was 48% lower than the other). However if this is what is reported and discussed, why would the relative variable importance vary in relation to the estimate? It seems that this should also be looked at but am not sure how the z and p values are calculated from a model. Therefore, I would like to know which is more important when trying to discuss the findings. I admit that my knowledge is limited, but I would like to grasp this in simple terms if I could. Also, I used a dredge command to find all the model combinations from a global model and found seven models (out of 1024) had a $\Delta AIC < 2$. I thought that this was a good cut-off for averaging, but could technically use any delta value as the cutoff. Thank you.


title: How to tell when factors "disagree" in linear regression to produce noisy predictions?

text: I use a regression as my predictor. Let's say my regression is $y = a_1 x_1 + a_2 x_2 + a_3 x_3$ I realized that in practice, when my prediction is way off, it's usually because one factor significantly skewed the prediction. For example, $x_1, x_2$ are both slightly negative, while $x_3$ is very positive. In this case (when factors "disagree" with each other), I would rather my predictor not do anything, than reporting a "controversal" prediction. In my particular application, I don't have to generate a prediction every single time (it can simply returns "no clue"), but when I do, I would rather it be correct. What is the best way of doing this? Thanks!


title: Proportion Estimates - Shrink to the Mean Based on Sample Size

text: **Set-up:** One of the fundamental elements of direct marketing is choosing which lists of prospects to select and send an offer to. There are hundreds of lists on the market to choose from and each list will contain a specific number of records you can send an offer to. I have a table of historical performance of list. For each, the number of times a promotion was sent to someone on that list, along with the number of responders is given. So, the table might look like this - it is a made up example with only a couple of lists versus the hundred actually but the very small response rates (p) are realistic.               lists<-data.frame(listName=c('A','B','C','D','E'),numberPromoted=c(1000,5000,25000,35000,1654),numberRespond=c(3,45,98,350,5))     lists$p<-lists$numberRespond/lists$numberPromoted      I am looking to use linear programming to optimally allocate list selection in future campaigns. So, given the current number of records available in each list at the preesent time and the best estimate of the response rate of the list (p), choose how many records from a list to select, given constraints. **Problem:** My question is how to best go about estimating p. Given these small response rates, I would like to shrink those without "adequate" information to a mean value. The naive way is to say that there has to be at least X number of records promoted before we trust the p estimated empirically. I am wondering if there is any other method that makes sense? For example, does it make sense to treat the list as a random effect and shrink the estimate? I **THINK** this code is shrinking the estimates of the response rate (p) to the mean rate for smaller cells.               install.packages("lme4")     library(lme4)     install.packages("arm")     library(arm)          #make individual 1/0 records     Lists<-data.frame(     listName<-c(replicate(1000,"List A"),replicate(5000,"List B"),replicate(25000,"List C"),replicate(35000,"List D"),replicate(1654,"List E")),     outcome<-c(replicate(997,0),replicate(3,1),replicate(4955,0),replicate(45,1),replicate(24902,0),replicate(98,1),replicate(34650,0),replicate(350,1),         replicate(1649,0),replicate(5,1))     )                   #random intercept model         mod<-lmer(outcome~1 + (1|listName),data=Lists,family=binomial)              invlogit(coef(mod)$listName[1,1])         invlogit(coef(mod)$listName[2,1])         invlogit(coef(mod)$listName[3,1])         invlogit(coef(mod)$listName[4,1])         invlogit(coef(mod)$listName[5,1])          List    Modeled P   "TRUE" P     A   0.004241    0.003000     B   0.008602    0.009000     C   0.003981    0.003920     D   0.009924    0.010000     E   0.003961    0.003023                                #compare raw versus shrunk modeled estimates          raw<-sqldf("select listName, avg(outcome) as raw_p from Lists group by listName")         meanRaw<-sqldf("select avg(outcome) as total_p from Lists")         modeled<-c(         invlogit(coef(mod)$listName[1,1]),         invlogit(coef(mod)$listName[2,1]),         invlogit(coef(mod)$listName[3,1]),         invlogit(coef(mod)$listName[4,1]),         invlogit(coef(mod)$listName[5,1]))              plot(x=as.numeric(raw$listName),y=raw$raw_p, col="blue")         points(x=raw$listName,y=modeled, col="red")         abline(h=meanRaw)


title: What can be gained by including the distribution of the dependent variable in a GLM?

text: In this discussion on probit vs. OLS models 'Conjugate Prior' mentions that probit and all GLMs allow to specify a DV distribution. How is the GLM improved by specifying a DV distribution and what is the implied DV distribution for OLS then?


title: The best way to analyze repeated measures from limited subjects?

text: I realized there are options to analyze repeated measures like anova, gee, mixed model. Here I have a challenge situation like only 4 measures from 10 subjects. I'm wondering, if I have to, what the best way to analyze it? Such that:               df<-data.frame(y=rnorm(40, 10, 2),     time=rep(1:4, 10),     age=rep(as.integer(rnorm(10, 20, 3)), each=4))      The goal is to see whether y is changing with time and age.


title: HMM-forward backward algorithm

text: Let $x_t, \, t=1, \dots ,T$ be a time series and suppose that $x_t | \xi_t \sim N(\mu_{\xi_t},\sigma^2_{\xi_t})$, where $\xi_t \in [1, \dots,K]$ is a group indicator (or regime or state), the parameter of the normal distribution depends on the value of $\xi_t$ and we suppose that $P(\xi_t=k|\xi_{t-1}=j) = \pi_{jk}$. I want to estimates this in a Bayesian framework and then the prior are $\pi_{j1},\pi_{j2},\dots, \pi_{jK} \sim Dirichlet(\alpha_1,\dots, \alpha_k)$, $\mu_{k} \sim N(M,V)$ and $\sigma_{k}^2 \sim InvGamma(a,b)$. I suppose $\xi_0=1$. I am trying to understand how sample from the full conditional of $\xi_t, t=1,\dots,T$. Can someone explain me if I have to use the forward, the backward or the Viterbi algorithm or something else? and how to use it? I am really confuse about all these stuff :-). A link where it is well explained is accepted as well.


title: How to explain the large p.value difference between these two models

text: I fit two gee models with two different correlation structures, exchangeable vs ar(1), which resulted in very different p.values. I'm wondering what reasons have led that. Would somebody offer an explanation?               > summary(fit1)          Call:     geeglm(formula = LAU ~ SAMPLENO, data = bd, id = MAGE, corstr = "exchangeable")           Coefficients:                  Estimate   Std.err   Wald Pr(>|W|)         (Intercept)  0.067349  0.004537 220.32   <2e-16 ***     SAMPLENO    -0.002800  0.000947   8.74   0.0031 **      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Estimated Scale Parameters:                 Estimate  Std.err     (Intercept) 0.000425 0.000166          Correlation: Structure = exchangeable  Link = identity           Estimated Correlation Parameters:           Estimate Std.err     alpha    0.206   0.249     Number of clusters:   9   Maximum cluster size: 8                          > summary(fit2)          Call:     geeglm(formula = LAU ~ SAMPLENO, data = bd, id = MAGE, corstr = "ar1")           Coefficients:                 Estimate  Std.err   Wald Pr(>|W|)         (Intercept)  0.06392  0.00588 118.25   <2e-16 ***     SAMPLENO    -0.00142  0.00219   0.42     0.52         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Estimated Scale Parameters:                 Estimate  Std.err     (Intercept) 0.000427 0.000172          Correlation: Structure = ar1  Link = identity           Estimated Correlation Parameters:           Estimate Std.err     alpha    0.659   0.221     Number of clusters:   9   Maximum cluster size: 8


title: PyMC: how can I use a custom sampler for one specific variable in a model?

text: Basically I have a function to efficiently sample one variable in my model that I would like PyMC to use, along with its normal samplers for everything else. Since I have this function I don't have a need to define a log likelihood for the variable either. Can I do this using the `Stochastic` class and making a `random` function? Or do I need to make a custom step method?


title: Is the policy function $\pi$ in Reinforcement learning a random variable?

text: I was reading Andrew Ng's lecture notes on reinforcement learning and on page 3 he defines the value function: $$V^{\pi}(s) = E[R(s_0) + \gamma R(s_1)+\gamma^2 R(s_2) + ... |s_0 = s, \pi]$$ Which means the expected total pay off, given that we start on state s and execute policy $\pi$. However, in his footnotes he says that this notation is a little "sloppy" because he says that $\pi$ isn't technically a random variable (even though the notation implies that both $\pi$ and s are random variables, though it makes sense that s is a r.v. since we don't always know for sure what state we will end up after doing some action a). My question is, **why isn't $\pi$ a random variable?** If its not a random variable is it in a reinforcement learning that we are just looking for some policy that is the "best" and was chosen by "nature" somehow? Does it mean that we are not allowed to have some prior belief on which $\pi$ might be true? Is reinforcement learning or at least, the value function here restrained to a frequentist point of view? Would a better notation for that equation be: $$V^{\pi}(s) = E[R(s_0) + \gamma R(s_1)+\gamma^2 R(s_2) + ... |s_0 = s; \pi]$$ These were my thoughts so far: $\pi$ is the policy function, its a function that maps states deterministically to actions $\pi(s) = a$. However, I didn't really see why reinforcement learning had to be restricted to a frequentist interpretation. I seemed reasonable to me that $\pi$ could be a r.v. and we instead try to execute the expected policy over all policies or something along those lines (I am not trying to make this idea too precise, but hopefully the idea/concept makes sense). Is it just that Andrew Ng is introducing the concepts of reinforcement learning first a frequentist, as it might be the easiest to understand?


title: Given the information, what is the probability that the mean number of drinks sold per week is between 700 and 720

text: The question . ![a busy cat](http://i.imgur.com/4qodjPL.png) I'm kinda a bit confused on how to make use of the information to answer the questions. I think I need to use the Bayesian Method to answer this question but I'm not sure if I'm right. would the prior density of the mean be $\mathcal{N}(738,13.4^2)$? Then what do I do with the $\mathcal{N}(692,42.5^2)$ Given in the second half of the question? I'm really clueless in using the Bayesian method and would like some guidance in understanding this question.


title: Cross validation MATLAB code for multi-output regression

text: I am using multi-output support vector regression (MSVR) which predicts multiple outputs at a time. I want to know how can I initially select parameters for my model, and then how to do cross validation. The input `x` is 100×78 (100 samples, 78 features) and the ouput `y` is 78-dimensional.   Please help me with MATLAB code.


title: What does Language Model look like?

text: I am working on a machine learning + NLP project. The corpora is from **a very specific domain**. Someone tells me I need a language model for that specific domain. So I decide to build one myself since there's none yet. My questions are:   1. **What does a language model look like?** Is it in the form of a bunch of n-gram words and probabilities? Maybe like this? term 1, probability 1 term 2, probability 2 term 3, probability 3 term 4, probability 4 ...   2. **How to obtain a language model?** Should I get a corpus of the domain and just summarize the term frequencies?   3. **How to use the language model?** Suppose I get a new document for the domain, and I want to vectorize the document, should I directly assign probabilities from the language model to the terms in the document? and generate. If so, the same terms from different documents will have identical probabilities.


title: Probability estimation and trapezoidal rule

text: proc iml;          call randseed(4545); * initialize the stream (like streaminit);          x = J(5000,1,.); * pre-allocate space for random numbers;          call randgen(x,'normal',0,1); * fill x with N(0,1) deviates;          y = y + (x**2 - 3*x**3 + 5x < 1);          p = y / 5000;  * MEAN acts on each matrix column;          se = sqrt(y*(1-y)/5000); * VAR, but not STD or STDERR, also acts on columns;          print "IML STEP: estimated mean of sqrt(abs(X)) is" p "with standard error" se;  * use PRINT, not PUT;      I'm trying to use monte carlo integration with proc iml to estimate the probability that x**2 - 3*x**3 + 5x is less than 1. What am I doing wrong? Do loops are not allowed by the way.               proc iml;          start tr(x,y); * create function called tr;     /* handle bad case */     /*if x<0 then do; return (.); end; * send back missing;*/               N = nrow(x);     dx = x[2:N] - x[1:N-1];     ymean = (y[2:N] + y[1:N-1]) / 2;     return(dx` * ymean );     finish tr;               x = do(-2,5,0.01);     print "Integral of y over x is" ( tr(x,sin(x##2)) );      I keep receiving the (execution) invalid subscript or subscript out of range. How do I solve this problem?


title: Train timeseries model in R and predict purely on out-of-sample data

text: I have two data files "train.txt' and 'test.txt' with single columns of data. I want to learn a model only on training-data and generate an output on test- data. I can't seem to find a way to do that. Most `predict` invocations seem to be just starting from the end of training data. I have always been c++ heavy, and am just learning R. I tried the `forecast` package in R.               train_data_ <- read.table ('train.txt');     train_data_ <- as.matrix ( train_data_ );          test_data_ <- read.table ('test.txt');     test_data_ <- as.matrix ( test_data_ );          fit_train_ <- ets(ts(train_data_));     fit_test_ <- ets(ts(test_data_),model=fit_train_);     onestep <-fitted ( fit_test_ );     fit_test10_ <- ets(ts (test_data_[1:10], model=fit_train_);     onestep10 <-fitted ( fit_test10_ );     head (onestep10); # print 10 lines     head (onestep); # prints 10 lines     # These are different.      This fails my goal of predicting the next step in test data in a way that the future testing data does not affect the prediction. Just to explain what I mean by not looking into the future for prediction this link might help : Problem Description


title: Statistical comparisons between large data sets

text: Currently, I am looking for the correct (or suitable) statistical method to compare 4 very large datasets (n = 31 million each), that are based on an experiment where a continuous variable was measured in response to 4 different discrete temperatures. Each data set has pretty much the same median, and the means and standard deviations do not vary much. However, each data set is not normally distributed, and are only vaguely lognormally distributed. I am trying to determine a way to test the null hypothesis that increased temperatures do not affect the mean of each data set. Am I right in my thinking that ANOVA would not be suitable, (nor any parametric test)? I have looked at the non parametric Kruskal-Wallis statistic, but are not sure if this is a suitable. So my questions are, am I correct in my assumptions above? and what would be a suitable statistical method to prove the null hypothesis as described above?


title: Probability usage in a customer event (modified question)

text: This question is slightly related to Probability usage in a customer event In a customer event we are planning users will choose tickets that contains one letter from alphabet {A,E,L,P}. The ticket is not given back and user can get (buy) as many tickets as he wants (but there will be some limitation - e.g 200 tickets per user). There is no maximum number `n` of involved users (but for simplification we can count that there will `n` = 1000 users). In the event I need to guarantee that there will be min. `x` winners (e.g. `x` = 20 people) who have collected whole word "APPLE". Question: I am looking for an equation how to count how many tickets of each letters do we need to print. And I would like to minimize the number of printed tickets. Is there a way how to count the number of each letter with respect to variables `x` and `n` and with probabilities of each letter - `P(A)`, `P(E)`, `P(L)`, `P(P)`? I'll want to transfer this equation as a model into Excel. And more modification - there can be a "joker" ticket. If you get a Joker, you can use this ticket as replacement of letter you want. Probability of getting a Joker ticket should be very low (e.g. 1-5%).


title: Shrinkage parameter in Adaboost?

text: I'm unclear how the shrinkage parameter works in Adaboost. I understand the concept of shrinkage in the theoretical sense related to ordinary least squares, but I'm not sure how to interpret this parameter in relation to Adaboost.


title: Mean interval estimation when data are not binomial

text: When a set of data are binomial I have a lot of methods available to me (Wilson score interval, Wilson score interval with continuity correction, Jeffreys interval, Clopper-Pearson interval, Agresti-Coull Interval) other than the standard Wald interval estimation method (see Interval Estimation for a Binomial Proportion Author(s): Lawrence D. Brown, T. Tony Cai and Anirban DasGupta). However, I can't use these methods for data that is itself proportional in the sense of it sitting in [0,1]. For example consider a sample of n respondents, giving a score for a question X as:               r_1 = 0.4     r_2 = 0.3     ...     r_n = 0.2      How would I calculate the confidence interval for a mean based on this data? Is there any way other than bootstrapping?


title: How to interpret the dendrogram of a hierarchical cluster analysis

text: I think I have some idea about hierarchical clustering but there a few subtle questions. I use the R example below: ` plot( hclust(dist(USArrests), "ave") ) `   1. What exactly does the y-axis "Height" mean?    2. Looking at North Carolina and California (rather on the left). Is California "closer" to North Carolina than Arizona? Can I make this interpretation?    3. Hawaii (right) joins the cluster rather late. I can see this as it is "higher" than other states. In general how can I interpret the fact that labels are "higher" or "lower" in the dendrogram correctly?  Thanks for any comments or references.![enter image description here](http://i.stack.imgur.com/tzoq6.png)


title: How to deal with heteroscedasticity in nonlinear system of equations?

text: What is the best way of dealing with unknown form of heteroscedasticity in nonlinear system of equations? Thanks in advance.


title: Statistical forensics: Benford and beyond

text: What broad methods are there to detect fraud, anomalies, fudging, etc. in scientific works produced by a third party? (I was motivated to ask this by the recent Marc Hauser affair.) Usually for election and accounting fraud, some variant of Benford's Law is cited. I am not sure how this could be applied to _e.g._ the Marc Hauser case, because Benford's Law requires numbers to be approximately log uniform. As a concrete example, suppose a paper cited the p-values for a large number of statistical tests. Could one transform these to log uniformity, then apply Benford's Law? It seems like there would be all kinds of problems with this approach ( _e.g._ some of the null hypotheses might legitimately be false, the statistical code might give p-values which are only approximately correct, the tests might only give p-values which are uniform under the null asymptotically, etc.)


title: Visualization of a multivariate function

text: This is somewhat vague, but suppose you have a black box function $f(x_1,x_2,\ldots,x_k)$, for which you have code, and you are interested in the behaviour of $f$ when the $x_i$ are i.i.d. standard Gaussian random variables. What are some good ways to visualize this function? To make it easier, we may assume that $k$ is smallish, say less than 10. One particular relationship of interest is how $f$ varies with one of the input, say $x_i$. An easy way to visualize this relationship would be to sample the function for fixed values of $x_i$ while varying the other input (either in a structured way, or randomly, say), then box-plotting, which could show how the mean trend is affected by $x_i$, but also whether the scatter is affected (i.e. heteroskedasticity). However, interaction between $x_i$ and the levels of the other input might be masked by this approach. What I am looking for is somewhat open-ended. I do not have a particular hypothesis that I am testing, but rather am looking for new ways of visualizing the response which might reveal peculiarities of the function.


title: Books on cluster algorithms

text: I'm searching for books on the basic k-means and divisive clustering algorithms. I'm interested in the pros and cons of both. It's a part of my bachelors thesis, I have implemented both and need books to create my used literature list for the theoretical part. Also, is there a book on "the curse of dimensionality"? Thanks!


title: Determination of effect size for a repeated measures ANOVA power analysis

text: This is a follow-up to the repeated measures sample size question. I am planning a repeated measures experiment. We record energy usage for 12 months, then give (a randomly assigned) half of the customers continuous information about their energy usage (perform the treatment), and record their energy usage for another 12 months. A similar study performed in the past showed a 5% reduction in energy usage. I want to estimate the required sample size using $\alpha=0.05, \beta=0.1$. G*Power 3 has a tool for repeated measures power analysis. However, it requires two inputs that I am not entirely familiar with:   * $\lambda$ - the noncentrality parameter (How do I estimate this?)   * $f$ - the effect size (I believe that this is the square root of Cohen's $f^2$) According to Wikipedia's effect size page: > Cohen's $f^2= {R^2_{AB} - R^2_A \over 1 - R^2_{AB}}$ where $R^2_A$ is the > variance accounted for by a set of one or more independent variables $A$, > and $R^2_{AB}$ is the combined variance accounted for by $A$ and another set > of one or more independent variables $B$. However, my expected 5% change in energy consumption does not tell me how much variability will be explained. Is there any way to make this conversion? If you know of a way to do this power analysis in R, I would love to hear it. I am planning to simulate some data and try using lmer from the lme4 package.


title: Papers on Bayesian factor analysis?

text: I am interested in fitting a factor analysis-like model on asset returns or other similar latent variable models. What are good papers to read on this topic? I am particularly interested in how to handle the fact that a factor analysis model is identical under a sign change for the "factor loadings".


title: Lumping in Markov process with absorbing states

text: I have a four-state, discrete time Markov process with time-dependent transition matrices such that after a given time T the matrices become constant. The idea is people in a program leaving the program in a variety of ways. Everyone starts in state 1, and states 2, 3 and 4 are absorbing, but state 4 represents the fairly small percentage of people who are 'lost in the system' - in other words state 4 represents our ignorance of what happens to people rather than a genuine outcome. I would like to use lumping to put those in state 4 in with those in state 1 and run this as a three-state system, and compare this with the naive approach of running this as a four-state system then apportioning those who are asymptotically in state 4 into states 2 and 3 according to their relative proportions. (in other words, p_2/(p_2 + p_3) of those in state 4 go into state 2 after the system is run to infinite time and similar for state 3) From some rough scribblings it doesn't seem that these two methods give the same results, so it would be good to get an idea on the error involved. To this end, here's my question: Can I have pointers to the literature on lumping in Markov chains (or related) that would apply - even roughly - in this example? Or otherwise some words of advice on how to approach this.


title: Lewandowski algorithm demand forecasting

text: I came across the Lewandowski method of demand forecasting in JDA Demand. Please help me understand at a high level the methodology it uses. I found a paper by Robert Hyndman titled "A state space framework for automatic forecasting using exponential smoothing methods" and it uses this method as one of methods they compare their algorithm to in the paper. Currently for us this is a black box, we want to get some high level understanding so that we can better fine tune the parameters they have provided as part of the software. It would be great if you can share some thoughts about the Lewandowski algorithm and point to some references that I could use for further research.


title: Interactions terms and higher order polynomials

text: If I were interested in fitting two-way interactions between a linear explanatory variable $a$ and another explanatory variable $b$ that has a quadratic relationship with the dependent variable $y$, would I have to include both the interaction with the quadratic component and the interaction with the linear component in the model? E.g.: $$ y\sim a+b+b^2+ab+ab^2 $$ In turn building on my previous thread: Curvature terms and model selection, if this were a model selection analysis using `MuMIn` in R, with many explanatory variables, would the outputted models containing an interaction term that included a quadratic term $a:b^2$ only be valid if the interaction term with the linear component $a:b$ was also present in that same model, as well as $a$, $b$ and $b^2$ as direct effects?


title: Analysing transformed data

text: Is it possible to carry out normal multiple linear regression when the dependent variable and one predictor variable have been transformed using square root transformation? (as they did not follow normal distribution). Is there any back transformation necessary for the R2 value, coefficients and confidence intervals? Thanks!


title: ARMA MLE - matrix setup

text: I am trying to estimate ARMA(2,1) series parameters via ML, but unfortunatelly I am confused of setting up Matrix "Gamma" . I have seen plenty of articles regarding this, but it seems like this an obvious thing, that is not mentioned or described anywhere, i.e.: "$L(\Gamma, X)= (2\pi)^{-\frac{n}{2}} \det(Γ)^{-\frac{1}{2}} \exp\left({-\frac{1}{2} X´ Γ^{-1} X}\right)$ where X n is the data, and $\Gamma$ depends in a complicated way on the parameters" Can anybody give me an advice please? Thank You.


title: Residual Value Prediction For Used Electronic Products

text: I am trying to predict the long term residual value of a product with only the releasing price. I have collected some data off the Internet related with one phone type, and it is pretty obvious that the residual value (say eBay and Amazon used pricing) of a used electronic product (say iphone) is decreasing exponential. Like the plot below: ![enter image description here](http://i.stack.imgur.com/l2Q9m.jpg) > >     par(mfrow=c(3,1)) >     x <- 1:100 >     y <- exp(-0.03 * x) * (2 + rnorm(100)) >     plot(x, y, main="raw data", ylim=c(0, 4)) >     plot(diff(y), main="diff", pch=4, ylim=c(-2, 2)) >     plot(y ~ as.factor(as.integer(x/10))) >   Clearly, the derivative is not constant and also decreasing. So I think exponential model might be a good fit. However, one interesting thing here is the longer the time goes since the release.. the less dispersive the data is. In another way, when phone first released. People started selling used phones, and since the value is the highest, say $700. They can possibly sell at 600 USD or even 300 USD if they want. But 3 years later, the phone value is much lower, say even a new phone is only $200, then when they sell old phones, the prices are probably at the range of 100~200, a much narrowed band. In my scenario, I am more interested in the residual value in **_a longer term, which is at the narrowed band._** Here is my question:   1. should I only look at the attenuation coefficient, the b from `Exp(b*x)` for different types of phones... and predict the new b for the phone that I am going to predict based on all the data that I have.  > >     PhoneModel   b >     iPhone3      -0.04 >     iPhone4      -0.05 >     iPhone5      -0.03 >     iPhone6      TBD >     2. or should I only look at the time period that I am interested in. say 2 ~ 3 years (or 36 months time point) and do the modeling in that range since the data is much narrowed.. To avoid the dispersion brought by the 0 ~ 2 years data.. > >     PhoneModel residualRatio >     iPhone3    20% >     iPhone4    25% >     iPhone5    30% >     iPhone6    TBD >   Again, I am trying to predict the long term residual value of a product with only the releasing price.


title: How to estimate Cohen's D using the independent groups definition but on a paired samples data?

text: I want to calculate cohen's d with confidence intervals for a paired samples designs. Some authors suggest that you use the paired t test value to adjust for the correlation between measures (Rosenthal, 1991). However, Dunlop et al 1996 suggest that the correlation between the paired samples should not be included. In particular, if such a correlation is used, then results are not readily comparable to between subjects effects. Instead, they recommend using independent samples formula: $$d= \frac{\mu_1-\mu_2}{\sigma_\textrm{pooled}}$$ ### Question **How do you estimate confidence intervals for Cohen's d in a paired samples design where Cohen's d uses the formula above?** The way I have estimated it so far is to use the R function `ci.smd` from the MBESS package. This function is designed for independent samples tests and takes the sizes of two independent groups as input. In the code below I use `i` to represent the number of subjects (equal, as the same participants are in each group) and use an effect size of 0.8 as an example.               library(MBESS)     x.ci <- ci.smd(smd=0.8, n.1=i,n.2=i)      **Is it appropriate to use such function`ci.smd` for confidence intervals for paired samples?**


title: how to calculate cost matrix for cost sensitive meta classifier?

text: i have an imbalanced data set with 12 classes and want to use cost sensitive ensemble for classification but i don't know how to calculate cost matrix for the classes. my question: what costs should i insert in cost matrix and how to calculate it for cost sensitive meta classifier in WEKA?


title: Question Regarding Nested-Cross Validation

text: I am slightly confused of the setup for tuning hyperparameters under gridsearch. Below is Python-style pseudo code for my understanding.               errors = []     for fold_training, fold_testing in 10-fold:              # Use 10-Fold Grid Search on the current fold          grid = GridSearchCV(fold_training, fold_testing, cv=10)         error.append( grid.error_of_best() ) # Add error to array          print "Cross validation score",error.average()      Is this proper practice for nested cross validation? This would train 100 models. Lastly, how do we get the optimal parameters from a nested cross-validation? Should we just use the highest scoring one out of all 100 models? I am quite confused on this.


title: Validating a manual classification

text: I am classifying some political articles by their bias w.r.t. a particular issue into 5 categories ( _strongly in favour, mildly in favour, neutral, mildly against, strongly against_ ). For most cases, the classification is quite simple, but I want some form of validation to be able to make sure that other subjects would classify the articles almost in the same way. As I have about 2,000 articles and no funding, I can't really ask other people to classify the entire dataset. An idea would consist of extracting a random sample and make other people classify the sample, and then look at the inter-subject agreement. If this approach makes sense, how do I determine a suitable sample size for this kind of validation w.r.t. to an error (e.g. 1%)? Cheers, Mulone


title: Bank acquisition effect

text: Given that I have 4 banks (bank A, B, C and D) and the data (stock value) I have is from 1998 to present year. At the year 2002, Bank A acquired Bank C. At 2005, Bank A acquired Bank D. Until there is only 1 surviving Bank (A) at present. I want to test whether each merge was effective (in terms of Bank A) or not. I consulted 2 professors and one recommended checking for structural change (and see change in level and volatility) while the other recommended intervention analysis. I'm torn to which one to use. Which do you think is best suited? Thanks in advance. :)


title: Gaussian Naive Bayes really equivalent to GMM with diagonal covariance matrices?

text: Murphy writes that a multivariate Gaussian used in a generative classifier ("Gaussian discriminant analysis"), i.e., $p(\mathbf x|y=c,\mathbf\theta) = \mathcal{N}(\mathbf x|\mathbf y_c,\mathbf \Sigma_c)$, is equivalent to Gaussian Naive Bayes if $\mathbf \Sigma_c$ is diagonal. This makes sense to me since GNB assumes independent features, each following a normal distribution. However, I've also read in multiple papers, for example this one, about an alleged equivalence of Naive Bayes and a mixture of Gaussians with diagonal covariance matrices. How can this be the case when a mixture can approximate any distribution given enough Gaussian base components? It seems to me that the GMM is more powerful than both Naive Bayes and Gaussian discriminant analysis in that it can provide better likelihoods.


title: How to recover linear probability model for binary classification?

text: In some binary classification problem I assume that the probability for a positive is exactly linearly dependent on the features $P(y=1|x_1,x_2,x_3,\ldots)=\beta_1x_1+\beta_2x_2+\beta_3x_3+\cdots$. Now for the training I only observe the realization 0 or 1 for $y$. Therefore the probability isn't directly observable, but would rather be some local density (which is hard to estimate for high dimensions). Is there some type of regression that can determine the coefficients $\beta_i$? I could try normal regression, but it doesn't seem obvious that any loss function will be able to converge to the same $\beta_i$. What type of regularization and loss function would be most appropriate?


title: Probability Calculation

text: Assume that a typical computer manufactured by a company lasts 10 Months and that the standard deviation is 50 days. Computer life follows a normal distribution. What is the probability that a computer made by this company will last at most 1 Year? Assumption is that one month has 30 days. Can you please explain how this is calculated ? Thanks


title: Combining treatment as usual and experimental placebo control groups in one meta analysis

text: I'm conducting a sr and meta analysis of interventions that target children with ADHD in school settings. The studies reviewed include both treatment as usual comparisons and experiments placebo controls. Can these types of studies be combined in a single meta analytic synthesis? I wondered about conducting a moderator analysis using design type as a predictor variable? Any comments welcome! thanks,


title: How do I create subscales based on sets of items in SPSS

text: I have a 37-question scale, where each item belongs to one of four scales. The items for each scale are interspersed. How can I compute a scale score for each of the four scales? I need to be able to analyse each subscale independently and then compare scale responses on and independent variable.


title: Assign attributes / categories to users based on their activity / likes

text: I have a very practical classification problem for which I need some help. I have a database of users along with their activity / likes for a number of car models. I also have the category each of these car models belong to (e.g. sport, coupe, family, etc). Ultimately, I want to either assign or not one or more of these categories to the users based on their overall activity. What I have done so far is 1) count the "likes" of each user on each category, 2) find the average and standard deviation of "likes" for each category considering all users together, and 3) perform some sort of classification using the following simple rule: if the number of "likes" of user X for category Y is bigger than mean + 2 standard deviations of the corresponding "likes" of the whole population then assign category Y to user X otherwise do nothing. My questions are: a) does the above process make any sense and 2) should I somehow take into account the total number of "likes" of each user as well as the total number of car models of each category - in other words perform some sort of data normalization - or this is somehow already incorporated in the already mentioned process?


title: Random effects are not normal

text: ![enter image description here](http://i.stack.imgur.com/LL3GL.jpg) I'm doing a neg binomial glmm in glmmADMB. The random effects are not quite normally distributed (see attached images). How concerned about this do I need to be? If this is too much of a violation, are there other options for analyzing longitudinal data?


title: Problem with mixture discriminant analysis in R returning NA for predictions

text: I'm trying to do a mixture discriminant analysis for a mid-sized data.frame, and bumped into a problem: all my predictions are NA. After tracing through way too much code, I figured it had something to do with the fact that some of the coefficients in the mda turn out to be NA. I've created a smaller data.frame that still has the problem:               dfr<-structure(list(min_GCs_last_3_bases = structure(c(1L, 2L, 1L,      1L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L,      1L), .Label = c("zero", "one"), class = "factor"), cq = c(-0.334138707578632,      -0.586150643906373, -0.474578712720667, -0.220268433139143, -0.486876877353103,      -0.0912154554410563, 0.00341593805213764, 0.713424582672338,      -0.448914652233824, 2.94156773625266, -0.0954835698859817, -0.238125375419562,      -0.448914652233824, -0.0290299100261503, -0.479666688891261,      -0.596316272919155, -0.0919665779183363, -0.394427887135795,      -0.396749521580081, -0.338156015653477), Coverage = structure(c(2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L,      1L, 2L, 2L), .Label = c("failure", "success"), class = "factor")), .Names = c("min_GCs_last_3_bases",      "cq", "Coverage"), class = "data.frame", row.names = c("SMURF2_16_1",      "MI0000437", "SMAD6_4_1", "SMAD2_8_1", "MI0003561", "SKP1_6_1",      "MI0003556", "LTBP1_19_1", "MI0006350", "MI0006351", "TGFBR1_6_1",      "MI0000772_1", "MI0003176", "MI0005758", "MI0003639", "MI0003641",      "CUL1_5_1", "APC_3_1", "MI0005769", "KPNB1_5_1"))      When I use this data.frame like so (after loading the package mda):               mda(Coverage~.,data=dfr,subclasses=2)$fit$coef      I get this result:                                               [,1]        [,2]     Intercept                -0.90574883 -0.42382446     min_GCs_last_3_baseszero  1.80667985  0.85786905     min_GCs_last_3_basesone  -0.02302343  0.04884007     cq                                NA          NA      What troubles me about this, is:   * Why are there coefficients for both the values of the factor min_GCs_last_3_bases?    * Why are the coefficients of cq NA? I'm willing to assume I have completely misunderstood mda. Can someone please enlighten me? Or did I just bump into a bug? (note: I've tried if the names of the columns had an influence: they don't)


title: Problem with year as a factor GLMM

text: So I need to do a GLMM, I do it this way, with package lme4 glmer(y~x1+x2+x3+year+(1|x4),family=binomial In my data, year is a factor (4 levels). So when I run my glmer, I have my result like this x1 x2 x3 year2 year3 year4 And I wish only have "year" globally. I can't transform year in a numeric variable, because it will be not correct because it is not a continuous variable. Any idea ?


title: Average minimum distance between two random vectors

text: Let $\mathbf{y_1} =\begin{bmatrix}g_1x_1 & g_2x_1 & \dots & g_Nx_1 \end{bmatrix}$ and $\mathbf{y_2} = \begin{bmatrix} f_1x_2 & f_2x_2 & \dots & f_Nx_2\end{bmatrix}$. All the elements of $\mathbf{g}=\begin{bmatrix} g_1&g_2 &\dots &g_N\end{bmatrix}$ and $\mathbf{f}= \begin{bmatrix} f_1 & f_2 & \dots f_N\end{bmatrix}$ are drawn from a rayleigh distribution. $x_1$ and $x_2$ are taken randomly from the set $\\{-1,1\\}$ with equal probability, i.e, $x_1$ can be $1$ or $-1$ with $0.5$ probability. There are four combinations for $(x_1, x_2)$: $(1,1),(1,-1),(-1,1),(-1,-1)$ which will lead to four combinations for $(\mathbf{y_1}, \mathbf{y_2})$ as well. Considering all combinations, how can I find the average $\it{minimum}$ euclidean distance between the vectors $\mathbf{y_1}$ and $\mathbf{y_2}$?


title: Bivariate Normal Distribution Mean

text: $X$ and $Y$ have bivariate normal distribution and have joint pdf   \begin{equation*} f\left( x,y\right) =a\exp \left( \frac{-1}{2}\omega \right) ,\text{where }% \omega =6x^{2}+12y^{2}-16xy-8x+24 \end{equation*} then what are the means of $X$ and $Y$? I have \begin{equation*} \frac{1}{\left( 1-\rho ^{2}\right) \sigma _{X}^{2}}=6,~\frac{1}{\left( 1-\rho ^{2}\right) \sigma _{Y}^{2}}=12~\text{and }\frac{-2\rho }{\left( 1-\rho ^{2}\right) \sigma _{X}\sigma _{Y}}=-16 \end{equation*} so i can find $\rho ,\sigma _{X}$ and $\sigma _{Y}$ using these equations. But how can i get $\mu _{X}$ and $\mu _{Y}$? this question asked in an exam and requested to be replied in 2 minutes. Is there a shortcut for this?


title: Kalman filter and Box-Cox

text: I'm interested in wind forecasting, which I have analyzed over some time by means of ARMA methods. Now I've being reading about Kalman filtering. Kalman filter is optimal when Gaussian assumption can be hold. However, wind distributions are far from normal (with Weibull the most popular). My question is, Is it correct to transform a wind speed time series by means of Box-Cox transformation before estimating the parameters of Kalman filter, so that the normality assumption holds? If it were so, why then using particle filters?


title: Complete Linkage Clustering of 3D data space coordinates

text: I have a large dataset of 3d points (XYZ coordinates) and I would like to perform hierarchical clustering using complete linkage method with Euclidean distance as clustering metric. Additionally, there are two constraints to the clustering process: The 1st constraint - maximum distance between points in the same cluster has to be lower than specified value (in this case 25), what will control the size of clusters The idea of the 2nd constraint is to use a threshold to determine whether marginal grouping of points will be identified as clusters. The clustering procedure should continue until no new clusters can be recognized - when there is no longer population of unclustered points larger than the threshold value within a spatial volume of less than maximum distance from the 1st constraint. I would like to change the number in the second constraint to find a solution where 75% to 90% of all the input points are clustered. Could anyone help me with this problem and suggest what language/software should I use? Any hints will be highly appreciated.


title: Fisher information for Poisson model

text: Assuming the data structure given below, $$ \left[\begin{array}{l|l} y & x \\\ \hline \hline y_{1} & x_{1} \\\ y_{2} & x_{2} \\\ \vdots & \vdots \\\ y_{n} & x_{n} \\\ \end{array}\right] $$ Suppose we want to fit a Poisson regression model such that $y_{i} \sim \text{Pois}(\mu_{i})$ for $i=1,2\ldots,n$. where: $\mu_{i} = e^{\beta_{0} + \beta_{1}x_{i}}$ The Fisher information can be found by: $$ I(\beta) = \sum_{i=1}^{n}\mu_{i}x_{i}x_{i}^{T} $$ Supposing we have the MLEs ($\hat{\beta}_{0}$ and $\hat{\beta}_{1}$) for $\beta_{0}$ and $\beta_{1}$, from the above, we should be able to find the Fisher information for both model parameters at the maximum: $$ I(\hat{\beta}) = \begin{bmatrix} I_{11} & I_{12} \\\ I_{21} & I_{22} \end{bmatrix} $$. I have the data and the resulting Fisher information matrix corresponding to the model above; however, when I try to code the Fisher information from the above manually, my result for the Fisher information differs from correct result. I have included the a snippet of R code. Perhaps somebody could tell me if the following code will achieve the Fisher information matrix based on the above model specification; or, more generally, if the model specification is correct.               mu <- rep(0, n) # Initialize data structure for mu vector.         for (i in 1:n) {               # Assign values to mu vector.             mu[i] <- exp(beta[1]+(beta[2]*x[i,1:2]))          }            # Initialize Fisher information matrix.     infor <- matrix(rep(0, 4), byrow=TRUE, nrow=2, ncol=2)      for (i in 1:n) {         temp_matrix <- mu[i]*x[i,1:2]%*%t(x[i,1:2]) # Hold current index value.         infor <- infor + temp_matrix # Sum to Fisher information matrix.     }          print(infor)      **NOTE:**   * Even though there is only one independent variable ($x$), I have augmented the matrix with a column of ones corresponding to the intercept parameter ($\beta_{0}$), as is usually the case with normal regression. Therefore, the independent (design) matrix is of the form: $$ X = \begin{bmatrix} 1 & x_{1,2} \\\ 1 & x_{2,2} \\\ \vdots & \vdots \\\ 1 & x_{n,2} \end{bmatrix} $$ Hence each row `i` iterates across `x[i,1:2]`. Should this column be included?   * In the preceding R code, it is assumed that data structures and values not initialized and/or assigned explicitly in the code have been previously initialized and/or assigned.   * Furthermore, the code is designed to be explicit, not efficient.


title: NIST randomness test fails in random generated bits

text: I am using the NIST randomness test suite which checks whether a bit sequence is random over 15 different tests. I have problems in verifying if this suite works fine. The situation is that if I input, for example, 10000 sequences of the same length (e.g., 1000) and generated by `rand` function of MATLAB, sometimes it passes the basic tests I use, i.e., frequency, runs, block frequency, cumulative sums. However, sometimes it fails the basic frequency test and cumulative sums test in terms of p-value of p-values but not in terms of proportion. I have two questions. The general one is how can I generate uniformly distributed random bits from MATLAB? The second is, why the problem mentioned above happens in NIST randomness suite, can anyone help me to verify if it works fine? For reference document : http://csrc.nist.gov/groups/ST/toolkit/rng/documents/SP800-22rev1a.pdf


title: Term used to describe two lines appearing closer to each other when lines are at a greater gradient

text: Forgive me if this is outside of the scope of Cross Validated, but it seemed relevant on the visualization front. As described here (page 3) there exists a phenomenon when gauging the distance between two lines on a plot to compare the shortest distance between the two lines, rather than the vertical distance. This leads to underestimation of the distance between lines at certain gradients. To illustrate, the following graphic is from the aforementioned link: ![enter image description here](http://i.stack.imgur.com/bKFfB.png) The difference on the y-axis between these two lines is uniform across all x values, and yet the difference appears to shrink. I am curious as to whether this phenomenon has an official name and if so, whether there is any kind of canonical reference paper/citation describing the phenomenon in more detail?


title: Can I repeat cross validation with a small dataset, and/or how can I improve my cross validation confidence?

text: For university we need to classify 3 cancer types and give an estimation of how well our model will perform. We received a dataset with 100 samples. We split the data up into a training and test set using stratified sampling with a ratio of 0.3 and 0.7. The resulting training set consists of 69 samples, and the test set out of 31 samples. We used 10-fold cross validation to calculate the accuracy of our models. When applying the same model on the test set for most models the accuracy on the test set is between 10-15% worse than with cross validation on the training set, except for one model where the accuracy on the test set was 2% better than during cross validation. The problem that we have now is that the two best scoring models on cross- validation are not significantly different, one has an accuracy of 88.57% +/- 12.45%, the other an accuracy of 88.00% +/- 7.92%. However, on the test set the first score 76%, and the second scores 90%. If we understood it correctly, we can't choose the second model as the best model based on the test set results, because then we would be using the test set as a training set. Instead, we would like to use repeated cross validation to improve my confidence in the cross validation results, and thereby hopefully being able to choose the best model. With the small dataset that we have, if we do repeated cross validation and take the average, would we run into the problem that the same folds would be used multiple times?


title: Sample problems on logit modeling and Bayesian methods

text: I'm looking for worked out solutions using Bayesian and/or logit analysis similar to a workbook or an annal. The worked out problems could be of any field; however, I'm interested in urban planning / transportation related fields.


title: Show that $E(x)=M'_X(0)$, where $M'_X(p)=\frac{dM_X(p)}{dp}$

text: I am going through some more exam practice questions for my stats class - Here is the complete question, it has two parts: > Let $X$ be a continuous random variable with known probability density > function (pdf) $f_X(x)$ for $-\infty < x < \infty$ and let $M_X(p)=E(e^{pX}) > $ be a real valued function of $p.$ > > (a) Write down $M_X(p)$ as an integral involving the pdf $f_X(x)$ and > $e^{pX}.$ > > (b) Show that $E(x)=M'_X(0)$, where $M'_X(p)=\frac{dM_X(p)}{dp}.$ Here is what I have so far: (i) Seems pretty straightforward from the definition of $E(g(x))$. I simply have: $$E(e^{pX}) = \int^{p}_{-\infty}e^{pX}f_X(x) dx.$$ (b) For part b, I am a little stuck. From (a), I have that $$M'_X(p)=\frac{d}{dp} M_X(p)=\frac{d}{dp}E(e^{pX})=\frac{d}{dp}\int^{p}_{-\infty}e^{pX}f_X(x) dx.$$ Is there a theorem similar to the fundamental theorem of calculus for taking the derivative of an improper integral? The way I have tried to evaluate this is: $$\frac{d}{dp}\left( \lim_{b \to -\infty}\int^{p}_{b}e^{pX}f_X(x) dx\right).$$ Now I am a little stuck. I know that the definition of $E(x) = \int^{\infty}_{-\infty}xf_X(x)dx$, but when I try to integrate the above improper integral I don't think I am getting anywhere but making it more complicated. Am I even on the right track? I feel like I am making it more complicated than it is supposed to be. Any help would be greatly appreciated!


title: Where can I find information about using SPSS for EFA and CFA? Is PCA (two samples) and reliability sufficient for scale development?

text: ### Context: I am in the process of developing a scale for my thesis. My advisor has guided me to using SPSS PCA to complete my analyses. Initially we reduced my scale to 3 factors (her insistence), and we collected a second sample with the 3 factor scale. I have had issues with the 3 factors loading appropriately. Participants respond differently to the third factor based on age and gender, and the loadings are coming in at .45 and then .42 in another factor. When I consider the two factor model (15 item scale vs. 10 item scale), all loading is .65 or higher, and responding is not effected by gender or age. ### Questions:   * What do you suggest I read in order to learn more about using SPSS PCA for scale development?   * Should I use SEM?   * Are there additional analyses I can conduct on this scale to determine the fit/factors?


title: Euclidean minimum spanning tree in R

text: I have recently been trying to build a euclidean minimum spanning tree. However, is it at all possible to plot the trees on scales, or do they always have to be scale-free? Thank you!


title: Compare modeled (fitted) paired data to actual data in forecasting problem (Excel sheet included)

text: In my endeavor to learn forecasting and improve my statistical knowledge, I've decided to forecast the population of a certain area. I've attached the excel file that I used. I know I'm using excel which is not very common around here. I'm familiar with R but prefer excel for fast analysis such as this. The collected data is for 1921 to 2009. The forecast period is between 2010 and 2019. I did:   * Arima over all the data   * Arima over the last 10 years or so of the data   * Linear regression over the whole time period   * Linear regression over the last 10 years   * Fitted model using Eureqa software   * Average of above models without the overall linear regression. The forecast worksheet includes all the raw data, while the forecasting simplified worksheet includes the simplified chart to better view the results. The question, what are some of the ways I can validate and/or calibrate my forecasts, numbers, assumptions, and improve my analysis from an amateur level to something that is more professional. Thanks Excel File


title: Use of Signal to Noise Ratio to compare Ridge Regression Estimators

text: On p170 of Gruber's "Regression Estimators" 2010, he shows that the Liu estimator has a smaller bias than the ordinary ridge regressor but also a larger variance. So when the signal to noise ratio is compared i.e. (Bias)^2/Variance the Liu estimator has a lower signal to noise ratio. However this seems odd to me, choosing the signal to noise ratio calculated this way will encourage you towards picking models with larger biases assuming they have same variance or smaller than the competing model? Is this the normal way to look at signal to noise ratios when comparing ridge regression estimators? I mean if you have an OLS estimate (and assuming no multicollinearity) there would be no bias and it would have the lowest variance? Yet it would have a signal to noise ratio of zero??


title: Information Geometry for Ridge Regression

text: A question has already been asked about the usefulness of informational geometry. Using information geometry to define distances and volumes…useful? I am reading the book "regression estimators" by Gruber 2010 where he uses this technique to compare Ridge Regressors, however he concentrates on deriving the mathematical results without giving any empirical results/examples. I am eager to read empirical type papers where these informational geometry techniques are used to compare ridge regressors on real data. It seems from reading the book that the regressor with the highest variance will have the shorter distance form OLS. So it would seem that the regressor with the biggest distance should be selected but this is not explicitly stated in the book. I'd also be very interested in hearing about alternate methods to compare ridge regressors and the relative pros/cons of those approaches ideally compared to the informational geometry approach.


title: Multivariate Time Series

text: I am trying to learn multivariate time series using R. I have two time series and I want to see if I could use one of those to predict the other one, and after that check if the model holds or there is no relationship between both variables. Does anyone have an example that I could use to understand this concept? Thank you and kind regards


title: Fitting time series models

text: I have two sets of time series data for 36 months. It contains seasonal trends with a 12-month cycle.   1. How to determine whether it is a good model? The smaller the AIC, the better the model?   2. Do I need to do any transformation before using auto.arima? As I find in google that auto.arima has already dealt with seasonal trends.   3. Do I need to re-model if the residuals of the forecast not following $\mathcal N(0,1)$?   4. Do I need more data to do the forecasting? (now only 3 sets of data to do forecasting, due to seasonal trends). Thanks.


title: Density of multivariate normal given linear condition

text: Given a variable $X \sim N(\mu, \Sigma)$ What is the density of $X$ given $HX = d$,   $H$ and $d$ both constant I.e. we observe $X$ with rank-deficient linear operator $H$ and obtain some value $d$. Given this knowledge how do we update our knowledge/distribution of $X$? This is just one of the steps in the Kalman filter. I'm having trouble separating it out. Ideally an answer would be analytic forms for both $\mu$ and $\Sigma$ after the update plus a rationale/derivation or link to such a derivation.


title: Linear regression parameters question

text: Are the slope and intercept of a simple linear regression model always normally distributed? Is there ever a difference between the distribution of the estimated slope and intercept and the actual ones? I have only just begun learning about the subject but I am still not clear on the details. A final question: is the least squares method the same as linear regression in that it gives information like the $R^2$? Thanks!


title: Role of coefficients in model selection for logistic regression

text: I have a model that I am using to predict mortality and it gives me an AUC of 0.799. The `R` code that I am using would look something like this:               fit <- glm(mortality ~ Age + Sex + ASA + Wndclass, data=sepsis, family=binomial)      I am calculating my AUC using `ROCR`:               library(ROCR)     fit.pred <- predict(fit)     sepsis.pred <- prediction(fit.pred, mortality[,1])     sepsis.perf <- performance(sepsis.pred, "auc")       Calling `sepsis.perf` gives me:               An object of class "performance"     Slot "x.name":     [1] "None"          Slot "y.name":     [1] "Area under the ROC curve"          Slot "alpha.name":     [1] "none"          Slot "x.values":     list()          Slot "y.values":     [[1]]     [1] 0.7990628               Slot "alpha.values":     list()      So far I have conluded that this model isn't too bad as a classifier. However when I look at the coefficients, especially the 95% CI, of the logistic regression model I get some strange results:               (Intercept)                   8.248484e-03  3.394696e-02     SEXfemale                     9.037957e-01  1.325291e+00     SEXNULL                      3.864079e-204 5.780619e+192     Age                           1.026437e+00  1.043482e+00     ASACLAS2-Mild Disturb        1.408200e-220 3.478857e+229     ASACLAS3-Severe Disturb      1.769693e-220 4.367441e+229     ASACLAS4-Life Threat         3.596686e-220 8.876357e+229     ASACLAS5-Moribund            8.561098e-220 2.113924e+230     ASACLASNULL                  2.702457e-219 6.747396e+230     ASACLASNone assigned         1.513924e-219 3.793387e+230     WNDCLAS2-Clean/Contaminated   9.281083e-01  1.897956e+00     WNDCLAS3-Contaminated         8.968408e-01  1.604651e+00     WNDCLAS4-Dirty/Infected       5.114258e-01  8.718848e-01      What bothers me here is the 95% CI for ASA class, which is ridiculous. However when I take it out of my model, the AUC drops to 0.78, and clinically, the ASA (American Society of Anesthesiologists) classification is often used to stratify patients into low to high risk for mortality in the perioperative period. So my question is, given these results of the 95% CI, should I exclude ASA class? If not, how should I interpret ASA class in view of my results? Thanks for the help!