Some Statistics Knowledge
What's the essence of probability? There are two views: Probability is related to sampling assumptions. Example: Bertrand Paradox : there are many ways to randoly select a chord on a circle, with different proability densities of chord. A distribution tells how likely a random variable will be what value: Independent means that two random variables don't affect each other. Knowing one doesn't affect the distribution of other. But there are dependent random variables that, when you know one, the distribution of another changes. P ( X = x ) P(X=x) P ( X = x ) means the probability of random variable X X X take value x x x . It can also be written as P X ( x ) P_X(x) P X ( x ) or P ( X ) P(X) P ( X ) . Sometimes the probability density function f f f is used to represent a distribution. A joint distribution tells how likely a combination of multiple variables will be what value. For a joint distribution of X and Y, each outcome is a pair of X and Y, denoted ( X , Y ) (X, Y) ( X , Y ) . If X and Y are independent, then P ( X = x , Y = y ) = P ( ( X , Y ) = ( x , y ) ) = P ( X = x ) ⋅ P ( Y = y ) P(X=x,Y=y)=P((X,Y)=(x,y))=P(X=x) \cdot P(Y=y) P ( X = x , Y = y ) = P (( X , Y ) = ( x , y )) = P ( X = x ) ⋅ P ( Y = y ) . For a joint distribution of ( X , Y ) (X, Y) ( X , Y ) , if we only care about X, then the distribution of X is called marginal distribution. You can only add probability when two events are mutually exclusive. You can only multiply probability when two events are independent, or multiplying a conditional probability with the condition's probability. P ( E ∣ C ) P(E \vert C) P ( E ∣ C ) means the probability of E E E happening if C C C happens. E and C both happen ) P ( E ∩ C ) = P ( E ∣ C ) ⋅ P ( C ) If E and C are independent, then P ( E ∩ C ) = P ( E ) P ( C ) P(E \cap C) = P(E)P(C) P ( E ∩ C ) = P ( E ) P ( C ) , then P ( E ∣ C ) = P ( E ) P(E \vert C)=P(E) P ( E ∣ C ) = P ( E ) . For example, there is a medical testing method of a disease. The test result can be positive (indicate having diesase) or negative. But that test is not always accurate. There are two random variables: whether test result is positive, whther the person actually has disease. This is a joint distribution. The 4 cases: a , b , c , d a, b, c, d a , b , c , d are four possibilities. a + b + c + d = 1 a + b + c + d = 1 a + b + c + d = 1 . For that distribution, there are two marginal distributions. If we only care about whether the person actually has disease and ignore the test result, then the marginal distribution is: Similarily there is also a marginal distribution of whether the test result is positive. False negative rate is P ( Test is negative ∣ Actually has disease ) P(\text{Test is negative } \vert \text{ Actually has disease}) P ( Test is negative ∣ Actually has disease ) , it means the rate of negative test when actually having disease. And false positive rate is P ( Test is positive ∣ Actually doesn’t have disease ) P(\text{Test is positive } \vert \text{ Actually doesn't have disease}) P ( Test is positive ∣ Actually doesn’t have disease ) . Some people may intuitively think false negative rate means P ( Test result is false ∣ Test is negative ) P(\text{Test result is false } \vert \text{ Test is negative}) P ( Test result is false ∣ Test is negative ) , which equals P ( Actually has disease ∣ Test is negative ) P(\text{Actually has disease } \vert \text{ Test is negative}) P ( Actually has disease ∣ Test is negative ) , which equals b b + d \frac{b}{b+d} b + d b . But that's not the official definition of false negative. Bayes theorem allow "reversing" P ( A ∣ B ) P(A \vert B) P ( A ∣ B ) as P ( B ∣ A ) P(B \vert A) P ( B ∣ A ) : The theoretical mean is the "weighted average" of all possible cases using theoretical probabilities. E [ X ] E[X] E [ X ] denotes the theoretical mean of random variable X X X , also called the expected value of X X X . It's also often denoted as μ \mu μ . For discrete case, E [ X ] E[X] E [ X ] is calculated by summing all theoretically possible values multiply by their theoretical probability. The mean for discrete case: x ∑ x ⋅ P ( X = x ) probability of that case The mean for continuous case: Some rules related to mean: (The constant k k k doesn't necessarily need to be globally constant. It just need to be a certain value that's not affected by the random outcome. It just need to be "constant in context".) Another important rule is that, if X X X and Y Y Y are independent, then Because when X X X and Y Y Y are independent, P ( X = x i , Y = y j ) = P ( X = x i ) ⋅ P ( Y = y j ) P(X=x_i, Y=y_j) = P(X=x_i) \cdot P(Y=y_j) P ( X = x i , Y = y j ) = P ( X = x i ) ⋅ P ( Y = y j ) , then: Note that E [ X + Y ] = E [ X ] + E [ Y ] E[X+Y]=E[X]+E[Y] E [ X + Y ] = E [ X ] + E [ Y ] always work regardless of independence, but E [ X Y ] = E [ X ] E [ Y ] E[XY]=E[X]E[Y] E [ X Y ] = E [ X ] E [ Y ] requires independence. For a sum, the common factor that's not related to sum index can be extraced out. So: j ∑ ( irrelevant to j f ( i ) ⋅ g ( j )) = i ∑ f ( i ) irrelevant to i j ∑ g ( j ) = ( i ∑ f ( i ) ) ( j ∑ g ( j ) ) Then: (That's for the discrete case. Continuous case is similar.) If we have n n n samples of X X X , denoted X 1 , X 2 , . . . X n X_1, X_2, ... X_n X 1 , X 2 , ... X n , each sample is a random variable, and each sample is independent to each other , and each sample are taken from the same distribution (independently and identically distributed, i.i.d ), then we can estimate the theoretical mean by calculating the average. The estimated mean is denoted as μ ^ \hat{\mu} μ ^ (Mu hat): Hat ^ \hat{} ^ means it's an empirical value calculated from samples, not the theoretical value. Some important clarifications: The mean of estimated mean equals the theoretical mean. Note that if the samples are not independent to each other, or they are taken from different distributions, then the estimation will be possibly biased. The theoretical variance, Var [ X ] \text{Var}[X] Var [ X ] , also denoted as σ 2 \sigma ^2 σ 2 , measures how "spread out" the samples are. If k k k is a constant: Standard deviation ( stdev ) σ \sigma σ is the square root of variance. Multiplying a random variable by a constant also multiplies the standard deviation. The covariance Cov [ X , Y ] \text{Cov}[X, Y] Cov [ X , Y ] measures the "joint variability" of two random variables X X X and Y Y Y . Some rules related to variance: If X X X and Y Y Y are indepdenent, as previouly mentioned E [ X Y ] = E [ X ] ⋅ E [ Y ] E[XY]=E[X]\cdot E[Y] E [ X Y ] = E [ X ] ⋅ E [ Y ] , then so Var [ X + Y ] = Var [ X ] + Var [ Y ] \text{Var}[X + Y]= \text{Var}[X] + \text{Var}[Y] Var [ X + Y ] = Var [ X ] + Var [ Y ] The mean is sometimes also called location. The variance is sometimes called dispersion. If we have some i.i.d samples but don't know the theoretical variance, how to estimate the variance? If we know the theoretical mean, then it's simple: However, the theoretical mean is different to the estimated mean. If we don't know the theoretical mean and use the estimated mean, it will be biased, and we need to divide n − 1 n-1 n − 1 instead of n n n to avoid bias: This is called Bessel's correction. note that the more i.i.d samples you have, the smaller the bias, so if you have many i.i.d samples, then the bias doesn't matter in practice. Originally, n samples have n degrees of freedom. If we keep the estimated mean fixed, then it will only have n-1 degrees of freedom. That's an intuitive explanation of the correction. The exact dedution of that correction is tricky: Firstly, the estimated mean itself also has variance As each sample is independent to other samples. As previously mentioned, if X X X and Y Y Y are independent, adding the variable also adds the variance: Var [ X + Y ] = Var [ X ] + Var [ Y ] \text{Var}[X + Y]= \text{Var}[X] + \text{Var}[Y] Var [ X + Y ] = Var [ X ] + Var [ Y ] . So: As previously mentioned E [ μ ^ ] = μ E[\hat{\mu}] = \mu E [ μ ^ ] = μ , then Var [ μ ^ ] = E [ ( μ ^ − E [ μ ^ ] ) 2 ] = E [ ( μ ^ − μ ) 2 ] = σ 2 n \text{Var}[\hat{\mu}] = E[(\hat{\mu} - E[\hat{\mu}])^2] = E[(\hat{\mu} - \mu)^2] = \frac{\sigma^2}{n} Var [ μ ^ ] = E [( μ ^ − E [ μ ^ ] ) 2 ] = E [( μ ^ − μ ) 2 ] = n σ 2 . This will be used later. A trick is to rewrite X i − μ ^ X_i - \hat{\mu} X i − μ ^ to ( X i − μ ) − ( μ ^ − μ ) (X_i - \mu) - (\hat{\mu} - \mu) ( X i − μ ) − ( μ ^ − μ ) and then expand: Then take mean of two sides: There are now three terms. The first one equals n σ 2 n\sigma^2 n σ 2 : So the second one becomes Now the above three things become E [ ( μ ^ − μ ) 2 ] E[(\hat{\mu}-\mu)^2] E [( μ ^ − μ ) 2 ] is also Var [ μ ^ ] \text{Var}[\hat{\mu}] Var [ μ ^ ] . As previously mentioned, it equals σ 2 n \frac{\sigma^2}{n} n σ 2 , so For a random variable X X X , if we know its mean μ \mu μ and standard deviation σ \sigma σ then we can "standardize" it so that its mean become 0 and standard deviation become 1: That's called Z-score or standard score. Often the theoretical mean and theoretical standard deviation is unknown, so z score is computed using sample mean and sample stdev: In deep learning, normalization uses Z score: Note that in layer normalization and batch normalization, the variance usually divides by n n n instead of n − 1 n-1 n − 1 . Computing Z-score for a vector can also be seen as a projection: (or n − 1 \sqrt{n-1} n − 1 ): σ 2 = 1 n ( y ) 2 \boldsymbol\sigma^2 = \frac 1 n (\boldsymbol y)^2 σ 2 = n 1 ( y ) 2 , σ = 1 n ∣ y ∣ \boldsymbol\sigma = \frac 1 {\sqrt{n}} \vert \boldsymbol y \vert σ = n 1 ∣ y ∣ . Dividing by standard deviation can be seen as projecting it onto unit sphere then multiply by n \sqrt n n (or n − 1 \sqrt{n-1} n − 1 ). So computing Z-score can be seen as firstly projecting onto a hyper-plane that's orthogonal to 1 \boldsymbol 1 1 and then projecting onto unit sphere then multiply by n \sqrt n n (or n − 1 \sqrt{n-1} n − 1 ). Skewness Skewness measures which side has more extreme values. A large positive skew means there is a fat tail on positive side (may have positive extreme values). A large negative skew means fat tail on negative side (may have negative extreme values). If two sides are symmetric, its skew is 0, regardless of how fat the tails are. Gaussian distributions are symmetric so they has zero skew. note that an asymmetric distribution can also has 0 skewness. There is a concept called moments that unify mean, variance, skewness and kurtosis: There is an unbiased way to estimate the thrid central moment μ 3 \mu_3 μ 3 . The deduction of unbiased third central moment estimator is similar to Bessel's correction, but more tricky. A common way of estimating skewness from i.i.d samples, is to use the unbiased third central moment estimator, to divide by cubic of unbiased estimator of standard deviation: But it's still biased, as E [ X Y ] E[\frac{X}{Y}] E [ Y X ] doesn't necessarily equal E [ X ] E [ Y ] \frac{E[X]}{E[Y]} E [ Y ] E [ X ] . Unfortunately, there is no completely unbiased way to estimate skewness from i.i.d samples (unless you have other assumptions about the underlying distribution). The bias gets smaller with more i.i.d samples. Larger kurtosis means it has a fatter tail. The more extreme values it has, the higher its kurtosis. Gaussian distributions have kurtosis of 3. Excess kurtosis is the kurtosis minus 3. A common way of estimating excess kurtosis from i.i.d samples, is to use the unbiased estimator of fourth cumulant ( E [ ( X − E [ X ] ) 4 ] − 3 V a r [ X ] 2 E[(X-E[X])^4]-3Var[X]^2 E [( X − E [ X ] ) 4 ] − 3 Va r [ X ] 2 ), to divide the square of unbiased estimator of variance: It's still biased. If we have some independent samples of X X X , can estimate mean E [ X ] E[X] E [ X ] by calculating average E ^ [ X ] = 1 n ∑ i X i \hat{E}[X]=\frac{1}{n}\sum_i X_i E ^ [ X ] = n 1 ∑ i X i . The variance of calculated average is 1 n Var [ X ] \frac{1}{n} \text{Var}[X] n 1 Var [ X ] , which will reduce by having more samples. However, if the variance of X X X is large and the amount of samples is few, the average will have a large variance, the estimated mean will be inaccurate. We can make the estimation more accurate by using control variate. Then we can estimate E [ X ] E[X] E [ X ] using E ^ [ X + λ ( Y − E [ Y ] ) ] \hat{E}[X+\lambda(Y-E[Y])] E ^ [ X + λ ( Y − E [ Y ])] , where λ \lambda λ is a constant. By choosing the right λ \lambda λ , the estimator can have lower variance than just calculating average of X. The Y here is called a control variate. Some previous knowledge: E [ E ^ [ A ] ] = E [ A ] E[\hat{E}[A]] = E[A] E [ E ^ [ A ]] = E [ A ] , Var [ E ^ [ A ] ] = 1 n Var [ A ] \text{Var}[\hat{E}[A]]=\frac{1}{n}\text{Var}[A] Var [ E ^ [ A ]] = n 1 Var [ A ] . The mean of that estimator is E [ X ] E[X] E [ X ] , meaning that the estimator is unbiased: E [ Y − E [ Y ]] ) = E [ X ] Then calculate the variance of the estimator: − λ E [ Y ] ] = 1 n Var [ X + λ Y ] = 1 n ( Var [ X ] + Var [ λ Y ] + 2 cov [ X , λ Y ] ) = 1 n ( Var [ X ] + λ 2 Var [ Y ] + 2 λ cov [ X , Y ] ) =\frac{1}{n}\text{Var}[X+\lambda Y] = \frac{1}{n}(\text{Var}[X]+\text{Var}[\lambda Y] +2\text{cov}[X,\lambda Y]) = \frac{1}{n}(\text{Var}[X]+\lambda^2 \text{Var}[Y]+2\lambda \text{cov}[X,Y]) = n 1 Var [ X + λY ] = n 1 ( Var [ X ] + Var [ λY ] + 2 cov [ X , λY ]) = n 1 ( Var [ X ] + λ 2 Var [ Y ] + 2 λ cov [ X , Y ]) We want to minimize the variance of estimator by choosing a λ \lambda λ . We want to find a λ \lambda λ that minimizes Var [ Y ] λ 2 + 2 cov [ X , Y ] λ \text{Var}[Y] \lambda^2 + 2\text{cov}[X,Y] \lambda Var [ Y ] λ 2 + 2 cov [ X , Y ] λ . Quadratic funciton knowledge tells a x 2 + b x + c ( a > 0 ) ax^2+bx+c \ \ (a>0) a x 2 + b x + c ( a > 0 ) minimizes when x = − b 2 a x=\frac{-b}{2a} x = 2 a − b , then the optimal lambda is: And by using that optimal λ \lambda λ , the variance of estimator is: If X and Y are correlated, then cov [ X , Y ] 2 Var [ Y ] > 0 \frac{\text{cov}[X,Y]^2}{\text{Var}[Y]} > 0 Var [ Y ] cov [ X , Y ] 2 > 0 , then the new estimator has smaller variance and is more accurate than the simple one. The larger the correlation, the better it can be. Information entropy measures: If we want to measure the amount of information of a specific event, an event E E E 's amount of information as I ( E ) I(E) I ( E ) , there are 3 axioms: Then according to the three axioms, the definition of I I I (self information) is: The base b b b is relative to the unit. We often use the amount of bits as the unit of amount of information. An event with 50% probability has 1 bit of information, then the base will be 2: Then, for a distribution, the expected value of information of one sample is the expected value of I ( E ) I(E) I ( E ) . That defines information entropy H H H : In discrete case: If there exists x x x where P ( x ) = 0 P(x) = 0 P ( x ) = 0 , then it can be ignored in entropy calculation, as lim x → 0 x log x = 0 \lim_{x \to 0} x \log x = 0 lim x → 0 x lo g x = 0 . Information entropy in discrete case is always positive. In continuous case, where f f f is the probability density function, this is called differential entropy: ( X \mathbb{X} X means the set of x x x where f ( x ) ≠ 0 f(x) \neq 0 f ( x ) = 0 , also called support of f f f .) In continuous case the base is often e e e rather than 2. Here log \log lo g by default means log e \log_e lo g e . In discrete case, 0 ≤ P ( x ) ≤ 1 0 \leq P(x) \leq 1 0 ≤ P ( x ) ≤ 1 , log 1 P ( x ) > 0 \log \frac{1}{P(x)} > 0 lo g P ( x ) 1 > 0 , so entropy can never be negative. But in continuous case, probability density function can take value larger than 1, so entropy may be negative. If X and Y are independent, then H ( ( X , Y ) ) = E [ I ( ( X , Y ) ) ] = E [ I ( X ) + I ( Y ) ] = E [ I ( X ) ] + E [ I ( Y ) ] = H ( X ) + H ( Y ) H((X,Y))=E[I((X,Y))]=E[I(X)+I(Y)]=E[I(X)]+E[I(Y)]=H(X)+H(Y) H (( X , Y )) = E [ I (( X , Y ))] = E [ I ( X ) + I ( Y )] = E [ I ( X )] + E [ I ( Y )] = H ( X ) + H ( Y ) . If one fair coin toss has 1 bit entropy, then n independent tosses has n bit entropy. If I split one case into two cases, entropy increases. If I merge two cases into one case, entropy reduces. Because p 1 log 1 p 1 + p 2 log 1 p 2 > ( p 1 + p 2 ) log 1 p 1 + p 2 p_1\log \frac{1}{p_1} + p_2\log \frac{1}{p_2} > (p_1+p_2) \log \frac{1}{p_1+p_2} p 1 lo g p 1 1 + p 2 lo g p 2 1 > ( p 1 + p 2 ) lo g p 1 + p 2 1 (if p 1 ≠ 0 , p 2 ≠ 0 p_1 \neq 0, p_2 \neq 0 p 1 = 0 , p 2 = 0 ), which is because that f ( x ) = log 1 x f(x)=\log \frac{1}{x} f ( x ) = lo g x 1 is convex, so p 1 p 1 + p 2 log 1 p 1 + p 2 p 1 + p 2 log 1 p 2 > log 1 p 1 + p 2 \frac{p_1}{p_1+p_2}\log\frac{1}{p_1}+\frac{p_2}{p_1+p_2}\log\frac{1}{p_2}>\log\frac{1}{p_1+p_2} p 1 + p 2 p 1 lo g p 1 1 + p 1 + p 2 p 2 lo g p 2 1 > lo g p 1 + p 2 1 , then multiply two sides by p 1 + p 2 p_1+p_2 p 1 + p 2 gets the above result. The information entropy is the theorecical minimum of information required to encode a sample. For example, to encode the result of a fair coin toss, we use 1 bit, 0 for head and 1 for tail (reversing is also fine). If the coin is biased to head, to compress the information, we can use 0 for two consecutive heads, 10 for one head, 11 for one tail, which require fewer bits on average for each sample. That may not be optimal, but the most optimal loseless compresion cannot be better than information entropy. In continuous case, if k k k is a positive constant, H ( k X ) = H ( X ) + log k H(kX) = H(X) + \log k H ( k X ) = H ( X ) + lo g k : Entropy is invariant to offset of random variable. H ( X + k ) = H ( X ) H(X+k)=H(X) H ( X + k ) = H ( X ) A joint distribution of X and Y is a distribution where each outcome is a pair of X and Y. Its entropy is called joint information entropy. Here I will use H ( ( X , Y ) ) H((X,Y)) H (( X , Y )) to denote joint entropy (to avoid confusing with cross entropy). If I fix the value of Y as y y y , then see the distribution of X: Take that mean over different Y, we get conditional entropy: Applying conditional probability rule: P ( ( X , Y ) ) = P ( X ∣ Y ) P ( Y ) P((X,Y)) = P(X \vert Y) P(Y) P (( X , Y )) = P ( X ∣ Y ) P ( Y ) So the conditional entropy is defined like this: P ( ( X , Y ) ) = P ( X ∣ Y ) P ( Y ) P((X, Y)) = P(X \vert Y) P(Y) P (( X , Y )) = P ( X ∣ Y ) P ( Y ) . Similarily, H ( ( X , Y ) ) = H ( X ∣ Y ) + H ( Y ) H((X,Y))=H(X \vert Y)+H(Y) H (( X , Y )) = H ( X ∣ Y ) + H ( Y ) . The exact deduction is as follows: If X X X and Y Y Y are not independent, then the joint entropy is smaller than if they are independent: H ( ( X , Y ) ) < H ( X ) + H ( Y ) H((X, Y)) < H(X) + H(Y) H (( X , Y )) < H ( X ) + H ( Y ) . If X and Y are not independent then knowing X will also give some information about Y. This can be deduced by mutual information which will be explained below. Here I A ( x ) I_A(x) I A ( x ) denotes the amount of information of value (event) x x x in distribution A A A . The difference of information of the same value in two distributions A A A and B B B : The KL divergence from A A A to B B B is the expected value of that regarding the probabilities of A A A . Here E A E_A E A means the expected value calculated using A A A 's probabilities: You can think KL divergence as: KL divergence is also called relative entropy. KL divergence is asymmetric, D K L ( A ∥ B ) D_{KL}(A\parallel B) D K L ( A ∥ B ) is different to D K L ( B ∥ A ) D_{KL}(B\parallel A) D K L ( B ∥ A ) . It's often that the first distribution is the real underlying distribution, and the second distribution is an approximation or the model output. If A and B are the same, the KL divergence betweem them are zero. Otherwise, KL divergence is positive. KL divergence can never be negative, will explain later. P B ( x ) P_B(x) P B ( x ) appears on denominator. If there exists x x x that P B ( x ) = 0 P_B(x) = 0 P B ( x ) = 0 and P A ( x ) ≠ 0 P_A(x) \neq 0 P A ( x ) = 0 , then it can be seen that KL divergence is infinite. It can be seen as "The model expect something to never happen but it actually can happen". If there is no such case, we say that A absolutely continuous with respect to B, written as A ≪ B A \ll B A ≪ B . This requires all outcomes from B to include all outcomes from A. Another concept is cross entropy . The cross entropy from A to B, denoted H ( A , B ) H(A, B) H ( A , B ) , is the entropy of A plus KL divergence from A to B: Information entropy H ( X ) H(X) H ( X ) can also be expressed as cross entropy of itself H ( X , X ) H(X, X) H ( X , X ) , similar to the relation between variance and covariance. (In some places H ( A , B ) H(A,B) H ( A , B ) denotes joint entropy. I use H ( ( A , B ) ) H((A,B)) H (( A , B )) for joint entropy to avoid ambiguity.) Cross entropy is also asymmetric. In deep learning cross entropy is often used as loss function. If each piece of training data's distribution's entropy H ( A ) H(A) H ( A ) is fixed, minimizing cross entropy is the same as minimizing KL divergence. Jensen's inequality states that for a concave function f f f : The reverse applies for convex. Here is a visual example showing Jensen's inequality. For example I have a discrete distribution with 5 cases X 1 , X 2 , X 3 , X 4 , X 5 X_1,X_2,X_3,X_4,X_5 X 1 , X 2 , X 3 , X 4 , X 5 (these are possible outcomes of distribution, not samples), corresponding to X coordinates of the red dots. The probabilities of the 5 cases are p 1 , p 2 , p 3 , p 4 , p 5 p_1, p_2, p_3, p_4, p_5 p 1 , p 2 , p 3 , p 4 , p 5 that sum to 1. E [ X ] = p 1 X 1 + p 2 X 2 + p 3 X 3 + p 4 X 4 + p 5 X 5 E[X] = p_1 X_1 + p_2 X_2 + p_3 X_3 + p_4 X_4 + p_5 X_5 E [ X ] = p 1 X 1 + p 2 X 2 + p 3 X 3 + p 4 X 4 + p 5 X 5 . E [ f ( x ) ] = p 1 f ( X 1 ) + p 2 f ( X 2 ) + p 3 f ( X 3 ) + p 4 f ( X 4 ) + p 5 f ( X 5 ) E[f(x)] = p_1 f(X_1) + p_2 f(X_2) + p_3 f(X_3) + p_4 f(X_4) + p_5 f(X_5) E [ f ( x )] = p 1 f ( X 1 ) + p 2 f ( X 2 ) + p 3 f ( X 3 ) + p 4 f ( X 4 ) + p 5 f ( X 5 ) . Then ( E [ X ] , E [ f ( x ) ] ) (E[X], E[f(x)]) ( E [ X ] , E [ f ( x )]) can be seen as an interpolation between five points ( X 1 , f ( X 1 ) ) , ( X 2 , f ( X 2 ) ) , ( X 3 , f ( X 3 ) ) , ( X 4 , f ( X 4 ) ) , ( X 5 , f ( X 5 ) ) (X_1, f(X_1)), (X_2, f(X_2)), (X_3, f(X_3)), (X_4, f(X_4)), (X_5, f(X_5)) ( X 1 , f ( X 1 )) , ( X 2 , f ( X 2 )) , ( X 3 , f ( X 3 )) , ( X 4 , f ( X 4 )) , ( X 5 , f ( X 5 )) , using weights p 1 , p 2 , p 3 , p 4 , p 5 p_1, p_2, p_3, p_4, p_5 p 1 , p 2 , p 3 , p 4 , p 5 . The possible area of the interpolated point correspond to the green convex polygon: For each point in green polygon ( E [ X ] , E [ f ( X ) ] ) (E[X], E[f(X)]) ( E [ X ] , E [ f ( X )]) , the point on function curve with the same X coordinate ( E [ X ] , f ( E [ X ] ) ) (E[X], f(E[X])) ( E [ X ] , f ( E [ X ])) is above it. So E [ f ( X ) ] ≤ f ( E [ X ] ) E[f(X)] \leq f(E[X]) E [ f ( X )] ≤ f ( E [ X ]) . The same applies when you add more cases to the discrete distribution, the convex polygon will have more points but still below the function curve. The same applies to continuous distribution when there are infinitely many cases. Jensen's inequality tells that KL divergence is non-negative: There is a trick that extracting -1 makes P A P_A P A be in denominator that will be cancelled later. The logarithm function is concave. Jensen's inequality gives: Multiplying -1 and flip: The right side equals 0 because: Then how to estimate the KL divergence D K L ( A , B ) D_{KL}(A, B) D K L ( A , B ) ? Reference: Approximating KL Divergence As KL divergence is E A [ log P A ( x ) P B ( x ) ] E_A\left[\log \frac{P_A(x)}{P_B(x)}\right] E A [ lo g P B ( x ) P A ( x ) ] , the simply way is to calculate the average of log P A ( x ) P B ( x ) \log \frac{P_A(x)}{P_B(x)} lo g P B ( x ) P A ( x ) : However it may to be negative in some cases. The true KL divergence can never be negative. This may cause issues. A better way to estimate KL divergence is: ( P A ( x ) = 0 P_A(x) = 0 P A ( x ) = 0 is impossible because it's sampled from A) It's always positive and has no bias. The P B ( x ) P A ( x ) − 1 \frac{P_B(x)}{P_A(x)}-1 P A ( x ) P B ( x ) − 1 is a control variate and is negatively correlated with log P A ( x ) P B ( x ) \log \frac{P_A(x)}{P_B(x)} lo g P B ( x ) P A ( x ) . Recall control variate: if we want to estimate E [ X ] E[X] E [ X ] from samples more accurately, we can find another variable Y Y Y that's correlated with X X X , and we must know its theoretical mean E [ Y ] E[Y] E [ Y ] , then we use E ^ [ X + λ Y ] − λ E [ Y ] \hat E[X+\lambda Y] - \lambda E[Y] E ^ [ X + λY ] − λ E [ Y ] to estimate E [ X ] E[X] E [ X ] . The parameter λ \lambda λ is choosed by minimizing variance. The mean of that control variate is zero, because E x ∼ A [ P B ( x ) P A ( x ) − 1 ] = ∑ x P A ( x ) ( P B ( x ) P A ( x ) − 1 ) = ∑ x ( P B ( x ) − P A ( x ) ) = ∑ x P B ( x ) − ∑ x P A ( x ) = 0 E_{x \sim A}\left[\frac{P_B(x)}{P_A(x)}-1\right]=\sum_x P_A(x) (\frac{P_B(x)}{P_A(x)}-1)=\sum_x (P_B(x) - P_A(x)) =\sum_x P_B(x) - \sum_x P_A(x)=0 E x ∼ A [ P A ( x ) P B ( x ) − 1 ] = ∑ x P A ( x ) ( P A ( x ) P B ( x ) − 1 ) = ∑ x ( P B ( x ) − P A ( x )) = ∑ x P B ( x ) − ∑ x P A ( x ) = 0 The λ = 1 \lambda=1 λ = 1 is not chosen by mimimizing variance, but chosen by making the estimator non-negative. If I define k = P B ( x ) P A ( x ) k=\frac{P_B(x)}{P_A(x)} k = P A ( x ) P B ( x ) , then log P A ( x ) P B ( x ) + λ ( P B ( x ) P A ( x ) − 1 ) = − log k + λ ( k − 1 ) \log \frac{P_A(x)}{P_B(x)} + \lambda(\frac{P_B(x)}{P_A(x)} - 1) = -\log k + \lambda(k-1) lo g P B ( x ) P A ( x ) + λ ( P A ( x ) P B ( x ) − 1 ) = − lo g k + λ ( k − 1 ) . We want it to be non-negative: − log k + λ ( k − 1 ) ≥ 0 -\log k + \lambda(k-1) \geq 0 − lo g k + λ ( k − 1 ) ≥ 0 for all k > 0 k>0 k > 0 , it can be seen that a line y = λ ( k − 1 ) y=\lambda (k-1) y = λ ( k − 1 ) must be above y = log k y=\log k y = lo g k , the only solution is λ = 1 \lambda=1 λ = 1 , where the line is a tangent line on log k \log k lo g k . If X and Y are independent, then H ( ( X , Y ) ) = H ( X ) + H ( Y ) H((X,Y))=H(X)+H(Y) H (( X , Y )) = H ( X ) + H ( Y ) . But if X and Y are not independent, knowing X reduces uncertainty of Y, then H ( ( X , Y ) ) < H ( X ) + H ( Y ) H((X,Y))<H(X)+H(Y) H (( X , Y )) < H ( X ) + H ( Y ) . Mutual information I ( X ; Y ) I(X;Y) I ( X ; Y ) measures how "related" X and Y are: For a joint distribution, if we only care about X, then the distribution of X is a marginal distribution, same as Y. If we treat X and Y as independent, consider a "fake" joint distribution as if X and Y are independent. Denote that "fake" joint distribution as Z Z Z , then P ( Z = ( x , y ) ) = P ( X = x ) P ( Y = y ) P(Z=(x,y))=P(X=x)P(Y=y) P ( Z = ( x , y )) = P ( X = x ) P ( Y = y ) . It's called "outer product of marginal distribution", because its probability matrix is the outer product of two marginal distributions, so it's denoted X ⊗ Y X \otimes Y X ⊗ Y . Then mutual information can be expressed as KL divergence between joint distribution ( X , Y ) (X, Y) ( X , Y ) and that "fake" joint distribution X ⊗ Y X \otimes Y X ⊗ Y : KL divergence is zero when two distributions are the same, and KL divergence is positive when two distributions are not the same. So: H ( ( X , Y ) ) = H ( X ) + H ( Y ) − I ( X ; Y ) H((X,Y))=H(X)+H(Y)-I(X;Y) H (( X , Y )) = H ( X ) + H ( Y ) − I ( X ; Y ) , so if X and Y are not independent then H ( ( X , Y ) ) < H ( X ) + H ( Y ) H((X,Y))<H(X)+H(Y) H (( X , Y )) < H ( X ) + H ( Y ) . Mutual information is symmetric, I ( X ; Y ) = I ( Y ; X ) I(X;Y)=I(Y;X) I ( X ; Y ) = I ( Y ; X ) . As H ( ( X , Y ) ) = H ( X ∣ Y ) + H ( Y ) H((X,Y)) = H(X \vert Y) + H(Y) H (( X , Y )) = H ( X ∣ Y ) + H ( Y ) , so I ( X ; Y ) = H ( X ) + H ( Y ) − H ( ( X , Y ) ) = H ( X ) − H ( X ∣ Y ) I(X;Y) = H(X) + H(Y) - H((X,Y)) = H(X) - H(X \vert Y) I ( X ; Y ) = H ( X ) + H ( Y ) − H (( X , Y )) = H ( X ) − H ( X ∣ Y ) . If knowing Y completely determines X, knowing Y make the distribution of X collapse to one case with 100% probability, then H ( X ∣ Y ) = 0 H(X \vert Y) = 0 H ( X ∣ Y ) = 0 , then I ( X ; Y ) = H ( X ) I(X;Y)=H(X) I ( X ; Y ) = H ( X ) . Some places use correlation factor Cov [ X , Y ] Var [ X ] Var [ Y ] \frac{\text{Cov}[X,Y]}{\sqrt{\text{Var}[X]\text{Var}[Y]}} Var [ X ] Var [ Y ] Cov [ X , Y ] to measure the correlation between two variables. But correlation factor is not accurate in non-linear cases. Mutual information is more accurate in measuring correlation. Information Bottleneck theory tells that the training of neural network will learn an intermediary representation that: If we have two independent random variablex X and Y, and consider the distribution of the sum Z = X + Y Z=X+Y Z = X + Y , then For each z, it sums over different x and y within the constraint z = x + y z=x+y z = x + y . The constraint z = x + y z=x+y z = x + y allows determining y y y from x x x and z z z : y = z − x y=z-x y = z − x , so it can be rewritten as: In continuous case The probability density function of the sum f Z f_Z f Z is denoted as convolution of f X f_X f X and f Y f_Y f Y : The convolution operator ∗ * ∗ can: Convolution can also work in 2D or more dimensions. If X = ( x 1 , x 2 ) X=(x_1,x_2) X = ( x 1 , x 2 ) and Y = ( y 1 , y 2 ) Y=(y_1,y_2) Y = ( y 1 , y 2 ) are 2D random variables (two joint distributions), Z = X + Y = ( z 1 , z 2 ) Z=X+Y=(z_1,z_2) Z = X + Y = ( z 1 , z 2 ) is convolution of X and Y: Convolution can also work on cases where the values are not probabilities. Convolutional neural network uses discrete version of convolution on matrices. Normally when talking about probability we mean the probability of an outcome under a modelled distribution: P ( outcome ∣ modelled distribution ) P(\text{outcome} \ \vert \ \text{modelled distribution}) P ( outcome ∣ modelled distribution ) . But sometimes we have some concrete samples from a distribution but want to know which model suits the best, so we talk about the probability that a model is true given some samples: P ( modelled distribution ∣ outcome ) P(\text{modelled distribution} \ \vert \ \text{outcome}) P ( modelled distribution ∣ outcome ) . If I have some samples, then some parameters make the samples more likely to come from the modelled distribution, and some parameters make the samples less likely to come from the modelled distribution. For example, if I model a coin flip using a parameter θ \theta θ , that and I observe 10 coin flips have 9 heads and 1 tail, then θ = 0.9 \theta=0.9 θ = 0.9 is more likely than θ = 0.5 \theta=0.5 θ = 0.5 . That's straightforward for a simple model. But for more complex models, we need to measure likelihood. Likelihood L ( θ ∣ x 1 , x 2 , . . . , x n ) L(\theta \vert x_1,x_2,...,x_n) L ( θ ∣ x 1 , x 2 , ... , x n ) measures: For example, if I model a coin flip distribution using a parameter θ \theta θ , the probability of head is θ \theta θ and tail is 1 − θ 1-\theta 1 − θ . If I observe 10 coin flip has 9 heads and 1 tail, then the likelihood of θ \theta θ : The more likely a parameter is, the higher its likelihood. If θ \theta θ equals the true underlying parameter then likelihood takes maximum. By taking logarithm, multiply becomes addition, making it easier to analyze. The log-likelihood function: The score function is the derivative of log-likelihood with respect to parameter, for one sample: If θ \theta θ equals true underlying parameter, then mean of likelihood E x [ L ( θ ∣ x ) ] E_x[L(\theta \vert x)] E x [ L ( θ ∣ x )] takes maximum, mean of log-likelihood E x [ log L ( θ ∣ x ) ] E_x[\log L(\theta \vert x)] E x [ lo g L ( θ ∣ x )] also takes maximum. A continuous function's maximum point has zero derivative, so when θ \theta θ is true, then the mean of score function E x [ s ( θ ; x ) ] = ∂ E x [ f ( x ∣ θ ) ] ∂ θ E_x[s(\theta;x)]= \frac{\partial E_x[f(x \vert \theta)]}{\partial \theta} E x [ s ( θ ; x )] = ∂ θ ∂ E x [ f ( x ∣ θ )] is zero. The Fisher information I ( θ ) \mathcal{I}(\theta) I ( θ ) is the mean of the square of score: (The mean is calculated over different outcomes, not different parameters.) We can also think that Fisher information is always computed under the assumption that θ \theta θ is the true underlying parameter, then E x [ s ( θ ; x ) ] = 0 E_x[s(\theta;x)]=0 E x [ s ( θ ; x )] = 0 , then Fisher information is the variance of score I ( θ ) = Var x [ s ( θ ; x ) ] \mathcal{I}(\theta)=\text{Var}_x[s(\theta;x)] I ( θ ) = Var x [ s ( θ ; x )] . Fisher informaiton I ( θ ) \mathcal{I}(\theta) I ( θ ) also measures the curvature of score function, in parameter space, around θ \theta θ . Fisher information measures how much information a sample can tell us about the underlying parameter. When the parameter is an offset and the offset is infinitely small, then the score function is called linear score. If the infinitely small offset is θ \theta θ . The offseted probability density is f 2 ( x ∣ θ ) = f ( x + θ ) f_2(x \vert \theta) = f(x+\theta) f 2 ( x ∣ θ ) = f ( x + θ ) , then ≈ 0 ) = d x d lo g f ( x ) In the places that use score function (and Fisher information) but doesn not specify which parameter, they usually refer to the linear score function. Recall that if we make probability distribution more "spread out" the entropy will increase. If there is no constraint, maximizing entropy of real-number distribution will be "infinitely spread out over all real numbers" (which is not well-defined). But if there are constraints, maximizing entropy will give some common and important distributions: There are other max-entropy distributions. See Wikipedia . We can rediscover these max-entropy distributions, by using Largrange multiplier and functional derivative. To find the distribution with maximum entropy under variance constraint, we can use Largrange multiplier. If we want to find maximum or minimum of f ( x ) f(x) f ( x ) under the constraint that g ( x ) = 0 g(x)=0 g ( x ) = 0 , we can define Largragian function L \mathcal{L} L : Its two partial derivatives have special properties: Then solving equation ∂ L ( x , λ ) ∂ x = 0 \frac{\partial \mathcal{L}(x,\lambda)}{\partial x}=0 ∂ x ∂ L ( x , λ ) = 0 and ∂ L ( x , λ ) ∂ λ = 0 \frac{\partial \mathcal{L}(x,\lambda)}{\partial \lambda}=0 ∂ λ ∂ L ( x , λ ) = 0 will find the maximum or minimum under constraint. Similarily, if there are many constraints, there are multiple λ \lambda λ s. Similar things also apply to functions with multiple arguments. The argument x x x can be a number or even a function, which involves functional derivative: A functional is a function that inputs a function and outputs a value. (One of) its input is a function rather than a value (it's a higher-order function). Functional derivative (also called variational derivative) means the derivative of a functional respect to its argument function. To compute functional derivative, we add a small "perturbation" to the function. f ( x ) f(x) f ( x ) becomes f ( x ) + ϵ ⋅ η ( x ) f(x)+ \epsilon \cdot \eta(x) f ( x ) + ϵ ⋅ η ( x ) , where epsilon ϵ \epsilon ϵ is an infinitely small value that approaches zero, and eta η ( x ) \eta(x) η ( x ) is a test function. The test function can be any function that satisfy some properties. The definition of functional derivative: Note that it's inside integration. For example, this is a functional: G ( f ) = ∫ x f ( x ) d x G(f) = \int x f(x) dx G ( f ) = ∫ x f ( x ) d x . To compute functional derivative ∂ G ( f ) ∂ f \frac{\partial G(f)}{\partial f} ∂ f ∂ G ( f ) , we firstly compute ∂ G ( f + ϵ η ) ∂ ϵ \frac{\partial G(f+\epsilon \eta)}{\partial \epsilon} ∂ ϵ ∂ G ( f + ϵη ) then try to make it into the form of ∫ ∂ G ∂ f ⋅ η ( x ) d x \int \boxed{\frac{\partial G}{\partial f}} \cdot \eta(x) dx ∫ ∂ f ∂ G ⋅ η ( x ) d x Then by pattern matching with the definition, we get ∂ G ∂ f = x \frac{\partial G}{\partial f}=x ∂ f ∂ G = x . Calculate functional derivative for G ( f ) = ∫ x 2 f ( x ) d x G(f)=\int x^2f(x)dx G ( f ) = ∫ x 2 f ( x ) d x : Then ∂ G ∂ f = x 2 \frac{\partial G}{\partial f}=x^2 ∂ f ∂ G = x 2 . Calculate functional derivative for G ( f ) = ∫ ( − f ( x ) log f ( x ) ) d x G(f) = \int (-f(x) \log f(x)) dx G ( f ) = ∫ ( − f ( x ) lo g f ( x )) d x : As log \log lo g is continuous, and ϵ η ( x ) \epsilon \eta(x) ϵη ( x ) is infinitely small, so log ( f ( x ) + ϵ η ( x ) ) = log ( f ( x ) ) \log(f(x)+\epsilon \eta(x))=\log (f(x)) lo g ( f ( x ) + ϵη ( x )) = lo g ( f ( x )) : If we constraint the variance range, a ≤ X ≤ b a \leq X \leq b a ≤ X ≤ b , then maximize its entropy using fuctional derivative We have constraint ∫ a b f ( x ) d x = 1 \int_a^b f(x)dx=1 ∫ a b f ( x ) d x = 1 , which is ∫ a b f ( x ) d x − 1 = 0 \int_a^b f(x)dx-1=0 ∫ a b f ( x ) d x − 1 = 0 . Compute derivatives Solve ∂ L ∂ f = 0 \frac{\partial \mathcal{L}}{\partial f}=0 ∂ f ∂ L = 0 : Solve ∂ L ∂ λ 1 = 0 \frac{\partial \mathcal{L}}{\partial \lambda_1}=0 ∂ λ 1 ∂ L = 0 : The result is f ( x ) = 1 b − a ( a ≤ x ≤ b ) f(x) = \frac 1 {b-a} \ \ \ (a \leq x \leq b) f ( x ) = b − a 1 ( a ≤ x ≤ b ) . The normal distribution, also called Gaussian distribution, is important in statistics. It's the distribution with maximum entropy if we constraint its variance σ 2 \sigma^2 σ 2 to be a finite value. It has two parameters: the mean μ \mu μ and the standard deviation σ \sigma σ . N ( μ , σ 2 ) N(\mu, \sigma^2) N ( μ , σ 2 ) denotes a normal distribution. Changing μ \mu μ moves the PDF alone X axis. Changing σ \sigma σ scales PDF along X axis. We can rediscover normal distribution by maximizing entropy under variance constraint. For a distribution's probability density function f f f , we want to maximize its entropy H ( f ) = ∫ f ( x ) log 1 f ( x ) d x H(f)=\int f(x) \log\frac{1}{f(x)}dx H ( f ) = ∫ f ( x ) lo g f ( x ) 1 d x under the constraint: We can simplify to make deduction easier: The Largragian function: ⎨ ⎧ ∫ − ∞ ∞ f ( x ) lo g f ( x ) 1 d x + λ 1 ( ∫ − ∞ ∞ f ( x ) d x − 1 ) + λ 2 ( ∫ − ∞ ∞ f ( x ) x 2 d x − σ 2 ) = ∫ − ∞ ∞ ( − f ( x ) log f ( x ) + λ 1 f ( x ) + λ 2 x 2 f ( x ) ) d x − λ 1 − λ 2 σ 2 =\int_{-\infty}^{\infty} (-f(x)\log f(x) + \lambda_1 f(x) + \lambda_2 x^2 f(x) ) dx - \lambda_1 - \lambda_2\sigma^2 = ∫ − ∞ ∞ ( − f ( x ) lo g f ( x ) + λ 1 f ( x ) + λ 2 x 2 f ( x )) d x − λ 1 − λ 2 σ 2 Then compute the functional derivative ∂ L ∂ f \frac{\partial \mathcal{L}}{\partial f} ∂ f ∂ L Then solve ∂ L ∂ f = 0 \frac{\partial \mathcal{L}}{\partial f}=0 ∂ f ∂ L = 0 : We get the rough form of normal distribution's probabilify density function. Then solve ∂ L ∂ λ 1 = 0 \frac{\partial \mathcal{L}}{\partial \lambda_1}=0 ∂ λ 1 ∂ L = 0 : That integration must converge, so λ 2 < 0 \lambda_2<0 λ 2 < 0 . A subproblem: solve ∫ − ∞ ∞ e − k x 2 d x \int_{-\infty}^{\infty} e^{-k x^2}dx ∫ − ∞ ∞ e − k x 2 d x ( k > 0 k>0 k > 0 ). The trick is to firstly compute its square ( ∫ − ∞ ∞ e − k x 2 d x ) 2 (\int_{-\infty}^{\infty} e^{-k x^2}dx)^2 ( ∫ − ∞ ∞ e − k x 2 d x ) 2 , turning the integration into two-dimensional, and then substitude polar coordinates x = r cos θ , y = r sin θ , x 2 + y 2 = r 2 , d x d y = r d r d θ x=r \cos \theta, \ y = r \sin \theta, \ x^2+y^2=r^2, \ dx\ dy = r \ dr \ d\theta x = r cos θ , y = r sin θ , x 2 + y 2 = r 2 , d x d y = r d r d θ : Then substitude u = − k r 2 , d u = − 2 k r d r , d r = − 1 2 k r d u u=-kr^2, \ du = -2kr\ dr, \ dr = -\frac{1}{2kr}du u = − k r 2 , d u = − 2 k r d r , d r = − 2 k r 1 d u : So ∫ − ∞ ∞ e − k x 2 d x = π k \int_{-\infty}^{\infty} e^{-kx^2}dx=\sqrt{\frac{\pi}{k}} ∫ − ∞ ∞ e − k x 2 d x = k π . Put − λ 2 = k -\lambda_2=k − λ 2 = k = 1 e − 1 + λ 1 = − λ 2 π e^{-1+\lambda_1} = \sqrt{\frac{-\lambda_2}{\pi}} e − 1 + λ 1 = π − λ 2 Then solve ∂ L ∂ λ 2 = 0 \frac{\partial \mathcal{L}}{\partial \lambda_2}=0 ∂ λ 2 ∂ L = 0 : It requires another trick. For the previous result ∫ − ∞ ∞ e − k x 2 d x = π k \int_{-\infty}^{\infty} e^{-kx^2}dx=\sqrt{\frac{\pi}{k}} ∫ − ∞ ∞ e − k x 2 d x = k π , take derivative to k k k on two sides: k − 2 1 take derivative to k ∫ − ∞ ∞ ( − x 2 ) e ( − x 2 ) k d x = − 2 1 π k − 2 3 So ∫ − ∞ ∞ x 2 e − k x 2 d x = 1 2 π k 3 \int_{-\infty}^{\infty} x^2e^{-kx^2}dx = \frac{1}{2}\sqrt{\frac{\pi}{k^3}} ∫ − ∞ ∞ x 2 e − k x 2 d x = 2 1 k 3 π = σ 2 By using e − 1 + λ 1 = − λ 2 π e^{-1+\lambda_1} = \sqrt{\frac{-\lambda_2}{\pi}} e − 1 + λ 1 = π − λ 2 , we get: ⋅ 2 1 − λ 2 3 π = σ 2 λ 2 2 1 = 2 σ 2 Previously we know that λ 2 < 0 \lambda_2<0 λ 2 < 0 , then λ 2 = − 1 2 σ 2 \lambda_2=-\frac{1}{2\sigma^2} λ 2 = − 2 σ 2 1 . Then e − 1 + λ 1 = 1 2 π σ 2 e^{-1+\lambda_1}=\sqrt{\frac{1}{2\pi\sigma^2}} e − 1 + λ 1 = 2 π σ 2 1 Then we finally deduced the normal distribution's probability density function (when mean is 0): e − 2 σ 2 1 x 2 When mean is not 0, substitute x x x as x − μ x-\mu x − μ , we get the general normal distribution: e − 2 σ 2 1 ( x − μ ) 2 = 2 π σ 1 e − 2 1 ( σ x − μ ) 2 Entropy of normal distribution We can then calculate the entropy of normal distribution: e 2 σ 2 ( x − μ ) 2 ) d x = ∫ f ( x ) ( 1 2 log ( 2 π σ 2 ) + ( x − μ ) 2 2 σ 2 ) d x = 1 2 log ( 2 π σ 2 ) ∫ f ( x ) d x ⏟ = 1 + 1 2 σ 2 ∫ f ( x ) ( x − μ ) 2 ⏟ = σ 2 d x =\int f(x) \left(\frac{1}{2}\log(2\pi\sigma^2)+\frac{(x-\mu)^2}{2\sigma^2}\right)dx=\frac{1}{2}\log(2\pi\sigma^2)\underbrace{\int f(x)dx} _ {=1}+ \frac{1}{2\sigma^2}\underbrace{\int f(x)(x-\mu)^2} _ {=\sigma^2}dx = ∫ f ( x ) ( 2 1 lo g ( 2 π σ 2 ) + 2 σ 2 ( x − μ ) 2 ) d x = 2 1 lo g ( 2 π σ 2 ) = 1 ∫ f ( x ) d x + 2 σ 2 1 = σ 2 ∫ f ( x ) ( x − μ ) 2 d x = 1 2 log ( 2 π σ 2 ) + 1 2 = 1 2 log ( 2 π e σ 2 ) =\frac{1}{2}\log(2\pi\sigma^2)+\frac{1}{2}=\frac{1}{2}\log(2\pi e \sigma^2) = 2 1 lo g ( 2 π σ 2 ) + 2 1 = 2 1 lo g ( 2 π e σ 2 ) If X follows normal distribution and Y's distribution that have the same mean and variance, the cross entropy H ( Y , X ) H(Y,X) H ( Y , X ) have the same value: 1 2 log ( 2 π e σ 2 ) \frac{1}{2}\log(2\pi e \sigma^2) 2 1 lo g ( 2 π e σ 2 ) , regardless of the exact probability density function of Y. The deduction is similar to the above: e 2 σ 2 ( x − μ ) 2 ) d x = ∫ f Y ( x ) ( 1 2 log ( 2 π σ 2 ) + ( x − μ ) 2 2 σ 2 ) d x = 1 2 log ( 2 π σ 2 ) ∫ f Y ( x ) d x ⏟ = 1 + 1 2 σ 2 ∫ f Y ( x ) ( x − μ ) 2 ⏟ = σ 2 d x =\int f_Y(x) \left(\frac{1}{2}\log(2\pi\sigma^2)+\frac{(x-\mu)^2}{2\sigma^2}\right)dx=\frac{1}{2}\log(2\pi\sigma^2)\underbrace{\int f_Y(x)dx} _ {=1}+ \frac{1}{2\sigma^2}\underbrace{\int f_Y(x)(x-\mu)^2} _ {=\sigma^2}dx = ∫ f Y ( x ) ( 2 1 lo g ( 2 π σ 2 ) + 2 σ 2 ( x − μ ) 2 ) d x = 2 1 lo g ( 2 π σ 2 ) = 1 ∫ f Y ( x ) d x + 2 σ 2 1 = σ 2 ∫ f Y ( x ) ( x − μ ) 2 d x = 1 2 log ( 2 π σ 2 ) + 1 2 = 1 2 log ( 2 π e σ 2 ) =\frac{1}{2}\log(2\pi\sigma^2)+\frac{1}{2}=\frac{1}{2}\log(2\pi e \sigma^2) = 2 1 lo g ( 2 π σ 2 ) + 2 1 = 2 1 lo g ( 2 π e σ 2 ) Central limit theorem We have a random variable X X X , which has meam 0 and (finite) variance σ 2 \sigma^2 σ 2 . If we add up n n n independent samples of X X X : X 1 + X 2 + . . . + X n X_1+X_2+...+X_n X 1 + X 2 + ... + X n , the variance of sum is n σ 2 n\sigma^2 n σ 2 . To make its variance constant, we can divide it by n \sqrt n n , then we get S n = X 1 + X 2 + . . . + X n n S_n = \frac{X_1+X_2+...+X_n}{\sqrt n} S n = n X 1 + X 2 + ... + X n . Here S n S_n S n is called the standardized sum , because it makes variance not change by sample count. Central limit theorem says that the standardized sum apporaches normal distribution as n n n increase. No matter what the original distribution of X X X is (as long as its variance is finite), the standardized sum will approach normal distribution. The information of distribution of X X X will be "washed out" during the process. This "washing out information" is also increasing of entropy . As n n n increase, the entopy of standardized sum always increase (except when X follows normal distribution the entropy stays at maximum). H ( S n + 1 ) > H ( S n ) H(S_{n+1}) > H(S_n) H ( S n + 1 ) > H ( S n ) if X X X is not normally distributed. Normal distribution has the maximum entropy under variance constraint. As the entropy of standardized sum increase, its entropy will approach maximum and it will approach normal distribution. This is similar to second law of theomodynamics. This is called Entropic Central Limit Theorem. Proving that is hard and requires a lot of prerequisite knowledges. See also: Solution of Shannon's problem on the monotonicity of entropy , Generalized Entropy Power Inequalities and Monotonicity Properties of Information In the real world, many things follow normal distribution, like height of people, weight of people, error in manufacturing, error in measurement, etc. The height of people is affect by many complex factors (nurtrition, health, genetic factors, exercise, environmental factors, etc.). The combination of these complex factors definitely cannot be similified to a standardized sum of i.i.d zero-mean samples X 1 + X 2 + . . . + X n n \frac{X_1+X_2+...+X_n}{\sqrt n} n X 1 + X 2 + ... + X n . Some factors have large effect and some factors have small effect. The factors are not necessarily independent. But the height of people still roughly follows normal distribution. This can be semi-explained by second law of theomodynamics. The complex interactions of many factors increase entropy of the height. At the same time there are also many factors that constraint the variance of height. Why is there a variance constraint? In some cases variance correspond to instability. A human that is 100 meters tall is impossible as it's physically unstable. Similarily a human that's 1 cm tall is impossible in maintaining normal biological function. The unstable things tend to collapse and vanish (survivorship bias), and the stable things remain. That's how the variance constraint occurs in nature. In some places, variance correspond to energy, and the variance is constrainted by conservation of energy. Although normal distribution is common, not all distributions are normal. There are also many things that follow fat-tail distributions. Also note that Central Limit Theorem works when n n n approaches infinity. Even if a distribution's standardized sum approach normal distribution, the speed of converging is important: some distribution converge to normal quickly, and some slowly. Some fat-tail distribution has finite variance but their standardized sum converge to normal distribution very slowly. In below, bold letter (like x \boldsymbol x x ) means column vector: x 1 x 2 ... x n Linear transform: for a (column) vector x \boldsymbol{x} x , muliply a matrix A A A on it: A x A\boldsymbol x A x is linear transformation. Linear transformation can contain rotation, scaling and shearing. For row vector it's x A \boldsymbol xA x A . Two linear transformations can be combined one, corresponding to matrix multiplication. Affine transform: for a (column) vector x \boldsymbol x x , multiply a matrix on it and then add some offset A x + b A\boldsymbol x + \boldsymbol b A x + b . It can move based on the result of linear transform. Two affine transformations can be combined into one. If y = A x + b , z = C y + d \boldsymbol y=A\boldsymbol x+\boldsymbol b, \boldsymbol z=C\boldsymbol y+\boldsymbol d y = A x + b , z = C y + d , then z = ( C A ) x + ( C b + d ) \boldsymbol z=(CA)\boldsymbol x +(C\boldsymbol b + \boldsymbol d) z = ( C A ) x + ( C b + d ) (in some places affine transformation is called "linear transformation".) Normal distribution has linear properties: Note that the elements of y \boldsymbol y y are no longer necessarily independent. What if I apply two or many affine transformations? Two affine transformations can be combined into one. So the result is still multivariate normal distribution. To describe a multivariate normal distribution, an important concept is covariance matrix . Recall covariance: Cov [ X , Y ] = E [ ( X − E [ X ] ) ( Y − E [ Y ] ) ] \text{Cov}[X,Y]=E[(X-E[X])(Y-E[Y])] Cov [ X , Y ] = E [( X − E [ X ]) ( Y − E [ Y ])] . Some rules about covariance: Covariance matrix: Here E [ x ] E[\boldsymbol x] E [ x ] taking mean of each element in x \boldsymbol x x and output a vector. It's element-wise. E [ x ] i = E [ x i ] E[\boldsymbol x]_i = E[\boldsymbol x_i] E [ x ] i = E [ x i ] . Similar for matrix. The covariance matrix written out: Cov [ x 1 , y 1 ] Cov [ x 2 , y 1 ] ⋮ Cov [ x n , y 1 ] Cov [ x 1 , y 2 ] Cov [ x 2 , y 2 ] ⋮ Cov [ x n , y 2 ] ... ... ⋱ ... Cov [ x 1 , y n ] Cov [ x 2 , y n ] ⋮ Cov [ x n , y n ] Recall that multiplying constant and addition can be "moved out of E [ ] E[] E [ ] ": E [ k X ] = k E [ X ] , E [ X + Y ] = E [ X ] + E [ Y ] E[kX] = k E[X], \ E[X+Y]=E[X]+E[Y] E [ k X ] = k E [ X ] , E [ X + Y ] = E [ X ] + E [ Y ] . If A A A is a matrix that contains random variable and B B B is a matrix that's not random, then E [ A ⋅ B ] = E [ A ] ⋅ B , E [ B ⋅ A ] = B ⋅ E [ A ] E[A\cdot B] = E[A]\cdot B, \ E[B\cdot A] = B\cdot E[A] E [ A ⋅ B ] = E [ A ] ⋅ B , E [ B ⋅ A ] = B ⋅ E [ A ] , because multiplying a matrix come down to multiplying constant and adding up, which all can "move out of E [ ] E[] E [ ] ". Vector can be seen as a special kind of matrix. So applying it to covariance matrix: Similarily, Cov ( x , B ⋅ y ) = Cov ( x , y ) ⋅ B T \text{Cov}(\boldsymbol x, B \cdot \boldsymbol y) = \text{Cov}(\boldsymbol x, \boldsymbol y) \cdot B^T Cov ( x , B ⋅ y ) = Cov ( x , y ) ⋅ B T . If x \boldsymbol x x follows multivariate normal distribution, it can be described by mean vector μ \boldsymbol \mu μ (the mean of each element of x \boldsymbol x x ) and covariance matrix Cov ( x , x ) \text{Cov}(\boldsymbol x,\boldsymbol x) Cov ( x , x ) . Initially, if I have some independent normal variables x 1 , x 2 , . . . x n x_1, x_2, ... x_n x 1 , x 2 , ... x n with mean values μ 1 , . . . , μ n \mu_1, ..., \mu_n μ 1 , ... , μ n and variances σ 1 2 , . . . , σ n 2 \sigma_1^2, ..., \sigma_n^2 σ 1 2 , ... , σ n 2 . If we treat them as a multivariate normal distribution, the mean vector μ x = ( μ 1 , . . . , μ n ) \boldsymbol \mu_x = (\mu_1, ..., \mu_n) μ x = ( μ 1 , ... , μ n ) . The covariance matrix will be diagonal as they are independent: σ 1 2 0 ⋮ 0 0 σ 2 2 ⋮ 0 ... ... ⋱ ... 0 0 ⋮ σ n 2 Then if we apply an affine transformation y = A x + b \boldsymbol y = A \boldsymbol x + \boldsymbol b y = A x + b , then μ y = A μ x + b \boldsymbol \mu_y = A \mu_x + \boldsymbol b μ y = A μ x + b . Cov ( y , y ) = Cov ( A x + b , A x + b ) = Cov ( A x , A x ) = A Cov ( x , x ) A T \text{Cov}(\boldsymbol y,\boldsymbol y) = \text{Cov}(A \boldsymbol x + \boldsymbol b,A \boldsymbol x + \boldsymbol b) = \text{Cov}(A \boldsymbol x, A \boldsymbol x) = A \text{Cov}(\boldsymbol x,\boldsymbol x) A^T Cov ( y , y ) = Cov ( A x + b , A x + b ) = Cov ( A x , A x ) = A Cov ( x , x ) A T . The industry standard of 3D modelling is to model the 3D object as many triangles, called mesh. It only models the visible surface object. It use many triangles to approximate curved surface. Gaussian splatting provides an alternative method of 3D modelling. The 3D scene is modelled by a lot of mutlivariate (3D) gaussian distributions, called gaussian. When rendering, that 3D gaussian distribution is projected onto a plane (screen) and approximately become a 2D gaussian distribution, now probability density correspond to color opacity. Note that the projection is perspective projection (near things big and far things small). Perspective projection is not linear. After perspective projection, the 3D Gaussian distribution is no longer strictly a 2D Gaussian distribution, can be approximated by a 2D Gaussian distribution. Triangle mesh is often modelled by people. But gaussian splatting scene is often trained from photos of different perspectives of a scene. A gaussian's color can be fixed or can change based on different view directions. Gaussian splatting also works in 4D by adding a time dimension. In diffusion model, we add gaussian noise to image (or other things). Then the diffusion model takes noisy input and we train it to output the noise added to it. There will be many steps of adding noise and the model should output the noise added in each step. Tweedie's formula shows that estimating the noise added is the same as computing the likelihood of image distribution. To simplify, here we only consider one dimension and one noise step (the same also applies to many dimensions and many noise steps). If the original value is x 0 x_0 x 0 , we add a noise ϵ ∼ N ( 0 , σ 2 ) \epsilon \sim N(0, \sigma^2) ϵ ∼ N ( 0 , σ 2 ) , the noise-added value is x 1 = x 0 + ϵ x_1 = x_0 + \epsilon x 1 = x 0 + ϵ , x 1 ∼ N ( x 0 , σ 2 ) x_1 \sim N(x_0, \sigma^2) x 1 ∼ N ( x 0 , σ 2 ) . The diffusion model only know x 1 x_1 x 1 and don't know x 0 x_0 x 0 . The diffusion model need to estimate ϵ \epsilon ϵ from x 1 x_1 x 1 . (I use p 1 ∣ 0 ( x 1 ∣ x 0 ) p_{1 \vert 0}(x_1 \vert x_0) p 1∣0 ( x 1 ∣ x 0 ) instead of shorter p ( x 1 ∣ x 0 ) p(x_1 \vert x_0) p ( x 1 ∣ x 0 ) is to reduce confusion between different distributions.) p 1 ∣ 0 ( x 1 ∣ x 0 ) p_{1 \vert 0}(x_1 \vert x_0) p 1∣0 ( x 1 ∣ x 0 ) is a normal distribution: σ 1 e − 2 1 ( σ x 1 − x 0 ) 2 Take log: σ 1 The linear score function under condition: Bayes rule: Take partial derivative to x 1 x_1 x 1 : ∂ x 1 ∂ lo g p 0 ( x 0 ) − ∂ x 1 ∂ lo g p 1 ( x 1 ) Using previous result ∂ log p 1 ∣ 0 ( x 1 ∣ x 0 ) ∂ x 1 = − x 1 − x 0 σ 2 \frac{\partial \log p_{1 \vert 0}(x_1 \vert x_0)}{\partial x_1} = - \frac{x_1-x_0}{\sigma^2} ∂ x 1 ∂ l o g p 1∣0 ( x 1 ∣ x 0 ) = − σ 2 x 1 − x 0 Now if we already know the noise-added value x 1 x_1 x 1 , but we don't know x 0 x_0 x 0 so x 0 x_0 x 0 is uncertain. We want to compute the expectation of x 0 x_0 x 0 under that condition that x 1 x_1 x 1 is known. x 1 ] = x 1 + E x 0 [ σ 2 ∂ log p 0 ∣ 1 ( x 0 ∣ x 1 ) ∂ x 1 ∣ x 1 ] + E x 0 [ σ 2 ∂ log p 1 ( x 1 ) ∂ x 1 ∣ x 1 ] = x_1 + E_{x_0}\left[\sigma^2 \frac{\partial\log p_{0 \vert 1}(x_0 \vert x_1)}{\partial x_1}\biggr\vert x_1\right] + E_{x_0}\left[ \sigma^2\frac{\partial \log p_1(x_1)}{\partial x_1} \biggr\vert x_1\right] = x 1 + E x 0 [ σ 2 ∂ x 1 ∂ lo g p 0∣1 ( x 0 ∣ x 1 ) x 1 ] + E x 0 [ σ 2 ∂ x 1 ∂ lo g p 1 ( x 1 ) x 1 ] Within it, E x 0 [ ∂ log p 0 ∣ 1 ( x 0 ∣ x 1 ) ∂ x 1 ∣ x 1 ] = 0 E_{x_0}\left[ \frac{\partial\log p_{0 \vert 1}(x_0 \vert x_1)}{\partial x_1} \biggr\vert x_1 \right]=0 E x 0 [ ∂ x 1 ∂ l o g p 0∣1 ( x 0 ∣ x 1 ) x 1 ] = 0 , because x 1 ] = ∫ p 0∣1 ( x 0 ∣ x 1 ) ⋅ ∂ x 1 ∂ lo g p 0∣1 ( x 0 ∣ x 1 ) d x 0 = ∫ p 0 ∣ 1 ( x 0 ) ⋅ 1 p 0 ∣ 1 ( x 0 ∣ x 1 ) ⋅ ∂ p 0 ∣ 1 ( x 0 ∣ x 1 ) ∂ x 1 d x 0 = ∫ ∂ p 0 ∣ 1 ( x 0 ∣ x 1 ) ∂ x 1 d x 0 = ∂ ∫ p 0 ∣ 1 ( x 0 ∣ x 1 ) d x 0 ∂ x 1 = ∂ 1 ∂ x 1 = 0 = \int p_{0 \vert 1}(x_0) \cdot \frac 1 {p_{0 \vert 1}(x_0 \vert x_1)} \cdot \frac{\partial p_{0 \vert 1}(x_0 \vert x_1)}{\partial x_1} dx_0 = \int \frac{\partial p_{0 \vert 1}(x_0 \vert x_1)}{\partial x_1} dx_0 = \frac{\partial \int p_{0 \vert 1}(x_0 \vert x_1) dx_0}{\partial x_1} = \frac{\partial 1}{\partial x_1}=0 = ∫ p 0∣1 ( x 0 ) ⋅ p 0∣1 ( x 0 ∣ x 1 ) 1 ⋅ ∂ x 1 ∂ p 0∣1 ( x 0 ∣ x 1 ) d x 0 = ∫ ∂ x 1 ∂ p 0∣1 ( x 0 ∣ x 1 ) d x 0 = ∂ x 1 ∂ ∫ p 0∣1 ( x 0 ∣ x 1 ) d x 0 = ∂ x 1 ∂ 1 = 0 And E x 0 [ σ 2 ∂ log p 1 ( x 1 ) ∂ x 1 ∣ x 1 ] = σ 2 ∂ log p 1 ( x 1 ) ∂ x 1 E_{x_0}\left[ \sigma^2\frac{\partial \log p_1(x_1)}{\partial x_1} \biggr\vert x_1\right] = \sigma^2\frac{\partial \log p_1(x_1)}{\partial x_1} E x 0 [ σ 2 ∂ x 1 ∂ l o g p 1 ( x 1 ) x 1 ] = σ 2 ∂ x 1 ∂ l o g p 1 ( x 1 ) because it's unrelated to random x 0 x_0 x 0 . σ 2 ∂ x 1 ∂ lo g p 1 ( x 1 ) That's Tweedie's formula (for 1D case). It can be generalized to many dimensions, where the x 0 , x 1 x_0, x_1 x 0 , x 1 are vectors, the distributions p 0 , p 1 , p 0 ∣ 1 , p 1 ∣ 0 p_0, p_1, p_{0 \vert 1}, p_{1 \vert 0} p 0 , p 1 , p 0∣1 , p 1∣0 are joint distributions where different dimensions are not necessarily independent. The gaussian noise added to different dimensions are still independent. The diffusion model is trained to estimate the added noise, which is the same as estimating the linear score. If we have constraint X ≥ 0 X \geq 0 X ≥ 0 and fix the mean E [ X ] E[X] E [ X ] to a specific value μ \mu μ , then maximizing entropy gives exponential distribution. It can also be rediscovered from Lagrange multiplier: ⎨ ⎧ ∫ 0 ∞ f ( x ) lo g f ( x ) 1 d x + λ 1 ( ∫ 0 ∞ f ( x ) d x − 1 ) + λ 2 ( ∫ 0 ∞ f ( x ) x d x − μ ) = ∫ 0 ∞ ( − f ( x ) log f ( x ) + λ 1 f ( x ) + λ 2 x f ( x ) ) d x − λ 1 − λ 2 μ =\int_{0}^{\infty} (-f(x)\log f(x) + \lambda_1 f(x) + \lambda_2 x f(x) ) dx - \lambda_1 - \lambda_2\mu = ∫ 0 ∞ ( − f ( x ) lo g f ( x ) + λ 1 f ( x ) + λ 2 x f ( x )) d x − λ 1 − λ 2 μ ∂ L ∂ f = − log f ( x ) − 1 + λ 1 + λ 2 x ∂ L ∂ λ 1 = ∫ 0 ∞ f ( x ) d x − 1 ∂ L ∂ λ 2 = ∫ 0 ∞ x f ( x ) d x − μ \frac{\partial \mathcal{L}}{\partial f} = -\log f(x) - 1 + \lambda_1 + \lambda_2 x \quad\quad\quad \frac{\partial \mathcal{L}}{\partial \lambda_1}=\int_0^{\infty}f(x)dx-1 \quad\quad\quad \frac{\partial \mathcal{L}}{\partial \lambda_2}=\int_0^{\infty} xf(x)dx-\mu ∂ f ∂ L = − lo g f ( x ) − 1 + λ 1 + λ 2 x ∂ λ 1 ∂ L = ∫ 0 ∞ f ( x ) d x − 1 ∂ λ 2 ∂ L = ∫ 0 ∞ x f ( x ) d x − μ Then solve ∂ L ∂ f = 0 \frac{\partial \mathcal{L}}{\partial f}=0 ∂ f ∂ L = 0 : Then solve ∂ L ∂ λ 1 = 0 \frac{\partial \mathcal{L}}{\partial \lambda_1}=0 ∂ λ 1 ∂ L = 0 : To make that integration finite, λ 2 < 0 \lambda_2 < 0 λ 2 < 0 . Let u = λ 2 x , d u = λ 2 d x , d x = 1 λ 2 d u u = \lambda_2 x, \ du = \lambda_2 dx, dx=\frac 1 {\lambda_2} du u = λ 2 x , d u = λ 2 d x , d x = λ 2 1 d u , Then solve ∂ L ∂ λ 2 = 0 \frac{\partial \mathcal{L}}{\partial \lambda_2}=0 ∂ λ 2 ∂ L = 0 : Now we have Solving it gives λ 2 = − 1 μ , e 1 − λ 1 = μ \lambda_2 = - \frac 1 {\mu}, \ e^{1-\lambda_1} = \mu λ 2 = − μ 1 , e 1 − λ 1 = μ . Then In the common definition of exponential distribution, λ = 1 μ \lambda = \frac 1 \mu λ = μ 1 , f ( x ) = λ e − λ x f(x) = \lambda e^{-\lambda x} f ( x ) = λ e − λ x . Its tail function: y = x y = ∞ = e − λ x If some event is happening in fixed rate ( λ \lambda λ ), exponential distribution measures how long do we need to wait for the next event , if how long we will need to wait is irrelevant how long we have aleady waited (memorylessness) . Exponential distribution can measure: How to understand memorlessness? For example, a kind of radioactive atom decays once per 5 minutes on average. If the time unit is minute, then λ = 1 5 \lambda = \frac 1 5 λ = 5 1 . For a specific atom, if we wait for it to decay, the time we need to wait is on average 5 minutes. However, if we have already waited for 3 minutes and it still hasn't decay, the expected time that we need to wait is still 5 minutes. If we have waited for 100 minutes and it still hasn't decay, the expected time that we need to wait is still 5 minutes. Because the atom doesn't "remember" how long we have waited. Memorylessness means the probability that we still need to wait needToWait \text{needToWait} needToWait amount of time is irrelevant to how long we have already waited: (We can also rediscover exponential distrbution from just memorylessness.) Memorylessness is related with its maximum entropy property. Maximizing entropy under constraints means maximizing uncertainty and minizing information other than the constraints. The only two constraints are X ≥ 0 X\geq 0 X ≥ 0 , the wait time is positive, and E [ X ] = 1 λ E[X]=\frac 1 \lambda E [ X ] = λ 1 , the average rate of the event. Other than the two constraints, there is no extra information. No information tells waiting reduces time need to wait, no information tells waiting increases time need to wait. So it's the most unbiased: waiting has no effect on the time need to wait. If the radioactive atom has some "internal memory" that changes over time and controls how likely it will decay, then the waiting time distribution encodes extra information other than the two constraints, which makes it no longer max-entropy. 80/20 rule: for example 80% of weallth are in the richest 20% (the real number may be different). It has fractal property : even within the richest 20%, 80% of wealth are in the richest 20% within. Based on this fractal-like property, we can naturally get Pareto distribution . If the total people count is N N N , the total wealth amount is W W W . Then 0.2 N 0.2N 0.2 N people have 0.8 W 0.8W 0.8 W wealth. Applying the same within the 0.2 N 0.2N 0.2 N people: 0.2 ⋅ 0.2 N 0.2 \cdot 0.2 N 0.2 ⋅ 0.2 N people have 0.8 ⋅ 0.8 W 0.8 \cdot 0.8W 0.8 ⋅ 0.8 W wealth. Applying again, 0.2 ⋅ 0.2 ⋅ 0.2 N 0.2 \cdot 0.2 \cdot 0.2 N 0.2 ⋅ 0.2 ⋅ 0.2 N people have 0.8 ⋅ 0.8 ⋅ 0.8 W 0.8 \cdot 0.8 \cdot 0.8 W 0.8 ⋅ 0.8 ⋅ 0.8 W welath. Generalize it, 0.2 k N 0.2^k N 0. 2 k N people have 0.8 k W 0.8^k W 0. 8 k W wealth ( k k k can be generalized to continuous real number). If the wealth variable is X X X (assume X > 0 X > 0 X > 0 ), its probability density function is f ( x ) f(x) f ( x ) , and porportion of people correspond to probability, the richest 0.2 k 0.2^k 0. 2 k porportion of people group have 0.8 k W 0.8^kW 0. 8 k W wealth, t t t is the wealth threshold (minimum wealth) of that group: Note that f ( x ) f(x) f ( x ) represents probability density function (PDF), which correspond to density of proportion of people. N ⋅ f ( x ) N\cdot f(x) N ⋅ f ( x ) is people amount density over wealth. Multiplying it with wealth x x x and integrate gets total wealth in range: We can rediscover Pareto distribution from these. The first thing to do is extract and eliminate k k k : Then we can take derivative to t t t on two sides: f ( t ) ≠ 0 f(t) \neq 0 f ( t ) = 0 . Divide two sides by − f ( t ) -f(t) − f ( t ) : Take derivative to t t t on two sides again: Now t t t is an argument and can be renamed to x x x . And do some adjustments: Now we get the PDF. We still need to make the total probability area to be 1 to make it a valid distribution. But there is no extra unknown parameter in PDF to change. The solution is to crop the range of X X X . If we set the minimum wealth in distribution to be m m m (but doesn't constraint the maximum wealth), creating constraint X ≥ m X \geq m X ≥ m , then using the previous result Now we rediscovered (a special case of) Pareto distribution from just fractal 80/20 rule. We can generalize it further for other cases like 90/10 rule, 80/10 rule, etc. and get Pareto (Type I) distribution. It has two parameters, shape parameter α \alpha α (correspond to − log 0.2 log 0.8 − log 0.2 = log 5 log 4 ≈ 1.161 -\frac {\log 0.2} {\log 0.8-\log 0.2} = \frac{\log 5}{\log 4} \approx 1.161 − l o g 0.8 − l o g 0.2 l o g 0.2 = l o g 4 l o g 5 ≈ 1.161 ) and minimum value m m m : Note that in real world the wealth of one can be negative (has debts more than assets). The Pareto distribution is just an approximation. m m m means the threshold where Pareto distribution starts to be good approximation. If α ≤ 1 \alpha \leq 1 α ≤ 1 then its theoretical mean is infinite. Of course if we have finite samples then the sample mean will be finite, but if the theoretical mean is infinite, the more sample we have, the larger the sample mean tend to be, and the trend won't stop. If α ≤ 2 \alpha \leq 2 α ≤ 2 then its theoretical variance is infinite. Recall that centrol limit theorem require finite variance. The standarized sum of values taken from Pareto distribution whose α ≤ 2 \alpha \leq 2 α ≤ 2 does not follow central limit theorem because it has infinite variance. Pareto distribution is often described using tail function (rather than probability density function): There are additive values, like length, mass, money. For additive values, we often compute arithmetic average 1 n ( x 1 + x 2 + . . + x n ) \frac 1 n (x_1 + x_2 + .. + x_n) n 1 ( x 1 + x 2 + .. + x n ) . There are also multiplicative values, like asset return rate, growth ratio. For multiplicative values, we often compute geometric average ( x 1 ⋅ x 2 ⋅ . . . ⋅ x n ) 1 n (x_1 \cdot x_2 \cdot ... \cdot x_n)^{\frac 1 n} ( x 1 ⋅ x 2 ⋅ ... ⋅ x n ) n 1 . For example, if an asset grows by 20% in first year, drops 10% in second year and grows 1% in third year, then the average growth ratio per year is ( 1.2 ⋅ 0.9 ⋅ 1.01 ) 1 3 (1.2 \cdot 0.9 \cdot 1.01)^{\frac 1 3} ( 1.2 ⋅ 0.9 ⋅ 1.01 ) 3 1 . Logarithm allows turning multiplication into addition, and turning power into multiplication. If y = log x y = \log x y = lo g x , then log of geometric average of x x x is arithmetic average of y y y : Pareto distribution maximizes entropy under geometric mean constraint E [ log X ] E[\log X] E [ lo g X ] . If we have constraints X ≥ m > 0 X \geq m > 0 X ≥ m > 0 , E [ log X ] = g E[\log X] = g E [ lo g X ] = g , using largrange multiplier to maximize entropy: ⎨ ⎧ ∫ m ∞ f ( x ) lo g f ( x ) 1 d x + λ 1 ( ∫ m ∞ f ( x ) d x − 1 ) + λ 2 ( ∫ m ∞ f ( x ) lo g x d x − g ) L ( f , λ 1 , λ 2 ) = ∫ m ∞ ( − f ( x ) log f ( x ) + λ 1 f ( x ) + λ 2 f ( x ) log x ) d x − λ 1 − g λ 2 \mathcal{L}(f, \lambda_1, \lambda_2) = \int_m^{\infty} (\ -f(x)\log f(x) + \lambda_1 f(x) + \lambda_2 f(x) \log x \ ) dx -\lambda_1 - g \lambda_2 L ( f , λ 1 , λ 2 ) = ∫ m ∞ ( − f ( x ) lo g f ( x ) + λ 1 f ( x ) + λ 2 f ( x ) lo g x ) d x − λ 1 − g λ 2 ∂ L ∂ f = − log f ( x ) − 1 + λ 1 + λ 2 log x \frac{\partial \mathcal{L}}{\partial f} = -\log f(x) - 1 + \lambda_1 + \lambda_2 \log x ∂ f ∂ L = − lo g f ( x ) − 1 + λ 1 + λ 2 lo g x ∂ L ∂ λ 1 = ∫ m ∞ f ( x ) d x − 1 \frac{\partial \mathcal{L}}{\partial \lambda_1} = \int_m^{\infty} f(x) dx -1 ∂ λ 1 ∂ L = ∫ m ∞ f ( x ) d x − 1 ∂ L ∂ λ 2 = ∫ m ∞ f ( x ) log x d x − g \frac{\partial \mathcal{L}}{\partial \lambda_2} = \int_m^{\infty} f(x) \log x \ dx-g ∂ λ 2 ∂ L = ∫ m ∞ f ( x ) lo g x d x − g Solve ∂ L ∂ f = 0 \frac{\partial \mathcal{L}}{\partial f}=0 ∂ f ∂ L = 0 : Solve ∂ L ∂ λ 1 = 0 \frac{\partial \mathcal{L}}{\partial \lambda_1}=0 ∂ λ 1 ∂ L = 0 : To make ∫ m ∞ x λ 2 d x \int_m^{\infty} x^{\lambda_2}dx ∫ m ∞ x λ 2 d x be finite, λ 2 < − 1 \lambda_2 < -1 λ 2 < − 1 . x = m x = ∞ = − λ 2 + 1 1 m λ 2 + 1 = e 1 − λ 1 m λ 2 + 1 λ 2 + 1 = − e 1 − λ 1 e − 1 + λ 1 = − λ 2 + 1 m λ 2 + 1 (1) \frac{m^{\lambda_2+1}}{\lambda_2+1} = -e^{1-\lambda_1} \tag{1}\quad\quad\quad e^{-1+\lambda_1}=-\frac{\lambda_2+1}{m^{\lambda_2+1}} λ 2 + 1 m λ 2 + 1 = − e 1 − λ 1 e − 1 + λ 1 = − m λ 2 + 1 λ 2 + 1 ( 1 ) Solve ∂ L ∂ λ 2 = 0 \frac{\partial \mathcal{L}}{\partial \lambda_2}=0 ∂ λ 2 ∂ L = 0 : If we temporarily ignore e − 1 + λ 1 e^{-1+\lambda_1} e − 1 + λ 1 and compute ∫ m ∞ x λ 2 log x d x \int_m^{\infty} x^{\lambda_2} \log x \ dx ∫ m ∞ x λ 2 lo g x d x . Let u = log x u=\log x u = lo g x , x = e u x=e^u x = e u , d x = e u d u dx = e^udu d x = e u d u : u = l o g m u = ∞ Then By using (1) e − 1 + λ 1 = − λ 2 + 1 m λ 2 + 1 e^{-1+\lambda_1}=-\frac{\lambda_2+1}{m^{\lambda_2+1}} e − 1 + λ 1 = − m λ 2 + 1 λ 2 + 1 : Let α = − 1 log m − g \alpha = -\frac 1 {\log m - g} α = − l o g m − g 1 , it become: Now we rediscovered Pareto (Type I) distribution by maximizing entropy. In the process we have λ 2 < − 1 \lambda_2 \lt -1 λ 2 < − 1 . From λ 2 + 1 = 1 log m − g \lambda_2+1 = \frac 1 {\log m - g} λ 2 + 1 = l o g m − g 1 we know log m − g < 0 \log m - g <0 lo g m − g < 0 , which is m < e g m < e^g m < e g . For example, if wealth follows Pareto distribution, how to compute the wealth share of the top 1%? Generally how to compute the share of the top p p p porpotion? We firstly need to compute the threshold value t t t of the top n n n : Then compute the share x = b x = ∞ = ( − α m − α − α + 1 1 ) b − α + 1 To make that integration finite, we need − α + 1 < 0 -\alpha+1< 0 − α + 1 < 0 , α > 1 \alpha > 1 α > 1 . The share porpotion is irrelevant to m m m . Some concrete numbers: A distribution is power law distribution if its tail function P ( X > x ) P(X>x) P ( X > x ) is roughly porpotional to x − α x^{-\alpha} x − α , where α \alpha α is called exponent. The "roughly" here means that it can have small deviations that is infinitely small when x x x is large enough. Rigorously speaking it's P ( X > x ) ∝ L ( x ) x − α P(X>x) \propto L(x) x^{-\alpha} P ( X > x ) ∝ L ( x ) x − α where L L L is a slow varying function that requires lim x → ∞ L ( r x ) L ( x ) = 1 \lim_{x \to \infty} \frac{L(rx)}{L(x)}=1 lim x → ∞ L ( x ) L ( r x ) = 1 for positive r r r . Note that in some places the power law is written as P ( X > x ) ∝ L ( x ) x − ( α − 1 ) P(X>x) \propto L(x) x^{-(\alpha-1)} P ( X > x ) ∝ L ( x ) x − ( α − 1 ) . In these places the α \alpha α is 1 larger than the α \alpha α in Pareto distribution. The same α \alpha α can have different meaning in different places. Here I will use the α \alpha α that's consistent with the α \alpha α in Pareto distribution. The lower the exponent α \alpha α , the more right-skewed it is, and the more extreme values it have. The power law parameter estimation according to Power laws, Pareto distributions and Zipf’s law : Book The Black Swan also provides some estimation of power law parameter in real world: Note that the estimation is not accurate because they are sensitive to rare extreme samples. Note that there are things whose estimated α < 1 \alpha < 1 α < 1 : intensity of solar flares, intensity of wars, frequency of family names. Recall that in Pareto (Type I) distribuion if α ≤ 1 \alpha \leq 1 α ≤ 1 then the theoretical mean is infinite. The sample mean tend to be higher and higher when we collect samples and the trend won't stop. If the intensity of war do follow power law and the real α < 1 \alpha < 1 α < 1 , then much larger wars exists in the future. Note that most of these things has estimated α < 2 \alpha < 2 α < 2 . In Pareto (Type I) distribution if α ≤ 2 \alpha \leq 2 α ≤ 2 then its theoretical variance is infinite. Not having a finite variance makes them not follow central limit theorem and should not be modelled using gaussian distribution. There are other distributions that can have extreme values: They all have less extreme values than power law distributions, but more extreme values than normal distribution and exponential distribution. If T T T follows exponential distribution, then a T a^T a T follows Pareto (Type I) distribution if a > 1 a>1 a > 1 . If T T T follows exponential distribution, its probability density f T ( t ) = λ e − λ t f_T(t) = \lambda e^{-\lambda t} f T ( t ) = λ e − λ t ( T ≥ 0 T\geq 0 T ≥ 0 ), its cumulative distribution function F T ( t ) = P ( T < t ) = 1 − e − λ t F_T(t) = P(T<t) = 1-e^{-\lambda t} F T ( t ) = P ( T < t ) = 1 − e − λ t If Y = a T Y=a^T Y = a T , a > 1 a>1 a > 1 , then Because T ≥ 0 T\geq 0 T ≥ 0 , Y ≥ a 0 = 1 Y \geq a^0=1 Y ≥ a 0 = 1 . Now Y Y Y 's tail function is in the same form as Pareto (Type I) distribution, where α = λ log a , m = 1 \alpha=\frac{\lambda}{\log a}, \ m =1 α = l o g a λ , m = 1 . If the lifetime of something follows power law distribution, then it has Lindy effect: the longer that it has existed, the longer that it will likely to continue existing. If the lifetime T T T follows Pareto distribution, if something keeps living at time t t t , then compute the expected lifetime under that condition. (The mean is weighted average. The conditional mean is also weighted average but under condition. But as the total integrated weight is not 1, it need to divide the total integrated weight.) (For that integration to be finite, − α + 1 < 0 -\alpha+1<0 − α + 1 < 0 , α > 1 \alpha>1 α > 1 ) The expected lifetime is α α − 1 t \frac{\alpha}{\alpha-1} t α − 1 α t under the condition that it has already lived to time t t t . The expected remaining lifetime is α α − 1 t − t = 1 α − 1 t \frac{\alpha}{\alpha-1} t-t= \frac{1}{\alpha-1}t α − 1 α t − t = α − 1 1 t . It increases by t t t . Lindy effect often doesn't apply to physical things. Lindy effect often applies to information, like technology, culture, art, social norm, etc. If some numbers spans multiple orders of magnitudes, Benford's law says that about 30% of numbers have leading digit 1, about 18% of numbers have leading digit of 2, ... The digit d d d 's porportion is log 10 ( 1 + 1 d ) \log_{10} \left(1 + \frac 1 d \right) lo g 10 ( 1 + d 1 ) . Pareto distribution is a distribution that spans many orders of magnitudes. Let's compute the distribution of first digit if the number follows Pareto distribution. If x x x starts with digit d d d then d 10 k ≤ x < ( d + 1 ) 10 k d 10^k \leq x < (d+1) 10^k d 1 0 k ≤ x < ( d + 1 ) 1 0 k , k = 0 , 1 , 2 , . . . k=0, 1, 2, ... k = 0 , 1 , 2 , ... Pareto distribution has a lower bound m m m . If we make m m m randomly distributed then analytically computing the probability of each starting digit become hard due to edge cases. In this case, doing a Monte Carlo simulation is easier. How to randomly sample numbers from a Pareto distribution? Firstly we know the cumulative distribution function F ( x ) = P ( X < x ) = 1 − P ( X > x ) = 1 − m α x − α F(x) = P(X<x) = 1-P(X>x) = 1- m^\alpha x^{-\alpha} F ( x ) = P ( X < x ) = 1 − P ( X > x ) = 1 − m α x − α . We can then get quantile function, which is the inverse of F F F : F ( x ) = p , Q ( p ) = x F(x)=p, \ \ Q(p) = x F ( x ) = p , Q ( p ) = x Now we can randomly sample p p p between 0 and 1 then Q ( p ) Q(p) Q ( p ) will follow Pareto distribution. Given x x x how to calculate its first digit? If 10 ≤ x < 100 10\leq x<100 10 ≤ x < 100 ( 1 ≤ log 10 x < 2 1 \leq \log_{10} x < 2 1 ≤ lo g 10 x < 2 ) then first digit is ⌊ x 10 ⌋ \lfloor {\frac x {10}} \rfloor ⌊ 10 x ⌋ . If 100 ≤ x < 1000 100 \leq x < 1000 100 ≤ x < 1000 ( 2 ≤ log 10 x < 3 2 \leq \log_{10}x < 3 2 ≤ lo g 10 x < 3 ) then the first digit is ⌊ x 100 ⌋ \lfloor {\frac x {100}} \rfloor ⌊ 100 x ⌋ . Generalize it, the first digit d d d is: Because Pareto distribution has a lot of extreme values, directly calculating the sample will likely to exceed floating-point range and give some . So we need to use log scale. Only calculate using log x \log x lo g x and avoid using x x x directly. Sampling in log scale: Calculating first digit in log scale: When α \alpha α approaches 0 0 0 it accurately follows Benford's law. The larger α \alpha α the larger deviation with Benford's law. If we fix the min value m m m as a specific number, like 3 3 3 , when α \alpha α is not very close to 0 0 0 it significantly deviates with Benford's law. However if we make m m m a random value between 1 and 10 then it will be close to Benford's law. We have a null hypothesis H 0 H_0 H 0 , like "the coin is fair", and an alternative hypothesis H 1 H_1 H 1 , like "the coin is unfair". We now need to test how likely H 1 H_1 H 1 is true using data. If you have some data and it's extreme if we assume null hypothesis H 0 H_0 H 0 , then P-value is the probability of getting the result that's as extreme or more extreme than the data if we assume null hypothesis H 0 H_0 H 0 is true. If p-value is small then the alternative hypothesis is likely true. If I do ten coin flips then get 9 heads and 1 tail, the probability that the coin flip is fair but still get 9 heads and 1 tail. P-value is the probability that we get as extreme or more extreme as the result, and the "extreme" is two sided, so p-value is P ( 9 heads 1 tail ) + P ( 10 heads 0 tail ) + P ( 1 heads 9 tail ) + P ( 0 heads 10 tail ) P(\text{9 heads 1 tail}) + P(\text{10 heads 0 tail}) + P(\text{1 heads 9 tail}) + P(\text{0 heads 10 tail}) P ( 9 heads 1 tail ) + P ( 10 heads 0 tail ) + P ( 1 heads 9 tail ) + P ( 0 heads 10 tail ) assume coin flip is fair. Can we swap the null hypothesis and alternative hypothesis? For two conflicting hypothesis, which one should be the null hypothesis? The key is burden of proof . The null hypothesis is the default that most people tend to agree and does not need proving. The alternative hypothesis is special and require you to prove using the data. The lower the p value, the higher your confidence that alternative hypothesis is true. But due to randomness you cannot be 100% sure. If you are doing an AB test, you keep collecting data, and when there is statistical significance (like p-value lower than 0.05) you make a conclusion, this is not statistically sound. A random fluctation in the process could lead to false positive results. A more rigorous approach is to determine required sample size before AB test. And the fewer data you have the stricter hypothesis test should be (lower p-value threshold). According to O'Brien-Fleming Boundary, the p-value threshold should be 0.001 when you have 25% data, 0.005 when you have 50% data, 0.015 when you have 75% data and 0.045 when you have 100% data. If I have some samples and I calculate values like mean, variance, median, etc. The calculated value is called statistic. The statistics themselves are also random. If you are sure "In 95% probability the real median is between 8.1 and 8.2" then [ 8.1 , 8.2 ] [8.1,8.2] [ 8.1 , 8.2 ] is a confidence interval with 95% confidence level. Confidence interval can measure how uncertain a statistics is. One way of computing confidence interval is called bootstrap . It doesn't require you to assume that the statistic is normally distributed. But it do require the samples to be i.i.d. It works by resample from the data and create many replacements of the data, then calculate the statistics of the replacement data, then get the confidence interval. For example if the original samples are [ 1.0 , 2.0 , 3.0 , 4.0 , 5.0 ] [1.0, 2.0, 3.0, 4.0, 5.0] [ 1.0 , 2.0 , 3.0 , 4.0 , 5.0 ] , resample means randomly select one from original data and repeat 5 times, giving things like [ 4.0 , 2.0 , 4.0 , 5.0 , 2.0 ] [4.0, 2.0, 4.0, 5.0, 2.0] [ 4.0 , 2.0 , 4.0 , 5.0 , 2.0 ] or [ 3.0 , 2.0 , 4.0 , 4.0 , 5.0 ] [3.0, 2.0, 4.0, 4.0, 5.0] [ 3.0 , 2.0 , 4.0 , 4.0 , 5.0 ] (they are likely to contain duplicates). Then compute the statistics for each resample. If the confidence level is 95%, then the confidence interval's lower bound is the 2.5% percentile number in these statistics, and the upper bound is the 97.5% percentile number in these statistics. When we train a model (including deep learning and linear regression) we want it to also work on new data that's not in training set. But the training itself is to change the model parameter to fit training data. Overfitting means the training make the model "memorize" the training data and does not discover the underlying rule in real world that generates training data. Reducing overfitting is a hard topic. The ways to reduce overfitting: Regularization. Force the model to be "simpler". Force the model to compress data. Weight sharing is also regularization (CNN is weight sharing comparing to MLP). Add inductive bias to limit the possibility of model. (The old way of regularization is to simply reduce parameter count, but in deep learning, there is deep double descent effect where more parameter is better.) Make the model more expressive. If the model is not exprssive enough to capture real underlying rule in real world that generates training data, it's simply unable to generalize. An example is that RNN is less expressive than Transformer due to fixed-size state. Make the training data more comprehensive. Reinforcement learning, if done properly, can provide more comprehensive training data than supervised learning, because of the randomness in interacting with environment. How to test how overfit a model is? Frequentist: Probability is an objective thing. We can know probability from the result of repeating a random event many times in the same condition. Bayesian: Probability is a subjective thing. Probability means how you think it's likely to happen based on your initial assumptions and the evidences you see. Probability is relative to the information you have. A discrete distribution can be a table, telling the probability of each possible outcome. A discrete distribuiton can be a function, where the input is a possible outcome and the output is probability. A discrete distribution can be a vector (an array), where i-th number is the probability of i-th outcome. A discrete distribution can be a histogram, where each pillar is a possible outcome, and the height of pillar is probability. A continuous distribution can be described by a probability density function (PDF) f f f . A continuous distribution has infinitely many outcomes, and the probability of each specific outcome is zero (usually). We care about the probability of a range: P ( a < X < b ) = ∫ a b f ( x ) d x P(a<X<b)=\int_a^b f(x)dx P ( a < X < b ) = ∫ a b f ( x ) d x . The integral of the whole range should be 1: ∫ − ∞ ∞ f ( x ) d x = 1 \int_{-\infty}^{\infty}f(x)dx=1 ∫ − ∞ ∞ f ( x ) d x = 1 . The value of PDF can be larger than 1. A distribution can be described by cumulative distribution function. F ( x ) = P ( X ≤ x ) F(x) = P(X \leq x) F ( x ) = P ( X ≤ x ) . It can be integration of PDF: F ( x ) = ∫ − ∞ x f ( x ) d x F(x) = \int_{-\infty}^x f(x)dx F ( x ) = ∫ − ∞ x f ( x ) d x . It start from 0 and monotonically increase then reach 1. Quantile function Q Q Q is the inverse of cumulative distribution function. Q ( p ) = x Q(p) = x Q ( p ) = x means F ( x ) = p F(x)=p F ( x ) = p and P ( X ≤ x ) = p P(X \leq x) = p P ( X ≤ x ) = p . The top 25% value is Q ( 0.75 ) Q(0.75) Q ( 0.75 ) . The bottom 25% value is Q ( 0.25 ) Q(0.25) Q ( 0.25 ) . Prior means what I assume the distribution is before knowing some new information. If I see some new information and improved my understanding of the distribution, then the new distribution that I assume is posterior . The mean of two random variables can add up E [ X + Y ] = E [ X ] + E [ Y ] E [ ∑ i X i ] = ∑ i E [ X i ] E[X + Y] = E[X] + E[Y]\quad \quad \quad E[\sum_iX_i] = \sum_iE[X_i] E [ X + Y ] = E [ X ] + E [ Y ] E [ ∑ i X i ] = ∑ i E [ X i ] Multiplying a random variable by a constant k k k multiplies its mean E [ k X ] = k ⋅ E [ X ] E[kX] = k \cdot E[X] E [ k X ] = k ⋅ E [ X ] A constant's mean is that constant E [ k ] = k E[k] = k E [ k ] = k The theoretical mean is weighted average using theoretical probabilities The estimated mean (empirical mean, sample mean) is non-weighted average over samples The theoretical mean is an accurate value, determined by the theoretical distribution The estimated mean is an inaccurate random variable, because it's calculated from random samples Layer normalization : it works on a vector. It treats each element in a vector as different samples from the same distribution, and then replace each element with their Z-score (using sample mean and sample stdev). Batch normalization : it works on a batch of vectors. It treats the elements in the same index in different vectors in batch as different samples from the same distribtion, and then compute Z-score (using sample mean and sample stdev). The input x = ( x 1 , x 2 , . . . , x n ) \boldsymbol{x} = (x_1,x_2,...,x_n) x = ( x 1 , x 2 , ... , x n ) The vector of ones: 1 = ( 1 , 1 , . . . , 1 ) \boldsymbol{1} = (1, 1, ..., 1) 1 = ( 1 , 1 , ... , 1 ) Computing sample mean can be seen as scaling 1 n \frac 1 n n 1 then dot product with the vector of ones: μ ^ = 1 n x ⋅ 1 {\hat \mu}= \frac 1 n \boldsymbol{x} \cdot \boldsymbol{1} μ ^ = n 1 x ⋅ 1 Subtracting the sample mean can be seen as subtracting μ ^ ⋅ 1 \hat {\mu} \cdot \boldsymbol{1} μ ^ ⋅ 1 , let's call it y \boldsymbol y y : y = x − μ ^ ⋅ 1 = x − 1 n ( x ⋅ 1 ) ⋅ 1 \boldsymbol y = \boldsymbol x - {\hat \mu} \cdot \boldsymbol{1} = \boldsymbol x- \frac 1 n (\boldsymbol{x} \cdot \boldsymbol{1}) \cdot \boldsymbol{1} y = x − μ ^ ⋅ 1 = x − n 1 ( x ⋅ 1 ) ⋅ 1 Recall projection: projecting vector a \boldsymbol a a onto b \boldsymbol b b is ( a ⋅ b b ⋅ b ) ⋅ b (\frac{\boldsymbol a \cdot \boldsymbol b}{\boldsymbol b \cdot \boldsymbol b}) \cdot \boldsymbol b ( b ⋅ b a ⋅ b ) ⋅ b . ( 1 ) 2 = n (\boldsymbol 1)^2 = n ( 1 ) 2 = n . So 1 n ( x ⋅ 1 ) ⋅ 1 \frac 1 n (\boldsymbol{x} \cdot \boldsymbol{1}) \cdot \boldsymbol{1} n 1 ( x ⋅ 1 ) ⋅ 1 is the projection of x \boldsymbol x x onto 1 \boldsymbol 1 1 . Subtracting it means removing the component in the direction of 1 \boldsymbol 1 1 from x \boldsymbol x x . So y \boldsymbol y y is orthogonal to 1 \boldsymbol 1 1 . y \boldsymbol y y is in a hyper-plane orthogonal to 1 \boldsymbol 1 1 . Standard deviation can be seen as the length of y \boldsymbol y y divide by n \sqrt{n} n (or n − 1 \sqrt{n-1} n − 1 ): σ 2 = 1 n ( y ) 2 \boldsymbol\sigma^2 = \frac 1 n (\boldsymbol y)^2 σ 2 = n 1 ( y ) 2 , σ = 1 n ∣ y ∣ \boldsymbol\sigma = \frac 1 {\sqrt{n}} \vert \boldsymbol y \vert σ = n 1 ∣ y ∣ . Dividing by standard deviation can be seen as projecting it onto unit sphere then multiply by n \sqrt n n (or n − 1 \sqrt{n-1} n − 1 ). So computing Z-score can be seen as firstly projecting onto a hyper-plane that's orthogonal to 1 \boldsymbol 1 1 and then projecting onto unit sphere then multiply by n \sqrt n n (or n − 1 \sqrt{n-1} n − 1 ). The n-th moment: E [ X n ] E[X^n] E [ X n ] . Mean is the first moment. The n-th central moment: E [ ( X − μ ) n ] E[(X-\mu)^n] E [( X − μ ) n ] . Variance is the second central moment. The n-th central standardized moment: E [ ( X − μ σ ) n ] E[(\frac{X-\mu}{\sigma})^n] E [( σ X − μ ) n ] . Skewness is the third central standardized moment. Kurtosis is the fourth central standardized moment. we have a random variable Y that's correlated with X we know the true mean of Y: E [ Y ] E[Y] E [ Y ] , How uncertain a distribution is. How much information a sample in that distribution carries. If that event always happens, then it carries zero information. I ( E ) = 0 I(E) = 0 I ( E ) = 0 if P ( E ) = 1 P(E) = 1 P ( E ) = 1 . The more rare an event is, the larger information (more surprise) it carries. I ( E ) I(E) I ( E ) increases as P ( E ) P(E) P ( E ) decreases. The information of two independent events happen together is the sum of the information of each event. Here I use ( X , Y ) (X, Y) ( X , Y ) to denote the combination of X X X and Y Y Y . That means I ( ( X , Y ) ) = I ( X ) + I ( Y ) I((X, Y)) = I(X) + I(Y) I (( X , Y )) = I ( X ) + I ( Y ) if P ( ( X , Y ) ) = P ( X ) ⋅ P ( Y ) P((X, Y)) = P(X) \cdot P(Y) P (( X , Y )) = P ( X ) ⋅ P ( Y ) . This implies the usage of logarithm. A fair coin toss with two cases has 1 bit of information entropy: 0.5 ⋅ log 2 ( 1 0.5 ) + 0.5 ⋅ log 2 ( 1 0.5 ) = 1 0.5 \cdot \log_2(\frac{1}{0.5}) + 0.5 \cdot \log_2(\frac{1}{0.5}) = 1 0.5 ⋅ lo g 2 ( 0.5 1 ) + 0.5 ⋅ lo g 2 ( 0.5 1 ) = 1 bit. If the coin is biased, for example the head has 90% probability and tail 10%, then its entropy is: 0.9 ⋅ log 2 ( 1 0.9 ) + 0.1 ⋅ log 2 ( 1 0.1 ) ≈ 0.47 0.9 \cdot \log_2(\frac{1}{0.9}) + 0.1 \cdot \log_2(\frac{1}{0.1}) \approx 0.47 0.9 ⋅ lo g 2 ( 0.9 1 ) + 0.1 ⋅ lo g 2 ( 0.1 1 ) ≈ 0.47 bits. If it's even more biased, having 99.99% probability of head and 0.01% probability of tail, then its entropy is: 0.9999 ⋅ log 2 ( 1 0.9999 ) + 0.0001 ⋅ log 2 ( 1 0.0001 ) ≈ 0.0015 0.9999 \cdot \log_2(\frac{1}{0.9999}) + 0.0001 \cdot \log_2(\frac{1}{0.0001}) \approx 0.0015 0.9999 ⋅ lo g 2 ( 0.9999 1 ) + 0.0001 ⋅ lo g 2 ( 0.0001 1 ) ≈ 0.0015 bits. If a coin toss is fair but has 0.01% percent of standing up on the table, having 3 cases each with probability 0.0001, 0.49995, 0.49995, then its entropy is 0.0001 ⋅ log 2 ( 1 0.0001 ) + 0.49995 ⋅ log 2 ( 1 0.49995 ) + 0.49995 ⋅ log 2 ( 1 0.49995 ) ≈ 1.0014 0.0001 \cdot \log_2(\frac{1}{0.0001}) + 0.49995 \cdot \log_2(\frac{1}{0.49995}) + 0.49995 \cdot \log_2(\frac{1}{0.49995}) \approx 1.0014 0.0001 ⋅ lo g 2 ( 0.0001 1 ) + 0.49995 ⋅ lo g 2 ( 0.49995 1 ) + 0.49995 ⋅ lo g 2 ( 0.49995 1 ) ≈ 1.0014 bits. (The standing up event itself has about 13.3 bits of information, but its probability is low so it contributed small in information entropy) The "distance" between two distributions. If I "expect" the distribution is B, but the distribution is actually A, how much "surprise" do I get on average. If I design a loseless compression algorithm optimized for B, but use it to compress data from A, then the compression will be not optimal and contain redundant information. KL divergence measures how much redundant information it has on average. We have two distributions: A A A is the target distribution, B B B is the output of our model We have n n n samples from A A A : x 1 , x 2 , . . . x n x_1, x_2, ... x_n x 1 , x 2 , ... x n We know the probablity of each sample in each distribution. We know P A ( x i ) P_A(x_i) P A ( x i ) and P B ( x i ) P_B(x_i) P B ( x i ) Mutual information I ( X ; Y ) I(X;Y) I ( X ; Y ) is zero if the joint distribution ( X , Y ) (X,Y) ( X , Y ) is the same as X ⊗ Y X\otimes Y X ⊗ Y , which means X and Y are independent. Mutual information I ( X ; Y ) I(X;Y) I ( X ; Y ) is positive if X and Y are not independent. Mutual information is never negative, because KL divergence is never negative. Minimize I ( Input ; IntermediaryRepresentation ) I(\text{Input} \ ; \ \text{IntermediaryRepresentation}) I ( Input ; IntermediaryRepresentation ) . Try to compress the intermediary representation and reduce unnecessary information related to input. Maximize I ( IntermediaryRepresentation ; Output ) I(\text{IntermediaryRepresentation} \ ; \ \text{Output}) I ( IntermediaryRepresentation ; Output ) . Try to keep the information in intermediary representation that's releveant to the output as much as possible. In continuous case, convolution takes two probability density functions, and give a new probability density function. In discrete case, convolution can take two functions and give a new function. Each function inputs an outcome and outputs the probability of that outcome. In discrete case, convolution can take two vectors and give a new vector. Each vector's i-th element correspond to the probability of i-th outcome. How likely that we get samples x 1 , x 2 , . . . , x n x_1, x_2, ... , x_n x 1 , x 2 , ... , x n from the modelled distribution using parameter θ \theta θ . how likely a parameter θ \theta θ is the real underlying parameter, given some independent samples x 1 , x 2 , . . . , x n x_1,x_2,...,x_n x 1 , x 2 , ... , x n . If I assume that the coin flip is fair, θ = 0.5 \theta=0.5 θ = 0.5 , then likelihood is about 0.000977. If I assume θ = 0.9 \theta=0.9 θ = 0.9 , then likelihood is about 0.387, which is larger. If I assume θ = 0.999 \theta=0.999 θ = 0.999 then likelihood is about 0.00099, which is smaller than when assuming θ = 0.9 \theta=0.9 θ = 0.9 . It's a valid probability density function: ∫ − ∞ ∞ f ( x ) d x = 1 \int_{-\infty}^{\infty} f(x)dx=1 ∫ − ∞ ∞ f ( x ) d x = 1 , and f ( x ) ≥ 0 f(x) \geq 0 f ( x ) ≥ 0 The mean: ∫ − ∞ ∞ x f ( x ) d x = μ \int_{-\infty}^{\infty} x f(x) dx = \mu ∫ − ∞ ∞ x f ( x ) d x = μ The variance constraint: ∫ − ∞ ∞ f ( x ) ( x − μ ) 2 d x = σ 2 \int_{-\infty}^{\infty} f(x) (x-\mu)^2 dx = \sigma^2 ∫ − ∞ ∞ f ( x ) ( x − μ ) 2 d x = σ 2 Moving the probability density function along X axis doesn't change entropy, so we can fix the mean as 0 (we can replace x x x as x − μ x-\mu x − μ after finishing deduction). log 1 f ( x ) \log\frac{1}{f(x)} lo g f ( x ) 1 already implicitly tells f ( x ) > 0 f(x)>0 f ( x ) > 0 It turns out that the mean constraint ∫ − ∞ ∞ x f ( x ) d x = 0 \int_{-\infty}^{\infty} x f(x) dx = 0 ∫ − ∞ ∞ x f ( x ) d x = 0 is not necessary to deduce the result, so we can not include it in Largrange multipliers. (Including it is also fine but will make it more complex.) if you multiply a constant, the result still follow normal distribution. X ∼ N → k X ∼ N X \sim N \rightarrow \ kX \sim N X ∼ N → k X ∼ N if you add a constant, the result still follow normal distribution. X ∼ N → ( X + k ) ∼ N X \sim N \rightarrow (X+k) \sim N X ∼ N → ( X + k ) ∼ N If you add up two independent normal random variables, the result still follows normal distribution. X ∼ N , Y ∼ N → ( X + Y ) ∼ N X \sim N, Y \sim N \rightarrow (X+Y) \sim N X ∼ N , Y ∼ N → ( X + Y ) ∼ N A linear combination of many independent normal distributions also follow normal distribution. X 1 ∼ N , X 2 ∼ N , . . . X n ∼ N → ( k 1 X 1 + k 2 X 2 + . . . + k n X n ) ∼ N X_1 \sim N, X_2 \sim N, ... X_n \sim N \rightarrow (k_1X_1 + k_2X_2 + ... + k_nX_n) \sim N X 1 ∼ N , X 2 ∼ N , ... X n ∼ N → ( k 1 X 1 + k 2 X 2 + ... + k n X n ) ∼ N We have a (row) vector x \boldsymbol x x of independent random variables x = ( x 1 , x 2 , . . . x n ) \boldsymbol x=(x_1, x_2, ... x_n) x = ( x 1 , x 2 , ... x n ) , each element in vector follows a normal distribution (not necessarily the same normal distribution), then, if we apply an affine transformation on that vector, which means multipling a matrix A A A and then adding an offset b \boldsymbol b b , y = A x + b \boldsymbol y=A\boldsymbol x+\boldsymbol b y = A x + b , then each element of y \boldsymbol y y is a linear combination of normal distributions, y i = x 1 A i , 1 + x 2 A i , 2 + . . . x n A i , n + b i y_i=x_1 A_{i,1} + x_2 A_{i, 2} + ... x_n A_{i,n} + b_i y i = x 1 A i , 1 + x 2 A i , 2 + ... x n A i , n + b i , so each element in y \boldsymbol y y also follow normal distribution. Now y \boldsymbol y y follows multivariate normal distribution. It's symmetric: Cov [ X , Y ] = Cov [ Y , X ] \text{Cov}[X,Y] = \text{Cov}[Y,X] Cov [ X , Y ] = Cov [ Y , X ] If X and Y are independent, Cov [ X , Y ] = 0 \text{Cov}[X,Y]=0 Cov [ X , Y ] = 0 Adding constant Cov [ X + k , Y ] = Cov [ X , Y ] \text{Cov}[X+k,Y] = \text{Cov}[X,Y] Cov [ X + k , Y ] = Cov [ X , Y ] . Variance is invariant to translation. Multiplying constant Cov [ k ⋅ X , Y ] = k ⋅ Cov [ X , Y ] \text{Cov}[k\cdot X,Y] = k \cdot \text{Cov}[X,Y] Cov [ k ⋅ X , Y ] = k ⋅ Cov [ X , Y ] Addition Cov [ X + Y , Z ] = Cov [ X , Z ] + Cov [ Y , Z ] \text{Cov}[X+Y,Z] = \text{Cov}[X,Z]+\text{Cov}[Y,Z] Cov [ X + Y , Z ] = Cov [ X , Z ] + Cov [ Y , Z ] p 0 ( x 0 ) p_0(x_0) p 0 ( x 0 ) is the probability density of original clean value (for image generation, it correspond to the probability distribution of images that we want to generate) p 1 ( x 1 ) p_1(x_1) p 1 ( x 1 ) is the probability density of noise-added value p 1 ∣ 0 ( x 1 ∣ x 0 ) p_{1 \vert 0}(x_1 \vert x_0) p 1∣0 ( x 1 ∣ x 0 ) is the probability density of noise-added value, given clean training data x 0 x_0 x 0 . It's a normal distribution given x 0 x_0 x 0 . It can also be seen as a function that take two arguments x 0 , x 1 x_0, x_1 x 0 , x 1 . p 0 ∣ 1 ( x 0 ∣ x 1 ) p_{0 \vert 1}(x_0 \vert x_1) p 0∣1 ( x 0 ∣ x 1 ) is the probability density of the original clean value given noise-added value. It can also be seen as a function that take two arguments x 0 , x 1 x_0, x_1 x 0 , x 1 . The lifetime of machine components. The time until a radioactive atom decays. The time length of phone calls. The time interval between two packets for a router. Log-normal distribution : If log X \log X lo g X is normally distributed, then X X X follows log-normal distribution. Put in another way, if Y Y Y is normally distributed, then e Y e^Y e Y follows log-normal distribution. Stretched exponential distribution : P ( X > x ) P(X>x) P ( X > x ) is roughly porpotional to e − k x β e^{-kx^\beta} e − k x β ( β < 1 \beta < 1 β < 1 ) Power law with exponential cutoff : P ( X > x ) P(X>x) P ( X > x ) is roughly porpotional to x − α e − λ x x^{-\alpha} e^{-\lambda x} x − α e − λ x Regularization. Force the model to be "simpler". Force the model to compress data. Weight sharing is also regularization (CNN is weight sharing comparing to MLP). Add inductive bias to limit the possibility of model. (The old way of regularization is to simply reduce parameter count, but in deep learning, there is deep double descent effect where more parameter is better.) Make the model more expressive. If the model is not exprssive enough to capture real underlying rule in real world that generates training data, it's simply unable to generalize. An example is that RNN is less expressive than Transformer due to fixed-size state. Make the training data more comprehensive. Reinforcement learning, if done properly, can provide more comprehensive training data than supervised learning, because of the randomness in interacting with environment. Separate the data into training set and test set. Only train using training set and check model performance on test set. Test sensitivity to random fluctation. We can add randomness to parameter, input, hyperparameter, etc., then see model performance. An overfit model is more prone to random perturbation because memorization is more "fragile" than real underlying rule. Survivorship bias and selection bias. Simpson's paradox and base rate fallacy. Confusing correlation with causalty. Try too many different hypothesis. Spurious correlations Collect data until significance. Wrongly remove outliers.
Programming
C#
C++
Go
JavaScript
Design
Tutorial
Rust
Frontend
Web Development
CSS
Sql
Java
HTML
API
Backend
Data Analysis
Statistics