Classification, Statistical Learning within R

This text contains:

• Logistic regression
• Linear Discriminant Analysis
• K-Nearest Neighbors

require(ISLR)
names(Smarket)

## [1] "Year"      "Lag1"      "Lag2"      "Lag3"      "Lag4"      "Lag5"     
## [7] "Volume"    "Today"     "Direction"

summary(Smarket)

##       Year           Lag1             Lag2             Lag3       
##  Min.   :2001   Min.   :-4.922   Min.   :-4.922   Min.   :-4.922  
##  1st Qu.:2002   1st Qu.:-0.640   1st Qu.:-0.640   1st Qu.:-0.640  
##  Median :2003   Median : 0.039   Median : 0.039   Median : 0.038  
##  Mean   :2003   Mean   : 0.004   Mean   : 0.004   Mean   : 0.002  
##  3rd Qu.:2004   3rd Qu.: 0.597   3rd Qu.: 0.597   3rd Qu.: 0.597  
##  Max.   :2005   Max.   : 5.733   Max.   : 5.733   Max.   : 5.733  
##       Lag4             Lag5            Volume          Today       
##  Min.   :-4.922   Min.   :-4.922   Min.   :0.356   Min.   :-4.922  
##  1st Qu.:-0.640   1st Qu.:-0.640   1st Qu.:1.257   1st Qu.:-0.640  
##  Median : 0.038   Median : 0.038   Median :1.423   Median : 0.038  
##  Mean   : 0.002   Mean   : 0.006   Mean   :1.478   Mean   : 0.003  
##  3rd Qu.: 0.597   3rd Qu.: 0.597   3rd Qu.:1.642   3rd Qu.: 0.597  
##  Max.   : 5.733   Max.   : 5.733   Max.   :3.152   Max.   : 5.733  
##  Direction 
##  Down:602  
##  Up  :648  

# ?Smarket
pairs(Smarket, col = Smarket$Direction)

Logistic regression

glm.fit = glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume, data = Smarket, 
    family = binomial)
summary(glm.fit)

## 
## Call:
## glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + 
##     Volume, family = binomial, data = Smarket)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
##  -1.45   -1.20    1.07    1.15    1.33  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.12600    0.24074   -0.52     0.60
## Lag1        -0.07307    0.05017   -1.46     0.15
## Lag2        -0.04230    0.05009   -0.84     0.40
## Lag3         0.01109    0.04994    0.22     0.82
## Lag4         0.00936    0.04997    0.19     0.85
## Lag5         0.01031    0.04951    0.21     0.83
## Volume       0.13544    0.15836    0.86     0.39
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1731.2  on 1249  degrees of freedom
## Residual deviance: 1727.6  on 1243  degrees of freedom
## AIC: 1742
## 
## Number of Fisher Scoring iterations: 3

glm.probs = predict(glm.fit, type = "response")
glm.probs[1:5]

##      1      2      3      4      5 
## 0.5071 0.4815 0.4811 0.5152 0.5108

# use training data ase test
glm.pred = ifelse(glm.probs > 0.5, "Up", "Down")
attach(Smarket)
table(glm.pred, Direction)

##         Direction
## glm.pred Down  Up
##     Down  145 141
##     Up    457 507

mean(glm.pred == Direction)

## [1] 0.5216

# Make training and test set
train = Year < 2005
glm.fit = glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume, data = Smarket, 
    family = binomial, subset = train)
glm.probs = predict(glm.fit, newdata = Smarket[!train, ], type = "response")
glm.pred = ifelse(glm.probs > 0.5, "Up", "Down")
Direction.2005 = Smarket$Direction[!train]
table(glm.pred, Direction.2005)

##         Direction.2005
## glm.pred Down Up
##     Down   77 97
##     Up     34 44

mean(glm.pred == Direction.2005)

## [1] 0.4802

# Fit smaller model
glm.fit = glm(Direction ~ Lag1 + Lag2, data = Smarket, family = binomial, subset = train)
glm.probs = predict(glm.fit, newdata = Smarket[!train, ], type = "response")
glm.pred = ifelse(glm.probs > 0.5, "Up", "Down")
table(glm.pred, Direction.2005)

##         Direction.2005
## glm.pred Down  Up
##     Down   35  35
##     Up     76 106

mean(glm.pred == Direction.2005)

## [1] 0.5595

106/(76 + 106)

## [1] 0.5824

Linear Discriminant Analysis

require(MASS)

lda.fit = lda(Direction ~ Lag1 + Lag2, data = Smarket, subset = Year < 2005)
lda.fit

## Call:
## lda(Direction ~ Lag1 + Lag2, data = Smarket, subset = Year < 
##     2005)
## 
## Prior probabilities of groups:
##  Down    Up 
## 0.492 0.508 
## 
## Group means:
##          Lag1     Lag2
## Down  0.04279  0.03389
## Up   -0.03955 -0.03133
## 
## Coefficients of linear discriminants:
##          LD1
## Lag1 -0.6420
## Lag2 -0.5135

plot(lda.fit)

Smarket.2005 = subset(Smarket, Year == 2005)
lda.pred = predict(lda.fit, Smarket.2005)
class(lda.pred)

## [1] "list"

data.frame(lda.pred)[1:5, ]

##      class posterior.Down posterior.Up      LD1
## 999     Up         0.4902       0.5098  0.08293
## 1000    Up         0.4792       0.5208  0.59114
## 1001    Up         0.4668       0.5332  1.16723
## 1002    Up         0.4740       0.5260  0.83335
## 1003    Up         0.4928       0.5072 -0.03793

table(lda.pred$class, Smarket.2005$Direction)

##       
##        Down  Up
##   Down   35  35
##   Up     76 106

mean(lda.pred$class == Smarket.2005$Direction)

## [1] 0.5595

K-Nearest Neighbors

library(class)
# ?knn
attach(Smarket)

Xlag = cbind(Lag1, Lag2)
train = Year < 2005
knn.pred = knn(Xlag[train, ], Xlag[!train, ], Direction[train], k = 10)
table(knn.pred, Direction[!train])

##         
## knn.pred Down Up
##     Down   51 62
##     Up     60 79

mean(knn.pred == Direction[!train])

## [1] 0.5159

---

Credit

Please note, this material is extracted from online Statistical Learning cource at Stanford University by Prof. T Hastie and Prof R. Tibshirani. It aims only for quick and future references in R and statistical learning. Please visit course page for more information and materials.

---

Linear Regression, Statistical Learning within R

This text contains:

• Simple linear regression
• Multiple linear regression
• Nonlinear terms and Interactions
• Qualitative predictors fix

# Libaries containing example data
library(MASS)
library(ISLR)

Simple linear regression

names(Boston)

##  [1] "crim"    "zn"      "indus"   "chas"    "nox"     "rm"      "age"    
##  [8] "dis"     "rad"     "tax"     "ptratio" "black"   "lstat"   "medv"

## ?Boston
plot(medv ~ lstat, Boston)
fit1 = lm(medv ~ lstat, data = Boston)
fit1

## 
## Call:
## lm(formula = medv ~ lstat, data = Boston)
## 
## Coefficients:
## (Intercept)        lstat  
##       34.55        -0.95

summary(fit1)

## 
## Call:
## lm(formula = medv ~ lstat, data = Boston)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -15.17  -3.99  -1.32   2.03  24.50 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  34.5538     0.5626    61.4   <2e-16 ***
## lstat        -0.9500     0.0387   -24.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.22 on 504 degrees of freedom
## Multiple R-squared:  0.544,  Adjusted R-squared:  0.543 
## F-statistic:  602 on 1 and 504 DF,  p-value: <2e-16

abline(fit1, col = "red")

names(fit1)

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"

confint(fit1)

##              2.5 % 97.5 %
## (Intercept) 33.448 35.659
## lstat       -1.026 -0.874

predict(fit1, data.frame(lstat = c(5, 10, 15)), interval = "confidence")

##     fit   lwr   upr
## 1 29.80 29.01 30.60
## 2 25.05 24.47 25.63
## 3 20.30 19.73 20.87

Multiple linear regression

fit2 = lm(medv ~ lstat + age, data = Boston)
summary(fit2)

## 
## Call:
## lm(formula = medv ~ lstat + age, data = Boston)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -15.98  -3.98  -1.28   1.97  23.16 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  33.2228     0.7308   45.46   <2e-16 ***
## lstat        -1.0321     0.0482  -21.42   <2e-16 ***
## age           0.0345     0.0122    2.83   0.0049 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.17 on 503 degrees of freedom
## Multiple R-squared:  0.551,  Adjusted R-squared:  0.549 
## F-statistic:  309 on 2 and 503 DF,  p-value: <2e-16

fit3 = lm(medv ~ ., Boston)
summary(fit3)

## 
## Call:
## lm(formula = medv ~ ., data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.594  -2.730  -0.518   1.777  26.199 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.65e+01   5.10e+00    7.14  3.3e-12 ***
## crim        -1.08e-01   3.29e-02   -3.29  0.00109 ** 
## zn           4.64e-02   1.37e-02    3.38  0.00078 ***
## indus        2.06e-02   6.15e-02    0.33  0.73829    
## chas         2.69e+00   8.62e-01    3.12  0.00193 ** 
## nox         -1.78e+01   3.82e+00   -4.65  4.2e-06 ***
## rm           3.81e+00   4.18e-01    9.12  < 2e-16 ***
## age          6.92e-04   1.32e-02    0.05  0.95823    
## dis         -1.48e+00   1.99e-01   -7.40  6.0e-13 ***
## rad          3.06e-01   6.63e-02    4.61  5.1e-06 ***
## tax         -1.23e-02   3.76e-03   -3.28  0.00111 ** 
## ptratio     -9.53e-01   1.31e-01   -7.28  1.3e-12 ***
## black        9.31e-03   2.69e-03    3.47  0.00057 ***
## lstat       -5.25e-01   5.07e-02  -10.35  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.75 on 492 degrees of freedom
## Multiple R-squared:  0.741,  Adjusted R-squared:  0.734 
## F-statistic:  108 on 13 and 492 DF,  p-value: <2e-16

par(mfrow = c(2, 2))
plot(fit3)

fit4 = update(fit3, ~. - age - indus)
summary(fit4)

## 
## Call:
## lm(formula = medv ~ crim + zn + chas + nox + rm + dis + rad + 
##     tax + ptratio + black + lstat, data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.598  -2.739  -0.505   1.727  26.237 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  36.34115    5.06749    7.17  2.7e-12 ***
## crim         -0.10841    0.03278   -3.31  0.00101 ** 
## zn            0.04584    0.01352    3.39  0.00075 ***
## chas          2.71872    0.85424    3.18  0.00155 ** 
## nox         -17.37602    3.53524   -4.92  1.2e-06 ***
## rm            3.80158    0.40632    9.36  < 2e-16 ***
## dis          -1.49271    0.18573   -8.04  6.8e-15 ***
## rad           0.29961    0.06340    4.73  3.0e-06 ***
## tax          -0.01178    0.00337   -3.49  0.00052 ***
## ptratio      -0.94652    0.12907   -7.33  9.2e-13 ***
## black         0.00929    0.00267    3.47  0.00056 ***
## lstat        -0.52255    0.04742  -11.02  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.74 on 494 degrees of freedom
## Multiple R-squared:  0.741,  Adjusted R-squared:  0.735 
## F-statistic:  128 on 11 and 494 DF,  p-value: <2e-16

Nonlinear terms and Interactions

# Note: * in formula means interaction, not multiply
fit5 = lm(medv ~ lstat * age, Boston)
summary(fit5)

## 
## Call:
## lm(formula = medv ~ lstat * age, data = Boston)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -15.81  -4.04  -1.33   2.08  27.55 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 36.088536   1.469835   24.55  < 2e-16 ***
## lstat       -1.392117   0.167456   -8.31  8.8e-16 ***
## age         -0.000721   0.019879   -0.04    0.971    
## lstat:age    0.004156   0.001852    2.24    0.025 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.15 on 502 degrees of freedom
## Multiple R-squared:  0.556,  Adjusted R-squared:  0.553 
## F-statistic:  209 on 3 and 502 DF,  p-value: <2e-16

# square of lstat
fit6 = lm(medv ~ lstat + I(lstat^2), Boston)
summary(fit6)

## 
## Call:
## lm(formula = medv ~ lstat + I(lstat^2), data = Boston)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -15.28  -3.83  -0.53   2.31  25.41 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 42.86201    0.87208    49.1   <2e-16 ***
## lstat       -2.33282    0.12380   -18.8   <2e-16 ***
## I(lstat^2)   0.04355    0.00375    11.6   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.52 on 503 degrees of freedom
## Multiple R-squared:  0.641,  Adjusted R-squared:  0.639 
## F-statistic:  449 on 2 and 503 DF,  p-value: <2e-16

attach(Boston)
par(mfrow = c(1, 1))
plot(medv ~ lstat)
points(lstat, fitted(fit6), col = "red", pch = 20)
# or directly polynomical function can be used
fit7 = lm(medv ~ poly(lstat, 4))
points(lstat, fitted(fit7), col = "blue", pch = 20)

Qualitative predictors

# fix(Carseats) # you may use fix command to edit a data
names(Carseats)

##  [1] "Sales"       "CompPrice"   "Income"      "Advertising" "Population" 
##  [6] "Price"       "ShelveLoc"   "Age"         "Education"   "Urban"      
## [11] "US"

summary(Carseats)

##      Sales         CompPrice       Income       Advertising   
##  Min.   : 0.00   Min.   : 77   Min.   : 21.0   Min.   : 0.00  
##  1st Qu.: 5.39   1st Qu.:115   1st Qu.: 42.8   1st Qu.: 0.00  
##  Median : 7.49   Median :125   Median : 69.0   Median : 5.00  
##  Mean   : 7.50   Mean   :125   Mean   : 68.7   Mean   : 6.63  
##  3rd Qu.: 9.32   3rd Qu.:135   3rd Qu.: 91.0   3rd Qu.:12.00  
##  Max.   :16.27   Max.   :175   Max.   :120.0   Max.   :29.00  
##    Population      Price      ShelveLoc        Age         Education   
##  Min.   : 10   Min.   : 24   Bad   : 96   Min.   :25.0   Min.   :10.0  
##  1st Qu.:139   1st Qu.:100   Good  : 85   1st Qu.:39.8   1st Qu.:12.0  
##  Median :272   Median :117   Medium:219   Median :54.5   Median :14.0  
##  Mean   :265   Mean   :116                Mean   :53.3   Mean   :13.9  
##  3rd Qu.:398   3rd Qu.:131                3rd Qu.:66.0   3rd Qu.:16.0  
##  Max.   :509   Max.   :191                Max.   :80.0   Max.   :18.0  
##  Urban       US     
##  No :118   No :142  
##  Yes:282   Yes:258  
##                     
##                     
##                     
## 

# use all fields and interaction of Income:Advertising and Age:Price
fit1 = lm(Sales ~ . + Income:Advertising + Age:Price, Carseats)
summary(fit1)

## 
## Call:
## lm(formula = Sales ~ . + Income:Advertising + Age:Price, data = Carseats)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.921 -0.750  0.018  0.675  3.341 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         6.575565   1.008747    6.52  2.2e-10 ***
## CompPrice           0.092937   0.004118   22.57  < 2e-16 ***
## Income              0.010894   0.002604    4.18  3.6e-05 ***
## Advertising         0.070246   0.022609    3.11  0.00203 ** 
## Population          0.000159   0.000368    0.43  0.66533    
## Price              -0.100806   0.007440  -13.55  < 2e-16 ***
## ShelveLocGood       4.848676   0.152838   31.72  < 2e-16 ***
## ShelveLocMedium     1.953262   0.125768   15.53  < 2e-16 ***
## Age                -0.057947   0.015951   -3.63  0.00032 ***
## Education          -0.020852   0.019613   -1.06  0.28836    
## UrbanYes            0.140160   0.112402    1.25  0.21317    
## USYes              -0.157557   0.148923   -1.06  0.29073    
## Income:Advertising  0.000751   0.000278    2.70  0.00729 ** 
## Price:Age           0.000107   0.000133    0.80  0.42381    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.01 on 386 degrees of freedom
## Multiple R-squared:  0.876,  Adjusted R-squared:  0.872 
## F-statistic:  210 on 13 and 386 DF,  p-value: <2e-16

# contrasts shows how R put quantitive variables in linear regression
contrasts(Carseats$ShelveLoc)

##        Good Medium
## Bad       0      0
## Good      1      0
## Medium    0      1

---

Credit

Please note, this material is extracted from online Statistical Learning cource at Stanford University by Prof. T Hastie and Prof R. Tibshirani. It aims only for quick and future references in R and statistical learning. Please visit course page for more information and materials.

---

Linear Algebra and Modern (Markowitz) Portfolio Theory within R

In this text, i demonstrate basic linear algebra, matrix factorization (LU and QR) to solve linear equations and regression, lagrange method, basic numerical methods (bisection and Newton-Raphson) in context of modern (Markowitz, efficient frontier) portfolio theory to estimate optimum weights for given constraints within R. Firstly, lets do some basic linear equation exercises in R.

Load matrixcalc library

library("matrixcalc")

Lets assume we have a linear equation such as: Ax=b where A is square matrix and x is unknown parameters column vector and b is also column vector.

A = matrix(c(2, 2, 4, 9), 2, 2, byrow = T)
A

##      [,1] [,2]
## [1,]    2    2
## [2,]    4    9

b = c(8, 21)

To solve this equation, solve command can be used such as:

solve(A, b)

## [1] 3 1

We can also use various factorization method to resolve this linear equation.

LU Factorization

LU factorization creates elimination matrices (lower and upper) which will be used to solve linear equations.

1. The upper triangular U has the pivots on its diagonal
2. The lower triangular L has ones on its diagonal
3. L has the multipliers \( l_{ij} \) below the diagonal

lu = lu.decomposition(A)
lu

## $L
##      [,1] [,2]
## [1,]    1    0
## [2,]    2    1
## 
## $U
##      [,1] [,2]
## [1,]    2    2
## [2,]    0    5

To solve x via L (lower triangle), U (upper triangle), we have two steps:

First step is to find a column vector c which solves Lc=b equation

c = solve(lu$L, b)
c

## [1] 8 5

In the last step, Ux=c equation is solved which yields x, unknown vector.

x = solve(lu$U, c)
x

## [1] 3 1

QR Factorization

QR factorization decomposites a matrix (A) into two matrices (Q, R) , where A=QR , Q is an orthogonal matrix and R is an upper triangular matrix.

aqr = qr(A)
# orthogonal matrix
aQ = qr.Q(aqr)
aQ

##         [,1]    [,2]
## [1,] -0.4472 -0.8944
## [2,] -0.8944  0.4472

# upper triangular matrix
aR = qr.R(aqr)
aR

##        [,1]   [,2]
## [1,] -4.472 -8.944
## [2,]  0.000  2.236

Since orthogonal matrix has properties of : \( QQ^T=I \), Ax=B equation can be written as: \( Rx=Q^Tb \) where , first A decomposited into QR and then mupltiplied by \( Q^T \) .

aQT = t(aQ)
qb = aQT %*% b
x = backsolve(aR, qb)
x

##      [,1]
## [1,]    3
## [2,]    1

# all of above steps can be done with a one command too
x = qr.solve(A, b)
x

## [1] 3 1

QR factorization for Linear regression

QR factorization is also useful for linear regression fit (\( y=\beta x+\alpha \)) For example, lets try to estimate coefficients (\( \beta \) and \( \alpha \)) of linear regression between Apple and Google's montly returns.

library(quantmod)

# download price of Google and apple
getSymbols(c("GOOG", "AAPL"))

## [1] "GOOG" "AAPL"

google = c(coredata(monthlyReturn(GOOG)))
apple = c(coredata(monthlyReturn(AAPL)))
# let plot returns
plot(google, apple)

Lets apply QR factorization now to find out linear regression (apple~google)

x = cbind(1, google)
xqr = qr(x)

xQ = qr.Q(xqr, complete = T)
xR = qr.R(xqr, complete = T)
# Compute u = QTy
u = t(xQ) %*% apple
# and finally solve for $ \beta $ and $ \alpha $
backsolve(xR[1:2, 1:2], u[1:2])

## [1] 0.01636 0.66522

# lets verify the result with lm command
lm(apple ~ google)

## 
## Call:
## lm(formula = apple ~ google)
## 
## Coefficients:
## (Intercept)       google  
##      0.0164       0.6652

Linear Algebra in Portfolio Theory

Linear equations are common in modern (markowitz, efficient frontier) portfolio theory in finance to estimate optimal weights subject to some constrains, such as maximize return or minimize variance. For example:

\[ maximize: \mu^{T}w \] \[ subject to: e^{T}w=1 \] \[ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ w^T\Sigma w=\sigma^2 \]

where we have following return vector, target risk and covariance matrix :

mu = c(0.08, 0.1, 0.13, 0.15, 0.2)
mu

## [1] 0.08 0.10 0.13 0.15 0.20

s = 0.25^2  # target risk 
s

## [1] 0.0625

Sigma = matrix(c(0.0196, -0.00756, 0.01288, 0.00875, -0.0098, -0.00756, 0.0324, 
    -0.00414, -0.009, 0.00945, 0.01288, -0.00414, 0.0529, 0.020125, 0.020125, 
    0.00875, -0.009, 0.020125, 0.0625, -0.013125, -0.0098, 0.00945, 0.020125, 
    -0.013125, 0.1225), 5, 5, byrow = T)
Sigma

##          [,1]     [,2]     [,3]     [,4]     [,5]
## [1,]  0.01960 -0.00756  0.01288  0.00875 -0.00980
## [2,] -0.00756  0.03240 -0.00414 -0.00900  0.00945
## [3,]  0.01288 -0.00414  0.05290  0.02013  0.02013
## [4,]  0.00875 -0.00900  0.02013  0.06250 -0.01312
## [5,] -0.00980  0.00945  0.02013 -0.01312  0.12250

To find optimum weights, above constraint problem needs to be solved. We will use Lagrangian method for this purpose. Lagrangian is:

\[ F(w,\lambda)= \mu^{T}w + \lambda_1(e^{T}w-1)+\lambda_2(w^T\Sigma w-\sigma^2) \]

To maximize return, we can take first partial derivative of $ F'(w, \lambda)=0$ equation:

\[ G(w, \lambda) = F'(w, \lambda) = \left| \begin{array}{c} \mu+\lambda_1 e + 2\lambda_2 w \Sigma \\ e^{T}w-1 \\ w^T\Sigma w-\sigma^2 \end{array} \right| = 0 \]

To solve this quadratic equation, we can deploy numerical methods such as bisection or Newton-Raphson method. Let me explain briefly, what are these methods:

Bisection Method

If a function f(x) is continuous between [a,b] and if f(a) and f(b) have different sign, then by using following algorithm f(x)=0 can be found:

1. Increase count= count+1
2. Compute f( c ) where c =(a + b)/2 is the midpoint of [a, b]
3. If sign(f( c )) == sign(f(a)), let a=c, otherwise b=c;
4. If |b-a| > tolerance and count < totalIteration then goto 1

R implementation would be:

bisection <- function(f, a, b, tol = 0.001, iteration = 100) {
    count = 0
    while (b - a > tol & count < iteration) {
        c <- (a + b)/2
        if (sign(f(c)) == sign(f(a))) 
            a <- c else b <- c
        count = count + 1
    }
    (a + b)/2
}

This is analogous to binary search algorithm in computer science.

Newton-Raphson method

Newton-Raphson method utilises fundemantal theory of derivatives in order to find root of continuous function f(x)=0 in a recursive function, with a given starting point:

\[ x_{k+1}=x_k- \frac{f(x_k)}{f'(x_k)} \]

Compare to bisection, Newton-Raphson method converge faster, if ever converges.

Lagrange's Method & Newton-Raphson Method

Armed with Newton-Raphson method, now we can solve \( G(w, \lambda) = F'(w, \lambda) \). Algorithm to solve it:

1. Compute G(x) and G'(x)
2. Pick a starting point ($x_0$)
3. Solve the linear system : $G'(x_k)u = G(x_k)$
4. Update $x_{k+1} = x_k − u$
5. Repeat steps 3 and 4 with a number of times or $x_{k+1}- x_k$ is less than a predefined tolerance threshold.

First, we need to get partial derivative of G which will be used in in Newton-Raphson method:

\[ G'(w, \lambda) = \left| \begin{array}{ccc} 2\lambda_2 \Sigma & e & 2 \Sigma w\\ e^{T} & 0 & 0 \\ 2\Sigma w & 0 & 0 \end{array} \right| = 0 \]

Function to compute $ G(w, \lambda) $

G <- function(x, mu, Sigma, sigmaP2) {
    n <- length(mu)
    c(mu + rep(x[n + 1], n) + 2 * x[n + 2] * (Sigma %*% x[1:n]), sum(x[1:n]) - 
        1, t(x[1:5]) %*% Sigma %*% x[1:5] - sigmaP2)
}

Derivative of G

DG <- function(x, mu, Sigma, sigmaP2) {
    n <- length(mu)
    grad <- matrix(0, n + 2, n + 2)
    grad[1:n, 1:n] <- 2 * x[n + 2] * Sigma
    grad[1:n, n + 1] <- 1
    grad[1:n, n + 2] <- 2 * (Sigma %*% x[1:5])
    grad[n + 1, 1:n] <- 1
    grad[n + 2, 1:n] <- 2 * t(x[1:5]) %*% Sigma
    grad
}

Initial weights and $ \lambda $:

x = c(rep(0.5, 5), 1, 1)
x

## [1] 0.5 0.5 0.5 0.5 0.5 1.0 1.0

Lets apply Newton-Raphson iteration now:

for (i in 1:100) {
    x <- x - solve(DG(x, mu, Sigma, s), G(x, mu, Sigma, s))
}

and numerical solution:

x

## [1] -0.39550  0.09606  0.04584  0.70988  0.54372 -0.09201 -0.85715

We need to verify that $ G'(x) $ is positive definite, in order to have a maximum critical points. Lets first see what is $ G'(x) $

DG(x, mu, Sigma, sigmaP2)[1:5, 1:5]

##          [,1]      [,2]      [,3]     [,4]    [,5]
## [1,] -0.03360  0.012960 -0.022080 -0.01500  0.0168
## [2,]  0.01296 -0.055543  0.007097  0.01543 -0.0162
## [3,] -0.02208  0.007097 -0.090686 -0.03450 -0.0345
## [4,] -0.01500  0.015429 -0.034500 -0.10714  0.0225
## [5,]  0.01680 -0.016200 -0.034500  0.02250 -0.2100

Since $ G'(x) $ is symetric and square, all negative eigen value of $ G'(x) $ confirms $G'(x) is $negative definite:

eigen(DG(x, mu, Sigma, s)[1:5, 1:5])$values

## [1] -0.02024 -0.05059 -0.05806 -0.14445 -0.22364

Since all eigen values are negative, we can conclude that x[1:5] is optimum weights for the constraint and maximum return is:

t(x[1:5]) %*% mu

##        [,1]
## [1,] 0.1992

Conjugate Prior, Poisson, Exponential and Gamma Distributions for Website statistics within R

Poisson , Exponential and Gamma distributions are ideal distributions to model number of visitors, arrival time between each visitors and number of visitors between an interval for Website statistics. Let say, if a website attracts 1000 visitors in a day, we can use Poisson distributions to find the probability of k number of visitors in a day. PMF of Poisson is distribution:

$$ pmf(k)= \frac{\lambda^k}{k!} e^{-\lambda} $$

Lets use R to estimate some probabilities with poisson distribution:

# lets assume we have 1000 visitors per day 
>rate=1000

# What is the prob of having exactly 700 visitors in a day
> dpois(700,lambda=rate)
[1] 2.095737e-24

# What is the prob of having no visitors at all in a day
> dpois(0,lambda=rate)
[1] 0

# What is the prob of having less than 700 visitors in a day
> ppois(700,lambda=rate,lower=TRUE)
[1] 6.93301e-24

# What is the prob of having more than 1100 visitors in a day
> ppois(1100,lambda=rate,lower=FALSE) # we are interested in upper tail now, 1100 > rate
[1] 0.000867641

# Maximum arrivals with at least 95% confidence
> qpois(0.95,lambda=rate)
[1] 1052

On the other hand, Exponential distribution is used to estimate probabilities of hit times for a given Poisson distribution. The probability density function (pdf) of an exponential distribution is

$$ f(x;\lambda) = \begin{cases} \lambda e^{-\lambda x} & x \ge 0, \\ 0 & x < 0. \end{cases} $$

The cumulative distribution function is given by

$$ F(x;\lambda) = \begin{cases} 1-e^{-\lambda x} & x \ge 0, \\ 0 & x < 0. \end{cases} $$

Lets use R to estimate some probabilities with Exponential distribution:

# lets assume we have 1000 visitors per day, then average number of visitors per minute time is 1.44:
>AvgVisitorNumberDay=1000
>avgMin=24*60/AvgVisitorNumberDay
> avgMin
[1] 1.44  # 1.44 visitor per minute
# then rate (lambda) for Exponential distribution is (lambda=1/E(X))
> rate=1/avgMin
> rate
[1] 0.6944444 

# lets find the probability that first hit happens at most 45 sec? (P(<45/60) = 1- e^(-rate*60/45))
> pexp(45/60,rate)
[1] 0.4059747

# How long do we have to wait maximum to observe first hit with 95% confidence
> qexp(0.95,rate)
[1] 4.313854 # 258.6 seconds

Another useful distribution is Gamma ($ \alpha, \beta $) distribution which can be used to model the time required for $ \alpha $ events to happen given a Poisson process with mean time $ \beta $. In another words, it is waiting times until a certain number of events happen. For example, Gamma(shape=5, rate=1/3) is the distribution of the length of time (in min) you’d expect to have to wait for 5 visitor hit, given that in average 3 visitors arrive per min. Another example is an insurance company observes that large commercial fire claims occur randomly in time with a mean of 0.7 years between claims. For its financial planning it would like to estimate how long it will be before it pays out the 5th such claim, The time is given by Gamma(5,0.7). Lets have a look at some examples in R.

> rate  # number of hits in a minute , estimated above
[1] 0.6944444 

# pdf of 1 events given rate and have to wait 5 
dgamma(1,rate=rate,shape=5)

# What is the probability that site admin have to wait between 2 to 4 minutes before 5 visitors arrive to website.  
> pgamma(4,rate=rate,shape=5)-pgamma(2,rate=rate,shape=5)
[1] 0.1350604

As you see above, it is trivial to estimate probabilities if parameters (rate) are known. If true parameters are not known, conjugate distribution can be used to estimate parameters. For example, to estimate $ \lambda $ in Poisson distribution, we can use Gamma ($ \alpha, \beta $) distribution as conjugate prior.

Lets assume, we want to estimate $ \lambda $ of Poisson distribution with real-time hit data, with number of visits in fixed time interval on the fly. These are the steps:

- Prior, hyperparameters initialize: Decide on $ \alpha $ and $ \beta $ of Gamma ($ \alpha, \beta $). Best estimator for these parameters are mean and standard deviation of historical arrivals. Set $ \alpha' = \mu $ and $ \beta' = \sigma $
- Posterior hyperparameters update : In each fixed time period $ t $, find number of web page hits $ x_t $ and then update hyperparameters such as $ \alpha' = \alpha' + x_t, \beta' = \beta'+1 $
- Posterior predictive: $\lambda$ is mean of $\operatorname{NB}(\tilde{x}|\alpha',\frac{1}{1+\beta'})$ where NM is Negative binomial distribution. Or more simply:

$$ \lambda= \frac{\alpha'}{\beta'}$$

Similarly we can estimate $ \lambda $ of Exponential distribution with real-time hit data by using number of visits and interval between each visit on the fly. These are the steps:

- Prior, hyperparameters initialize: Decide on $ \alpha $ and $ \beta $ of Gamma ($ \alpha, \beta $). Best estimator for these parameters are mean and standard deviation of historical arrivals. Set $ \alpha' = \mu $ and $ \beta' = \sigma $
- Posterior hyperparameters update : With each page visit and estimate time interval ($ \delta$), and then update parameters such as $ \alpha' = \alpha' + 1, \beta' = \beta'+ \delta $
- Posterior predictive: $\lambda$ is mean of $\operatorname{Lomax}(\tilde{x}|\beta',\alpha')$ where Lomax is Lomax distribution. Or more simply:

$$ \lambda= \frac{\beta'}{\alpha -1'}$$