I am trying to calculate variances via np.std(array,ddof = 0). The problem emerges if I happen to have a long delta array, i.e., all values in the array are the same. Instead of returning std = 0, it gives some small value which in turn causes further estimation errors. The mean is returned correctly... Example:  np.std([0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1],ddof = 0)   gives 1.80411241502e-16  but   np.std([0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1],ddof = 0)   gives std = 0  Is there a way to overcome this except for now checking for uniqueness of data for every iteration without calculating std at all?  Thanks  P.S. Following marking as duplication of , copy-pasting the reply by @kxr on why it's a different question:   "The current duplicate marking is wrong. Its not just about simple float comparison, but about internal aggregation of small errors for near-zero outcome by using the np.std on long arrays - as the questioner indicated extra. Compare e.g. &gt;&gt;&gt; np.std([0.1, 0.1, 0.1, 0.1, 0.1, 0.1]*200000) -&gt; 2.0808632594793153e-12 . So he can e.g. solve by: &gt;&gt;&gt; mean = a.mean(); xmean = round(mean, int(-log10(mean)+9)); std = np.sqrt(((a - xmean) ** 2).sum()/ a.size)"  The problem certainly starts with floating representation but it does not stop there. @kxr - I appreciate the comment and the example
27
Consider the following code:   0.1 + 0.2 == 0.3  -&gt;  false   0.1 + 0.2         -&gt;  0.30000000000000004   Why do these inaccuracies happen?
46479
How can I show that the variance of local polynomial regression is increasing with the degree of the polynomial (Exercise 6.3 in Elements of Statistical Learning, second edition)?  This question has been asked  but the answer just states it follows easliy.  More precisely, we consider $y_{i}=f(x_{i})+\epsilon_{i}$ with $\epsilon_{i}$ being independent with standard deviation $\sigma.$  The estimator is given by  $$ \hat{f}(x_{0})=\left(\begin{array}{ccccc} 1 &amp; x_{0} &amp; x_{0}^{2} &amp; \dots &amp; x_{0}^{d}\end{array}\right)\left(\begin{array}{c} \alpha\\ \beta_{1}\\ \vdots\\ \beta_{d} \end{array}\right) $$ for $\alpha,\beta_{1},\dots,\beta_{d}$ solving the following weighted least squares problem $$ \min\left(y_{d}-\underbrace{\left(\begin{array}{ccccc} 1 &amp; x_{1} &amp; x_{1}^{2} &amp; \dots &amp; x_{1}^{d}\\ \vdots\\ 1 &amp;  &amp;  &amp;  &amp; x_{n}^{d} \end{array}\right)}_{X}\left(\begin{array}{c} \alpha\\ \beta_{1}\\ \vdots\\ \beta_{d} \end{array}\right)\right)^{t}W\left(y-\left(\begin{array}{ccccc} 1 &amp; x_{1} &amp; x_{1}^{2} &amp; \dots &amp; x_{1}^{d}\\ \vdots\\ 1 &amp;  &amp;  &amp;  &amp; x_{n}^{d} \end{array}\right)\left(\begin{array}{c} \alpha\\ \beta_{1}\\ \vdots\\ \beta_{d} \end{array}\right)\right) $$ for $W=\text{diag}\left(K(x_{0},x_{i})\right)_{i=1\dots n}$ with $K$ being the regression kernel. The solution to the weighted least squares problem can be written as $$ \left(\begin{array}{cccc} \alpha &amp; \beta_{1} &amp; \dots &amp; \beta_{d}\end{array}\right)=\left(X^{t}WX\right)^{-1}X^{t}WY. $$ Thus, for $l(x_{0})=\left(\begin{array}{ccccc} 1 &amp; x_{0} &amp; x_{0}^{2} &amp; \dots &amp; x_{0}^{d}\end{array}\right)\left(X^{t}WX\right)^{-1}X^{t}W$ we obtain $$ \hat{f}(x_{0})=l(x_{0})Y $$ implying that  $$ \text{Var }\hat{f}(x_{0})=\sigma^{2}\left\Vert l(x_{0})\right\Vert ^{2}=\left(\begin{array}{ccccc} 1 &amp; x_{0} &amp; x_{0}^{2} &amp; \dots &amp; x_{0}^{d}\end{array}\right)\left(X^{t}WX\right)^{-1}X^{t}W^{2}X\left(X^{t}WX\right)^{-1}\left(\begin{array}{ccccc} 1 &amp; x_{0} &amp; x_{0}^{2} &amp; \dots &amp; x_{0}^{d}\end{array}\right)^{t}. $$ My approach: An induction using the formula for the inverse of a block matrix but I did not succeed.  The paper  by D. Ruppert and M. P. Wand derives an asymptotic expression for the variance for $n\rightarrow\infty$ in Theorem 4.1 but it is not clear that is increasing in the degree.
101154
Suppose I have a discrete set of (finite) data (values are only positive integers) (which always remains discrete whenever the observation/survey is taken), divided into some identical categories . Now, I want to calculate the mean values of some random variable for each category, hence getting a sequence of mean values. Considering mean as a random variable (let's call it $M$), it's taking decimal values also (obviously possible). But I am not sure if $M$ is a discrete or a continuous random variable.   My intuition suggest me that it should be discrete. Since, my original data from which I calculated the sets of mean is discrete so the mean can take finitely many values which I suppose is not the case with a continuous random variable where the mean can take infinitely many values within an interval. So is my intuition right or wrong?  Is it true that a discrete random variable takes finitely many values and continuous R.V. can take infinitely many possible values?
27333
In a scientific paper, I've seen the following $$\frac{\delta K^{-1}}{\delta p} = -K^{-1}\frac{\delta K}{\delta p}K^{-1}$$ where $K$ is a $n \times n$ matrix that depends on $p$. In my calculations I would have done the following $$\frac{\delta K^{-1}}{\delta p} = -K^{-2}\frac{\delta K}{\delta p}=-K^{-T}K^{-1}\frac{\delta K}{\delta p}$$ Is my calculation wrong? Note: I think $K$ is symmetric.
1061
How can I prove that a random variable taking values in $[0,1]$ has variance no larger than $\frac{1}{4}$? If it matters, discrete and continuous proofs are both welcome.
10131
I need to compute the mean and variance of the truncated normal distribution. For simplicity, let us focus on a standard normal, since the general case can be reduced to this. The PDF is given by:  $$f\left(x;a,b\right)=\begin{cases} \frac{\phi\left(x\right)}{\Phi\left(b\right)-\Phi\left(a\right)} &amp; a\le x\le b\\ 0 &amp; \text{otherwise} \end{cases}$$  where  $$\phi\left(x\right)=\frac{1}{\sqrt{2\pi}}\mathrm{e}^{-x^{2}/2},\quad\Phi\left(x\right)=\frac{1}{2}\left(1+\mathrm{erf}\left(x/\sqrt{2}\right)\right)$$  The analytical formulas for the mean and variance of $f$ are:  $$\left\langle x\right\rangle  =\frac{\phi\left(a\right)-\phi\left(b\right)}{\Phi\left(b\right)-\Phi\left(a\right)}$$  $$\mathrm{var}\,x =1+\frac{a\phi\left(a\right)-\beta\phi\left(b\right)}{\Phi\left(b\right)-\Phi\left(a\right)}-\left(\frac{\phi\left(a\right)-\phi\left(b\right)}{\Phi\left(b\right)-\Phi\left(a\right)}\right)^{2}$$  However, if you try to evaluate those expressions numerically, when $a,b$ are sufficiently away from the mean, you will encounter multiple numerical issues. I have tried Mathematica, Python scipy, Julia Distributions.jl, and they all give NaN.  For example, this is what happens in Mathematica:    Has anyone encountered this issue before? Anyone can suggest a solution, or an alternative package?
88239
The mean absolute deviation is:   $$\dfrac{\sum_{i=1}^{n}|x_i-\bar x|}{n}$$   The variance is: $$\dfrac{\sum_{i=1}^{n}(x_i-\bar x)^2}{n-1}$$   So the mean deviation and the variance are measuring the same thing, yet variance requires squaring the difference. Why? Squaring always gives a non-negative value, but the absolute value is also a non-negative value. Why isn't it $|x_i-\bar x|^2$, then? Squaring just enlarges, why do we need to do this?  A similar question is , but mine is a little different. Thanks.
50690
I am writing an algorithm where, given a model, I compute likelihoods for a list of datasets and then need to normalize (to probability) each one of the likelihood. So something like [0.00043, 0.00004, 0.00321] might be converted to may be like [0.2, 0.03, 0.77].   My problem is that the log likelihoods, I am working with, are quite small (for instance, in log space, the values are like -269647.432, -231444.981 etc). In my C++ code, when I try to add two of them (by taking their exponent) I get an answer of "Inf". I tried to add them in log-space , but again stumbled upon the same problem.  Can anybody share his/her expert opinion on this?
37989
I have following simple X and Y vectors:   &gt; X [1] 1.000 0.063 0.031 0.012 0.005 0.000 &gt; Y [1] 1.000 1.000 1.000 0.961 0.884 0.000 &gt;  &gt; plot(X,Y)     I want to do regression using log of X. To avoid getting log(0), I try to put +1 or +0.1 or +0.00001 or +0.000000000000001 :   &gt; summary(lm(Y~log(X))) Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :    NA/NaN/Inf in 'x' &gt; summary(lm(Y~log(1+X)))  Call: lm(formula = Y ~ log(1 + X))  Residuals:        1        2        3        4        5        6  -0.03429  0.22189  0.23428  0.20282  0.12864 -0.75334   Coefficients:             Estimate Std. Error t value Pr(&gt;|t|)   (Intercept)   0.7533     0.1976   3.812   0.0189 * log(1 + X)    0.4053     0.6949   0.583   0.5910   --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Residual standard error: 0.4273 on 4 degrees of freedom Multiple R-squared:  0.07838,   Adjusted R-squared:  -0.152  F-statistic: 0.3402 on 1 and 4 DF,  p-value: 0.591  &gt; summary(lm(Y~log(0.1+X)))  Call: lm(formula = Y ~ log(0.1 + X))  Residuals:        1        2        3        4        5        6  -0.08099  0.20207  0.23447  0.21870  0.15126 -0.72550   Coefficients:              Estimate Std. Error t value Pr(&gt;|t|)   (Intercept)    1.0669     0.3941   2.707   0.0537 . log(0.1 + X)   0.1482     0.2030   0.730   0.5058   --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Residual standard error: 0.4182 on 4 degrees of freedom Multiple R-squared:  0.1176,    Adjusted R-squared:  -0.103  F-statistic: 0.5331 on 1 and 4 DF,  p-value: 0.5058  &gt; summary(lm(Y~log(0.00001+X)))  Call: lm(formula = Y ~ log(1e-05 + X))  Residuals:        1        2        3        4        5        6  -0.24072  0.02087  0.08796  0.13872  0.14445 -0.15128   Coefficients:                Estimate Std. Error t value Pr(&gt;|t|)     (Intercept)     1.24072    0.12046  10.300 0.000501 *** log(1e-05 + X)  0.09463    0.02087   4.534 0.010547 *   --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Residual standard error: 0.1797 on 4 degrees of freedom Multiple R-squared:  0.8371,    Adjusted R-squared:  0.7964  F-statistic: 20.56 on 1 and 4 DF,  p-value: 0.01055  &gt;  &gt; summary(lm(Y~log(0.000000000000001+X)))  Call: lm(formula = Y ~ log(1e-15 + X))  Residuals:         1         2         3         4         5         6  -0.065506  0.019244  0.040983  0.031077 -0.019085 -0.006714   Coefficients:                Estimate Std. Error t value Pr(&gt;|t|)     (Intercept)     1.06551    0.02202   48.38 1.09e-06 *** log(1e-15 + X)  0.03066    0.00152   20.17 3.57e-05 *** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Residual standard error: 0.04392 on 4 degrees of freedom Multiple R-squared:  0.9903,    Adjusted R-squared:  0.9878  F-statistic: 406.9 on 1 and 4 DF,  p-value: 3.565e-05   The output is different in all cases. What is the correct value to put to avoid log(0) in regression? What is the correct method for such situations.  Edit: my main aim is to improve prediction of the regression model by adding log term, i.e.: lm(Y ~ X + log(X))
34219
In the vein of , I'm now at a roadblock on question 3 of this sheet:    (note: it's not my intention to ask every question I get stuck on here, merely the ones which have interesting general results; it just so happens that the two I've struggled with so far fit this criteria)  It goes as follows:  Let $N$ be a non-negative integer-valued random variable with mean $\mu_1$ and variance $\sigma_1^2$, and let $X_1, X_2, ...$ be random variables, each with mean $\mu_2$ and variance $\sigma_2^2$; furthermore, assume that $N, X_1, X_2,  . . .$ are independent. Without using generating functions, calculate the mean and variance of the random variable $S_N = X_1 + ... + X_N$ (when $N=0$ interpret $S_N$ as $0$).  I have the answer to calculating the mean, which I'll write up as an answer for future reference. Here's where I've got to with the variance bit:  $\mathbb{E}(S_N^2) = 0.P(N=0) + \sum_{r=1}^\infty\mathbb{E}((\sum_{i=1}^rX_i)^2).P(N=r)$  For $r=1$, we have the inner expectation as: $\mathbb{E}(X_1^2)$ which is equal to $\sigma_2^2+\mu_2^2$.  For $r &gt; 1$, we have the inner expectation equal to: $r(\sigma_2^2 + 2\mu_2^2)$  So $\mathbb{E}(S_N^2) = (\sigma_2^2+\mu_2^2)P(N=1) + (\sigma_2^2 + 2\mu_2^2)\sum_{r=2}^\infty rP(N=r)$  which is almost in the form where I can use the same trick as in calculating the mean, but not quite.  Any help is much appreciated.
14594
Say that the raw data is $N$-dimensional, where $N$ is a large positive integer. If we apply PCA to the dataset, and compute the reconstruction error (in $\ell^2$-norm) using the first $d \leq N$ components, will the reconstruction error monotonically decrease to $0$ as $d \rightarrow N$?
94976
I've been testing PCA via SVD to decompose a simple time series data matrix, $X$.  I have two signals $x_1(t)$ and $x_2(t)$ in a data matrix where $M$ rows represents each timepoint sample and each column represents $x_1$ and $x_2$.  The mean signal, $\hat{x}$, is defined as the mean along the row axis (average of $x_1$ and $x_2$ along each timepoint).  I normalize each column of $X$ by subtracting its mean and dividing by the standard deviation.  When I use [U S V] = svd(Xz) in matlab, regardless of how the variables are distributed (whether they are correlated or uncorrelated), one of the columns of the right singular matrix, V, always points in the same direction (to a multiplicative constant) as the mean vector $(1/2, 1/2)$.  But when add an additional third time series, this is never the case (where the mean vector is $(1/3, 1/3, 1/3)$.  Because I normalize the standard deviation for each vector, it does make sense that the direction of most variance given by PCA would be the diagonal 45 degree line.  But if both variables $x_1$ and $x_2$ are independent gaussians, couldn't the PCA direction be any direction since the distribution is radially symmetric?  MATLAB Code:  s = RandStream('mcg16807', 'Seed', 0); RandStream.setDefaultStream(s);  G = zeros(1000,2); G(:,1) = 40*randn(1000,1)-100; G(:,2) = tan(G(:,1)) + randn(1000,1); X = G - repmat(mean(G),[size(G,1), 1]); Xz = X./repmat(std(G),[size(G,1), 1]); [U S V] = svd(Xz); G_mean = mean(Xz,2); corrcoef(G_mean, U(:,1)) corrcoef(G_mean, U(:,2))
99542
I am running a GLM, using the function glm.nb (pscl package) trying to figure how what could influence a particular trait in several locations and years. The output as follow (with slight modification and removing the things beyond this question)            Estimate Std. Error z value Pr(&gt;|z|)               state2        -0.163190   0.807307  -0.202  0.00081               state3        -0.530588   0.783023  -0.678  0.00302            state4         0.942953   0.737328   1.279  0.00094               year2         -1.759102   0.733214  -2.399  0.01643            year3         -0.633870   0.662903  -0.956  0.33897       I understand the output as such: compared to year1, year 2 is statistically different, but year3 is not different from year1. but what about year2 and year3?  For the state factor, state2, state3 and state4 are all statistically different than state1. so I need a post-hoc test between state2 state3 and state4, am I right?  thanks
42890
We have ground-truth data $\mathbf{x}^* = (0,0)^T \in \mathbb{R}^{2}$. Furthermore, we have $N$ measurements $\mathbf{x}_i \in \mathbb{R}^{2}, i\in \{1,\ldots,N\}$. We calculate the $N$ error vectors $\mathbf{e}_i = \mathbf{x}_i - \mathbf{x}^*$. We assume that the two components of $\mathbf{e}_i = (e_{x,i}, e_{y,i})^T $ are independent of each other and normally distributed with equal mean $\mathbf{\mu} = (0,0)^T$ and variance $\sigma^2$. The absolute errors $e_i = \sqrt {e_{x,i}^2 + {e}_{y,i}^2} $ follow a Rayleigh distribution, is that correct?  The probability density function (PDF) of the Rayleigh function is zero at a value of zero. Why is this the case (intuitively)? I would assume that our errors should most likely be 0 (i.e. PDF should be highest at 0)...
117149
I am trying to do the following exercise:  Let $\xi_i,\dots, \xi_n$ be independent Bernoulli random variables such that $$ P\{\xi_i=0\}=1-\lambda_i \Delta\\ P\{\xi_i=1\}=\lambda_i \Delta $$ where $\Delta$ is a small number, $\Delta&gt;0,\lambda_i&gt;0$.  Show that $$ P\{\xi_1+\dots+\xi_n=1\}=(\sum_{i=1}^n\lambda_i)\Delta + O(\Delta^2) $$    I already wrote out (as already shown in part of @Marcin Malogrosz answer): $$ P\{\xi_1+\dots+\xi_n) = \sum_i \lambda_i \Delta \prod_{i\neq j} (1-\lambda_j\Delta) \\ \Bigl|\sum_i \lambda_i \Delta \prod_{i\neq j} (1-\lambda_j\Delta) - \sum_i \lambda_i \Delta\Bigr| = \sum_i \lambda_i \Delta \Bigl(1-\prod_{i\neq j} (1-\lambda_j\Delta) \Bigr) $$ But I am stuck here    I am able to show a simplified version with $P\{\xi_i=1\}=\lambda_i\Delta=\lambda\Delta$ (when the $n$ random variables are also identical) as following: $$ P\{\xi_1+\dots+\xi_n=1\}=\sum_i \lambda\Delta (1-\lambda\Delta)^{n-1} $$ using Taylor approximation when $\Delta$ is close to $0$ as $$ \lambda\Delta(1-\lambda\Delta)^{n-1}\Bigr|_{\Delta=0}+\frac{\Delta}{1!}\Bigl\{\lambda(1-\lambda\Delta)^{n-1}-\lambda^2\Delta^2(n-1)(1-\lambda\Delta)^{n-2}\Bigl\}\Bigr|_{\Delta=0} + O(\Delta^2)\\ =\lambda\Delta+O(\Delta^2) $$ and then sum over $i$.    My problem is to show:  $$ \Bigl(1-\prod_{i\neq j} (1-\lambda_j\Delta) \Bigr)=O(\Delta) $$  Thank you in advance
95976
Unbiased weighted variance was already addressed  and  but there still seems to be a surprising amount of confusion.  There appears to be a consensus toward the formula presented in the first link as well as in the .  This also looks like the formula used by R, Mathematica, and GSL (but not MATLAB).  However, the Wikipedia article also contains the following line which looks like a great sanity check for a weighted variance implementation:     For example, if values {2,2,4,5,5,5} are drawn from the same   distribution, then we can treat this set as an unweighted sample, or   we can treat it as the weighted sample {2,4,5} with corresponding   weights {2,1,3}, and we should get the same results.   My calculations give the value of 2.1667 for variance of the original values and 2.9545 for the weighted variance.  Should I really expect them to be the same?  Why or why not?
16546
Suppose, $X\sim N(\mu,\sigma^2)$.   Can anyone help in finding the following : $\text{Var } \bigg(\Phi\big(\frac{X + c}{d}\big) \bigg)$ ? Here, c and d are positive.  Here, $\Phi(x)$ is the "Cumulative Distribution Function" of the above-mentioned normal distribution.  Thanking you in advance.