text:
# Summer School, Big Data and Economics

This week I will be giving a lecture at the  2018 edition of the Summer School at the UB School of Economics, in Barcelona. It will be a four day crash course, starting on Tuesday (morning).

Lecture 1: Introduction : Why Big Data brings New Questions
Lecture 2: Simulation Based Techniques & Bootstrap
Lecture 3: Loss Functions : from OLS to Quantile Regression
Lecture 4: Nonlinearities and Discontinuities
Lecture 5: Cross-Validation and Out-of-Sample diagnosis
Lecture 6: Variable and model selection
Lecture 7: New Tools for Classification Problems
Lecture 8: New Tools for Time Series & Forecasting

Some slides are available on github, and probably more interesting, I will upload a R markdown with all the codes.

# Classification from scratch, logistic with splines 2/8

Today, second post of our series on classification from scratch, following the brief introduction on the logistic regression.

## Piecewise linear splines

To illustrate what’s going on, let us start with a “simple” regression (with only one explanatory variable). The underlying idea is natura non facit saltus, for “nature does not make jumps”, i.e. process governing equations for natural things are continuous. That seems to be a rather strong assumption, because we can assume that there is a fixed threshold to explain death. For instance, if patients die (for sure) if the “stroke index” exceeds a threshold, we might expect some discontinuity. Exceept that if that threshold is an heterogeneous (non-observable continuous) variable, then we get back to the continuity assumption.

The most simple model we can think of to extend the linear model we’ve seen in the previous post is to consider a piecewise linear function, with two parts : small values of $x$, and larger values of $x$. The most convenient way to do so is to use the positive part function $(x-s)_+$ which is the difference between $x$ and $s$ if that difference is positive, and $0$ otherwise. For instance $$\beta_1 x+\beta_2(x-s)_+$$ is the following piecewise linear function, continuous, with a “rupture” at knot $s$.

Observe also the following interpretation: for small values of $x$, there is a linear increase, with slope $\beta_1$, and for lager values of $x$, there is a linear decrease, with slope $\beta_1+\beta_2$. Hence, $\beta_2$ is interpreted as a change of the slope.

And of course, it is possible to consider more than one knot. The function to get the positive value is the following

pos = function(x,s) (x-s)*(x>=s)

then we can use it direcly in our regression model

reg = glm(PRONO~INSYS+pos(INSYS,15)+ pos(INSYS,25),data=myocarde,family=binomial)

The output of the regression is here

summary(reg)   Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.1109 3.2783 -0.034 0.9730 INSYS -0.1751 0.2526 -0.693 0.4883 pos(INSYS, 15) 0.7900 0.3745 2.109 0.0349 * pos(INSYS, 25) -0.5797 0.2903 -1.997 0.0458 *

Hence, the original slope, for very small values is not significant, but then, above 15, it become significantly positive. And above 25, there is a significant change again. We can plot it to see what’s going on

u = seq(5,55,length=201) v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,type="l") points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=c(5,15,25,55),lty=2)

## Using bs() linear splines

Using the GAM function, things are slightly different. We will use here so called b-splines,

library(splines)

We can define spline functions with support $(5,55)$ and with knots $\{15,25\}$

clr6 = c("#1b9e77","#d95f02","#7570b3","#e7298a","#66a61e","#e6ab02") x = seq(0,60,by=.25) B = bs(x,knots=c(15,25),Boundary.knots=c(5,55),degre=1) matplot(x,B,type="l",lty=1,lwd=2,col=clr6)

as we can see, the functions defined here are different from the one before, but we still have (piecewise) linear functions on each segment $(5,15)$, $(15,25)$ and $(25,55)$. But linear combinations of those functions (the two sets of functions) will generate the same space. Said differently, if the interpretation of the output will be different, predictions should be the same

reg = glm(PRONO~bs(INSYS,knots=c(15,25), Boundary.knots=c(5,55),degre=1), data=myocarde,family=binomial) summary(reg)   Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.9863 2.0555 -0.480 0.6314 bs(INSYS,..)1 -1.7507 2.5262 -0.693 0.4883 bs(INSYS,..)2 4.3989 2.0619 2.133 0.0329 * bs(INSYS,..)3 5.4572 5.4146 1.008 0.3135

Observe that there are three coefficients, as before, but again, the interpretation is here more complicated…

v=predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red") points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=c(5,15,25,55),lty=2)

Nevertheless, the prediction is the same… and that’s nice.

Let us go one step further… Can we have also the continuity of the derivative ? Yes, and that’s easy actually, considering parabolic functions. Instead of using a decomposition on $x,(x-s_1)_+$ and $(x-s_2)_+$ consider now a decomposition on $x,x^{\color{red}{2}},(x-s_1)^{\color{red}{2}}_+$ and $(x-s_2)^{\color{red}{2}}_+$.

 pos2 = function(x,s) (x-s)^2*(x>=s) reg = glm(PRONO~poly(INSYS,2)+pos2(INSYS,15)+pos2(INSYS,25), data=myocarde,family=binomial) summary(reg)   Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 29.9842 15.2368 1.968 0.0491 * poly(INSYS, 2)1 408.7851 202.4194 2.019 0.0434 * poly(INSYS, 2)2 199.1628 101.5892 1.960 0.0499 * pos2(INSYS, 15) -0.2281 0.1264 -1.805 0.0712 . pos2(INSYS, 25) 0.0439 0.0805 0.545 0.5855

As expected, there are here five coefficients: the intercept and two for the part on the left (three parameters for the parabolic function), and then two additional terms for the part in the center – here $(15,25)$ – and for the part on the right. Of course, for each portion, there is only one degree of freedom since we have a parabolic function (three coefficients) but two constraints (continuity, and continuity of the first order derivative).

On a graph, we get the following

v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red",lwd=2,xlab="INSYS",ylab="") points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=c(5,15,25,55),lty=2)

Of course, we can do the same with our R function. But as before, the basis of function is expressed here differently

 x = seq(0,60,by=.25) B=bs(x,knots=c(15,25),Boundary.knots=c(5,55),degre=2) matplot(x,B,type="l",xlab="INSYS",col=clr6)

If we run R code, we get

reg = glm(PRONO~bs(INSYS,knots=c(15,25), Boundary.knots=c(5,55),degre=2),data=myocarde, family=binomial) summary(reg)   Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 7.186 5.261 1.366 0.1720 bs(INSYS, ..)1 -14.656 7.923 -1.850 0.0643 . bs(INSYS, ..)2 -5.692 4.638 -1.227 0.2198 bs(INSYS, ..)3 -2.454 8.780 -0.279 0.7799 bs(INSYS, ..)4 6.429 41.675 0.154 0.8774

But that’s not really a big deal since the prediction is exactly the same

v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red") points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=c(5,15,25,55),lty=2)

## Cubic splines

Last, but not least, we can reach the cubic splines. With our previous notions, we would consider a decomposition on (guess what) $x,x^2,x^{\color{red}{3}},(x-s_1)^{\color{red}{3}}_+,(x-s_2)^{\color{red}{3}}_+$, to get this time continuity, as well as continuity of the first two derivatives (and to get a very smooth function, since even variations will be smooth). If we use the bs function, the basis is the followin

B=bs(x,knots=c(15,25),Boundary.knots=c(5,55),degre=3) matplot(x,B,type="l",lwd=2,col=clr6,lty=1,ylim=c(-.2,1.2)) abline(v=c(5,15,25,55),lty=2)

and the prediction will now be

reg = glm(PRONO~bs(INSYS,knots=c(15,25), Boundary.knots=c(5,55),degre=3), data=myocarde,family=binomial) u = seq(5,55,length=201) v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red",lwd=2) points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=c(5,15,25,55),lty=2)

Two last things before concluding (for today), the location of the knots, and the extension to additive models.

## Location of knots

In many applications, we do not want to specify the location of the knots. We just want – say – three (intermediary) knots. This can be done using

reg = glm(PRONO~1+bs(INSYS,degree=1,df=4),data=myocarde,family=binomial)

We can actually get the locations of the knots by looking at

attr(reg$terms, "predvars")[[3]] bs(INSYS, degree = 1L, knots = c(15.8, 21.4, 27.15), Boundary.knots = c(8.7, 54), intercept = FALSE) which provides us with the location of the boundary knots (the minumun and the maximum from from our sample) but also the three intermediary knots. Observe that actually, those five values are just (empirical) quantiles quantile(myocarde$INSYS,(0:4)/4) 0% 25% 50% 75% 100% 8.70 15.80 21.40 27.15 54.00

If we plot the prediction, we get

v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red",lwd=2) points(myocarde$INSYS,myocarde$PRONO,pch=19) abline(v=quantile(myocarde$INSYS,(0:4)/4),lty=2) If we get back on what was computed before the logit transformation, we clealy see ruptures are the different quantiles B = bs(x,degree=1,df=4) B = cbind(1,B) y = B%*%coefficients(reg) plot(x,y,type="l",col="red",lwd=2) abline(v=quantile(myocarde$INSYS,(0:4)/4),lty=2)

Note that if we do specify anything about knots (number or location), we get no knots…

reg = glm(PRONO~1+bs(INSYS,degree=2),data=myocarde,family=binomial) attr(reg$terms, "predvars")[[3]] bs(INSYS, degree = 2L, knots = numeric(0), Boundary.knots = c(8.7,54), intercept = FALSE) and if we look at the prediction u = seq(5,55,length=201) v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red",lwd=2) points(myocarde$INSYS,myocarde$PRONO,pch=19) actually, it is the same as a quadratic regression (as expected actually) reg = glm(PRONO~1+poly(INSYS,degree=2),data=myocarde,family=binomial) v = predict(reg,newdata=data.frame(INSYS=u),type="response") plot(u,v,ylim=0:1,type="l",col="red",lwd=2) points(myocarde$INSYS,myocarde$PRONO,pch=19) ## Additive models Consider now the second dataset, with two variables. Consider here a model like $$\mathbb{P}[Y|X_1=x_1,X_2=x_2]=\frac{\exp[\eta(x_1,x_2)]}{1+\exp[\eta(x_1,x_2)]}$$ where $$\exp[\eta(x_1,x_2)]=\beta_0+\color{red}{s_1(x_1)}+\color{blue}{s_2(x_2)}$$ $$\color{red}{s_1(x_1)}=\beta_{1,0}x_1+\beta_{1,1}(x_1-s_{11})_++\beta_{1,2}(x_1-s_{12})_+$$ and $$\color{blue}{s_2(x_2)}=\beta_{2,0}x_2+\beta_{2,1}(x_2-s_{21})_++\beta_{2,2}(x_2-s_{22})_+$$ It might seem a little bit restrictive, but that’s actually the idea of additive models. reg = glm(y~bs(x1,degree=1,df=3)+bs(x2,degree=1,df=3),data=df,family=binomial(link = "logit")) u = seq(0,1,length=101) p = function(x,y) predict.glm(reg,newdata=data.frame(x1=x,x2=y),type="response") v = outer(u,u,p) image(u,u,v,xlab="Variable 1",ylab="Variable 2",col=clr10,breaks=(0:10)/10) points(df$x1,df$x2,pch=19,cex=1.5,col="white") points(df$x1,df$x2,pch=c(1,19)[1+(df$y=="1")],cex=1.5) contour(u,u,v,levels = .5,add=TRUE)

Now, if think about is, we’ve been able to get a “perfect” model, so, somehow, it seems no longer continuous…

persp(u,u,v,theta=20,phi=40,col="green"

Of course, it is… it is piecewise linear, with hyperplane, some being almost vertical.

And one can also consider piecewise quadratic functions

reg = glm(y~bs(x1,degree=2,df=3)+bs(x2,degree=2,df=3),data=df,family=binomial(link = "logit")) u = seq(0,1,length=101) p = function(x,y) predict.glm(reg,newdata=data.frame(x1=x,x2=y),type="response") v = outer(u,u,p) image(u,u,v,xlab="Variable 1",ylab="Variable 2",col=clr10,breaks=(0:10)/10) points(df$x1,df$x2,pch=19,cex=1.5,col="white") points(df$x1,df$x2,pch=c(1,19)[1+(df$y=="1")],cex=1.5) contour(u,u,v,levels = .5,add=TRUE) Funny thing, we now have two “perfect” models, with different areas for the white and the black dots… Don’t ask me how to choose on that one. In R, it is possible to use the mgcv package to run a gam regression. It is used for generalized additive models, but here, we have only one variable, so it is difficult to see the “additive” part, actually. And to be more specific, mgcv is using penalized quasi-likelihood from the nlme package (but we’ll get back on penalized routines later on). But maybe I should also mention another smoothing tool before, kernels (and maybe also $k$-nearest neighbors). To be continued # Classification from scratch, logistic regression 1/8 Let us start today our series on classification from scratch The logistic regression is based on the assumption that given covariates $\mathbf{x}$, $Y$ has a Bernoulli distribution,$$Y|\mathbf{X}=\mathbf{x}\sim\mathcal{B}(p_{\mathbf{x}}),~~~~p_\mathbf{x}=\frac{\exp[\mathbf{x}^T\mathbf{\beta}]}{1+\exp[\mathbf{x}^T\mathbf{\beta}]}$$The goal is to estimate parameter $\mathbf{\beta}$. Recall that the heuristics for the use of that function for the probability is that$$\log[\text{odds}(Y=1)]=\log\frac{\mathbb{P}[Y=1]}{\mathbb{P}[Y=0]}=\mathbf{x}^T\mathbf{\beta}$$ ## Maximimum of the (log)-likelihood function The log-likelihood is here$$\log\mathcal{L} = \sum_{i=1}^n y_i\log p_i+(1-y_i)\log (1-p_i)$$ where $p_{i}=(1+\exp[-\mathbf{x}_i^T\mathbf{\beta}])^{-1}$. Numerical techniques are based on (numerical) gradient descent to compute the maximum of the likelihood function. The (negative) log-likelihood is the following function y = myocarde$PRONO X = cbind(1,as.matrix(myocarde[,1:7])) negLogLik = function(beta){ -sum(-y*log(1 + exp(-(X%*%beta))) - (1-y)*log(1 + exp(X%*%beta))) }

We use the minus sign since standard optimization routines compute minima, not maxima. Now, to find the minimum of that function, we need a starting point to initiate the algorithm

beta_init = lm(PRONO~.,data=myocarde)$coefficients Why not start with the parameter of the OLS. Somehow, we might think that at least, sign should be ok for instance. Anyway, we need a starting point, and let us use that one. logistic_opt = optim(par = beta_init, negLogLik, hessian=TRUE, method = "BFGS", control=list(abstol=1e-9)) Here, we obtain  logistic_opt$par (Intercept) FRCAR INCAR INSYS 1.656926397 0.045234029 -2.119441743 0.204023835 PRDIA PAPUL PVENT REPUL -0.102420095 0.165823647 -0.081047525 -0.005992238

Let us verify here that this output is valid. For instance, what if we change the value of the starting point (randomly)

simu = function(i){ logistic_opt_i = optim(par = rnorm(8,0,3)*beta_init, negLogLik, hessian=TRUE, method = "BFGS", control=list(abstol=1e-9)) logistic_opt_i$par[2:3] } v_beta = t(Vectorize(simu)(1:1000)) plot(v_beta) par(mfrow=c(1,2)) hist(v_beta[,1],xlab=names(myocarde)[1]) hist(v_beta[,2],xlab=names(myocarde)[2]) Ooops. There is a problem here. Clearly, we cannot rely on numerical optimization here. We can think about using another optimization routine library(optimx) logit = function(mX, vBeta) { exp(mX %*% vBeta)/(1+ exp(mX %*% vBeta)) } logLikelihoodLogitStable = function(vBeta, mX, vY) { -sum(vY*(mX %*% vBeta - log(1+exp(mX %*% vBeta))) + (1-vY)*(-log(1 + exp(mX %*% vBeta)))) } likelihoodScore = function(vBeta, mX, vY) { return(t(mX) %*% (logit(mX, vBeta) - vY) ) } optimLogitLBFGS = optimx(beta_init, logLikelihoodLogitStable, method = 'L-BFGS-B', gr = likelihoodScore, mX = X, vY = y, hessian=TRUE) The optimum is here attr(optimLogitLBFGS, "details")[[2]] [,1] 0.066680272 FRCAR 0.003080542 INCAR 0.079031364 INSYS -0.001586194 PRDIA 0.040500697 PAPUL -0.041870705 PVENT -0.014162756 REPUL 0.195632244 Let’s be honest here, I do not feel confortable with those techniques. So, what happened here ? Here, the technique we use is based on the following idea,$$\mathbf{\beta}_{new}=\mathbf{\beta}_{old} -\left(\frac{\partial^2\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}\partial\mathbf{\beta}^T}\right)^{-1}\cdot \frac{\partial\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}}$$The problem is that my computer does not know this first and second derivatives. So it will compute them using approximation techniques. Actually, it is possible to use functions dedicated to such computation library(numDeriv) library(MASS) logit = function(x){1/(1+exp(-x))} logLik = function(beta, X, y){ -sum(y*log(logit(X%*%beta)) + (1-y)*log(1-logit(X%*%beta))) } optim_second = function(beta, num_iter){ LL = vector() for(i in 1:num_iter){ grad = (t(X)%*%(logit(X%*%beta) - y)) H = hessian(logLik, beta, method = "complex", X = X, y = y) beta = beta - ginv(H)%*%grad LL[i] = logLik(beta, X, y) } result = list(beta, H) return(result) } With our OLS starting point, we obtain opt0 = optim_second(beta_init,500) opt0[[1]] [,1] [1,] 0.951074420 [2,] 0.018860280 [3,] 0.275428978 [4,] 0.144803636 [5,] -0.058535606 [6,] 0.001182178 [7,] -0.108651776 [8,] -0.002940315 But if we try with another starting point opt1 = optim_second(beta_init*runif(8),500) opt1[[1]] [,1] [1,] 0.052894794 [2,] 0.024718435 [3,] 0.167953661 [4,] 0.171662947 [5,] -0.057458066 [6,] -0.011361034 [7,] -0.107532114 [8,] -0.002679064 Clearly, some coefficients are rather close. But other aren’t. From my point of viezw, that is a major problem (keep in mind that we do not deal here with massive data ! There are only 7 explanatory variables, and only 71 observations). Why not try to be clever, and use the analytical values of those derivatives ? Even if some people claim the oppositive, sometimes, it can actually be usefull to do the maths, instead of considering only numerical values. ## Newton (or Fisher) Algorithm If you open any Econometrics textbooks (one can also try to derive it), you will get $$\frac{\partial\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}}=\mathbf{X}^T(\mathbf{y}-\mathbf{p}_{old})$$ while$$\frac{\partial^2\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}\partial\mathbf{\beta}^T}=-\mathbf{X}^T\mathbf{\Delta}_{old}\mathbf{X}$$ Y=myocarde$PRONO X=cbind(1,as.matrix(myocarde[,1:7])) colnames(X)=c("Inter",names(myocarde[,1:7])) beta=as.matrix(lm(Y~0+X)$coefficients,ncol=1) for(s in 1:9){ pi=exp(X%*%beta[,s])/(1+exp(X%*%beta[,s])) gradient=t(X)%*%(Y-pi) omega=matrix(0,nrow(X),nrow(X));diag(omega)=(pi*(1-pi)) Hessian=-t(X)%*%omega%*%X beta=cbind(beta,beta[,s]-solve(Hessian)%*%gradient)} Observe that here, I use only ten iterations of the algorithm !  beta[,8:10] [,1] [,2] [,3] XInter -10.187641685 -10.187641696 -10.187641696 XFRCAR 0.138178119 0.138178119 0.138178119 XINCAR -5.862429035 -5.862429037 -5.862429037 XINSYS 0.717084018 0.717084018 0.717084018 XPRDIA -0.073668171 -0.073668171 -0.073668171 XPAPUL 0.016756506 0.016756506 0.016756506 XPVENT -0.106776012 -0.106776012 -0.106776012 XREPUL -0.003154187 -0.003154187 -0.003154187 The thing is that is seems to converge extremely fast. And it is rather robust ! Look at what we get if we change our starting point beta=as.matrix(lm(Y~0+X)$coefficients,ncol=1)*runif(8) for(s in 1:9){ pi=exp(X%*%beta[,s])/(1+exp(X%*%beta[,s])) gradient=t(X)%*%(Y-pi) omega=matrix(0,nrow(X),nrow(X));diag(omega)=(pi*(1-pi)) Hessian=-t(X)%*%omega%*%X beta=cbind(beta,beta[,s]-solve(Hessian)%*%gradient)} beta[,8:10] [,1] [,2] [,3] XInter -10.187641586 -10.187641696 -10.187641696 XFRCAR 0.138178118 0.138178119 0.138178119 XINCAR -5.862429017 -5.862429037 -5.862429037 XINSYS 0.717084013 0.717084018 0.717084018 XPRDIA -0.073668172 -0.073668171 -0.073668171 XPAPUL 0.016756508 0.016756506 0.016756506 XPVENT -0.106776012 -0.106776012 -0.106776012 XREPUL -0.003154187 -0.003154187 -0.003154187

Nice, isn’t it? Looks like we got our winner, don’t we? And one can use the inverse of the Hessian matrix to get standard deviations.

## Weighted Least-Squares

Let us go one step further. We’ve seen that we want to compute something like$$\mathbf{\beta}_{new} =(\mathbf{X}^T\mathbf{\Delta}_{old}\mathbf{X})^{-1}\mathbf{X}^T\mathbf{\Delta}_{old}\mathbf{z}$$(if we do substitute matrices in the analytical expressions) where $\mathbf{z}=\mathbf{X}\mathbf{\beta}_{old}+\mathbf{\Delta}_{old}^{-1}[\mathbf{y}-\mathbf{p}_{old}]$. But actually, that’s simply a standard least-square problem$$\mathbf{\beta}_{new} = \text{argmin}\left\lbrace(\mathbf{z}-\mathbf{X}\mathbf{\beta})^T\mathbf{\Delta}_{old}^{-1}(\mathbf{z}-\mathbf{X}\mathbf{\beta})\right\rbrace$$The only problem here is that weights $\mathbf{\Delta}_{old}$ are functions of unknown $\mathbf{\beta}_{old}$. But actually, if we keep iterating, we should be able to solve it : given the $\mathbf{\beta}$ we got the weights, and with the weights, we can use weighted OLS to get an updated $\mathbf{\beta}$. That’s the idea of iteratively reweighted least squares.

The algorithm will be

df = myocarde beta_init = lm(PRONO~.,data=df)$coefficients X = cbind(1,as.matrix(myocarde[,1:7])) beta = beta_init for(s in 1:1000){ p = exp(X %*% beta) / (1+exp(X %*% beta)) omega = diag(nrow(df)) diag(omega) = (p*(1-p)) df$Z = X %*% beta + solve(omega) %*% (df$PRONO - p) beta = lm(Z~.,data=df[,-8], weights=diag(omega))$coefficients }

and the output is here

 beta (Intercept) FRCAR INCAR INSYS PRDIA -10.187641696 0.138178119 -5.862429037 0.717084018 -0.073668171 PAPUL PVENT REPUL 0.016756506 -0.106776012 -0.003154187

which is almost what we’ve obtained before. Nice isn’t it ? Actually, here we also have standard deviations of estimators

summary( lm(Z~.,data=df[,-8], weights=diag(omega)))   Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -10.187642 10.668138 -0.955 0.343 FRCAR 0.138178 0.102340 1.350 0.182 INCAR -5.862429 6.052560 -0.969 0.336 INSYS 0.717084 0.503527 1.424 0.159 PRDIA -0.073668 0.261549 -0.282 0.779 PAPUL 0.016757 0.306666 0.055 0.957 PVENT -0.106776 0.099145 -1.077 0.286 REPUL -0.003154 0.004386 -0.719 0.475

## The standard glm function

Of course, it is possible to use an R built-in function to get our estimate

summary(glm(PRONO~.,data=myocarde,family=binomial(link = "logit")))   Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -10.187642 11.895227 -0.856 0.392 FRCAR 0.138178 0.114112 1.211 0.226 INCAR -5.862429 6.748785 -0.869 0.385 INSYS 0.717084 0.561445 1.277 0.202 PRDIA -0.073668 0.291636 -0.253 0.801 PAPUL 0.016757 0.341942 0.049 0.961 PVENT -0.106776 0.110550 -0.966 0.334 REPUL -0.003154 0.004891 -0.645 0.519

## Application and visualisation

Let us visualize the prediction obtained from the logistic regression, on our second dataset

x = c(.4,.55,.65,.9,.1,.35,.5,.15,.2,.85) y = c(.85,.95,.8,.87,.5,.55,.5,.2,.1,.3) z = c(1,1,1,1,1,0,0,1,0,0) df = data.frame(x1=x,x2=y,y=as.factor(z)) reg = glm(y~x1+x2,data=df,family=binomial(link = "logit")) u = seq(0,1,length=101) p = function(x,y) predict.glm(reg,newdata=data.frame(x1=x,x2=y),type="response") v = outer(u,u,p) image(u,u,v,xlab="Variable 1",ylab="Variable 2",col=clr10,breaks=(0:10)/10) points(x,y,pch=19,cex=1.5,col="white") points(x,y,pch=c(1,19)[1+z],cex=1.5) contour(u,u,v,levels = .5,add=TRUE)

Here level curves – or iso-probabilities – are linear, so the space is divided in two (0 and 1, survival and death, white and black) by a straight line (or an hyperplane in higher dimension). Furthermore, since we have a linear model, if we change the cutoff (the threshold used to create the two classes), we obtain another straight line (or hyperplane) parallel to the first one.

Next time, we will introduce splines to smooth those continuous covariates… to be continued.

# Classification from scratch, overview 0/8

Before my course on « big data and economics » at the university of Barcelona in July, I wanted to upload a series of posts on classification techniques, to get an insight on machine learning tools.

According to some common idea, machine learning algorithms are black boxes. I wanted to get back on that saying. First of all, isn’t it the case also for regression models, like generalized additive models (with splines) ? Do you really know what the algorithm is doing ? Even the logistic regression. In textbooks, we can easily find math formulas. But what is really done when I run it, in R ?

When I started working on academia, someone told me something like « if you really want to understand a theory, teach it ». And that has been my moto for more than 15 years. I wanted to add a second part to that statement: « if you really want to understand an algorithm, recode it ». So let’s try this… My ambition is to recode (more or less) most of the standard algorithms used in predictive modeling, from scratch, in R. What I plan to mention, within the next two weeks, will be

I will use two datasets to illustrate. The first one is inspired by the cover of « Foundations of Machine Learning » by Mehryar Mohri, Afshin Rostamizadeh and Ameet Talwalkar. At least, with this dataset, it will be possible to plot predictions (since there are only two – continuous – features)

x = c(.4,.55,.65,.9,.1,.35,.5,.15,.2,.85) y = c(.85,.95,.8,.87,.5,.55,.5,.2,.1,.3) z = c(1,1,1,1,1,0,0,1,0,0) df = data.frame(x1=x,x2=y,y=as.factor(z)) plot(x,y,pch=c(1,19)[1+z])

Here is some code to get a visualization of the prediction (here the probability to be a black point)

rmatrix_model = function(model){ u = seq(0,1,length=101) p = function(x,y) predict(model,newdata=data.frame(x1=x,x2=y),type="response") v = outer(u,u,p) return(v)} nice_graph=function(v){ u = seq(0,1,length=101) image(u,u,v,xlab="Variable 1",ylab="Variable 2",col=clr10[c(1,10)],breaks=c(0,5,10)/10) points(x,y,pch=19,cex=1.5,col="white") points(x,y,pch=c(1,19)[1+z],cex=1.5) contour(u,u,v,levels = .5,add=TRUE) } reg = glm(y~x1+x2,data=df,family=binomial) nice_graph(rmatrix_model(reg))

Note that colors are defined here as

clr10= c("#ffffff","#f7fcfd","#e5f5f9","#ccece6","#99d8c9","#66c2a4","#41ae76","#238b45","#006d2c","#00441b")

or with some nonlinear model

The second one is a dataset I got from Gilbert Saporta, about heart attacks and decease (our binary variable).

myocarde = read.table("http://freakonometrics.free.fr/myocarde.csv",head=TRUE, sep=";") myocarde$PRONO = (myocarde$PRONO=="SURVIE")*1 y = myocarde$PRONO X = as.matrix(cbind(1,myocarde[,1:7])) So far, I do not plan to talk (too much) on the choice of tunning parameters (and cross-validation), on comparing models, etc. The goal here is simply to understand what’s going on when we call either glm, glmnet, gam, random forest, svm, xgboost, or any function to get a predict model. # Machine Learning Just a short joke that explains what machine learning is about interviewer: what’s your biggest strength? me: I’m an expert in machine learning. interviewer: What’s 9+10? me: It’s 3. interviewer: Not even close. It’s 19. me: It’s 16, interviewer: Wrong. It’s still 19. me: It’s 18. interviewer: No it’s 19. me: It’s 19. interviewer: You’re hired! # Econometrics vs. Machine Learning with Temporal Patterns A few months ago, I did publish a (long) post entitled ‘some thoughts on economics, mathematics, econometrics, machine learning, etc‘. In that post, I was discussing possible differences between foundations of econometrics, and machine learning. I wanted to get back today on an important point, related to training/sampling datasets, when we have temporal data. I was discussing this morning, with a student of the Data Science for Actuaries program, an interesting point related to claim frequency models, for insurance ratemaking. Since the goal is to predict claims frequency (to assess the level of the insurance premium), he suggested to use old data to train the model, and more recent one to test it. The problem is that the model did not incorporate any temporal pattern, and we got surprising results. Consider here a simple dataset, > set.seed(1) > n=50000 > X1=runif(n) > T=sample(2000:2015,size=n,replace=TRUE) > L=exp(-3+X1-(T-2000)/20) > E=rbeta(n,5,1) > Y=rpois(n,L*E) > B=data.frame(Y,X1,L,T,E) Claims frequency is driven by a Poisson process, with one covariate, X1, and we assume that the intensity decreases (with an exponential rate). Consider here a standard linear regression, without any time effect > reg=glm(Y~X1+offset(log(E)),data=B, + family=poisson) We can also compute the empirical annualized claims frequency > u=seq(0,1,by=.01) > v=predict(reg,newdata=data.frame(X1=u,E=1)) > p=function(x){ + B=B[abs(B$X1-x)<.1,]
+   sum(B$Y)/sum(B$E)
+ }
> vp=Vectorize(p)(seq(.05,.95,by=.1))

and plot the two curves on the same graph,

> plot(seq(.05,.95,by=.1),vp,type="b")
> lines(u,exp(v),lty=2,col="red")

This is what we usually do in econometrics. In machine learning, and more specifically to assess the quality of the model, and for model selection, it is common to split the dataset in two parts. A training sample, and a validation sample. Consider some randomized training/validation samples, then fit a model on the training sample, and finally use it to get a prediction,

> idx=sample(1:nrow(B),size=nrow(B)*7/8)
> B_a=B[idx,]
> B_t=B[-idx,]
> reg=glm(Y~X1+offset(log(E)),data=B_a,
+ family=poisson)
> u=seq(0,1,by=.01)
> v=predict(reg,newdata=data.frame(X1=u,E=1))
> p=function(x){
+   B=B_a[abs(B_a$X1-x)<.1,] + sum(B$Y)/sum(B$E) + } > vp_a=Vectorize(p)(seq(.05,.95,by=.1)) > plot(seq(.05,.95,by=.1),vp_a,col="blue") > lines(u,exp(v),lty=2) > p=function(x){ + B=B_t[abs(B_t$X1-x)<.1,]
+   sum(B$Y)/sum(B$E)
+ }
> vp_t=Vectorize(p)(seq(.05,.95,by=.1))
> lines(seq(.05,.95,by=.1),vp_t,col="red")

The blue curve is the prediction on the training sample (as we usually do in econometrics), but then the red curve is the prediction on the testing sample. Here, volatility probably comes from the small size of the testing sample (1 observation out of 8).

Now, what if we use the year as a splitting criteria : we fit a model on old years to fit a model, and we test it on recent years,

> B_a=subset(B,T<2014)
> B_t=subset(B,T>=2014)
> reg=glm(Y~X1+offset(log(E)),data=B_a,family=poisson)
> u=seq(0,1,by=.01)
> v=predict(reg,newdata=data.frame(X1=u,E=1))
> p=function(x){
+   B=B_a[abs(B_a$X1-x)<.1,] + sum(B$Y)/sum(B$E) + } > vp_a=Vectorize(p)(seq(.05,.95,by=.1)) > plot(seq(.05,.95,by=.1),vp_a,col="blue") > lines(u,exp(v),lty=2) > p=function(x){ + B=B_t[abs(B_t$X1-x)<.1,]
+   sum(B$Y)/sum(B$E)
+ }
> vp_t=Vectorize(p)(seq(.05,.95,by=.1))
> lines(seq(.05,.95,by=.1),vp_t,col="red")

Clearly, we miss something here…

We were looking at such a graph this morning, and it took me some time to understand how training and validation samples were designed, and that there was a possible temporal effect (actually, this morning, it was based on a 3 year training sample, and a 1 year validation sample).

Since there is a temporal pattern, let us capture it. As an econometrician, let me use a regression model

> reg=glm(Y~X1+T+offset(log(E)),data=B,
+ family=poisson)
> C=coefficients(reg)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3)
> plot(2000:2015,exp(C[1]+C[3]*(2000:2015)))
> lines(u,v,lty=2,col="red")

(I focus only on the evolution of the temporal variate on that graph).

Here, we use a linear model, but there are usually no reason to assume linearity. So we might consider splines

> library(splines)
> reg=glm(Y~X1+bs(T)+offset(log(E)),
+ data=B,family=poisson)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3)
> v2=predict(reg,newdata=data.frame(X1=0,
+ T=2000:2015,E=1))
> plot(2000:2015,exp(v2),type="b")
> lines(u,v,lty=2,col="red")

But here again, why should we assume that there is an underlying smooth function? There might be some ruptures… So let us consider a regression on factors

> reg=glm(Y~0+X1+as.factor(T)+offset(log(E)),
+ data=B,family=poisson)
> C=coefficients(reg)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3)
> plot(2000:2015,exp(C[2:17]),type="b")
> lines(u,v,lty=2,col="red")

An alternative might be to consider some more general model, like a regression tree

> library(rpart)
> reg=rpart(Y~X1+T+offset(log(E)),data=B,
+ method="poisson",cp=1e-4)
> p=function(t){
+   B=B[B$T==t,] + B$E=1
+   mean(predict(reg,newdata=B))
+ }
> y_m=Vectorize(function(t) p(t))(2000:2015)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3+.5)
> plot(2000:2015,y_m,ylim=c(.02,.085),type="b")
> lines(u,v,lty=2,col="red")

Here, it seems that something went wrong. I guess it’s coming from the exposure. So consider a simplier model, on the annualized frequency, and with weights that are related to the exposure

> reg=rpart(Y/E~X1+T,data=B,weights=B$E,cp=1e-4) > p=function(t){ + B=B[B$T==t,]
+   B$E=1 + mean(predict(reg,newdata=B)) + } > y_m=Vectorize(function(t) p(t))(2000:2015) > u=seq(1999,2016,by=.1) > v=exp(-(u-2000)/20-3+.5) > plot(2000:2015,y_m,ylim=c(.02,.085),type="b") > lines(u,v,lty=2,col="red") That was for the econometrician perspective. With a machine learning perspective, consider a training sample (here based on old data) and a validation sample (based on more recent ones) > B_a=subset(B,T<2014) > B_t=subset(B,T>=2014) If we consider a model, it is easy to get a prediction on recent years, even if the model was designed to model older ones, > reg_a=glm(Y~X1+T+offset(log(E)), + data=B_a,family=poisson) > C=coefficients(reg_a) > u=seq(1999,2016,by=.1) > v=exp(-(u-2000)/20-3) > plot(2000:2015,exp(C[1]+C[3]*c(2000:2013, + NA,NA)),type="b") > lines(u,v,lty=2,col="red") > points(2014:2015,exp(C[1]+C[3]*2014:2015), + pch=19,col="blue") But if we use years as factors, things are more complicated. > reg_a=glm(Y~0+X1+as.factor(T)+offset(log(E)), + data=B_a,family=poisson) > C=coefficients(reg_a) > RMSE=function(A){ + L=exp(C[1]*B_t$X1+ A[1]*(B_t$T==2014) + A[2]*(B_t$T==2015))
+   Y_t=L*B_t$E + sum( (Y_t - B_t$Y )^2)}
> i=optim(c(.4,.4),RMSE)\$par
> plot(2000:2015,c(exp(C[2:15]),NA,NA),)
> u=seq(1999,2016,by=.1)
> v=exp(-(u-2000)/20-3)
> lines(u,v,lty=2,col="red")
> points(2014:2015,exp(i),pch=19,col="blue")

becase we need to get a prediction on levels that were not in our training sample. Here, we minimize the RMSE to quantify factor levels for recent years. And the output is not that bad.

So yes, it is possible to get a training dataset on older data, and test it on recent years. But one should be careful, and take into account, properly, temporal patterns.

# Big Data and Machine Learning with an Actuarial Perspective

Next week, I will give a (one day) crash course at the IA | BE Sumer School, in Louvain-la-Neuve. Slides are available on the blog.

# Some thoughts on Economics, Mathematics, Econometrics, Statistics, Machine Learning, etc

There were a lot of posts, recently, related to those topics, starting with Noah Smith ‘s piece entitled “Economics has a Math Problem” and more recently “Econometrics, Math, and Machine Learning…what?” by Matt Bogard. I don’t have (yet) a clear mind on those issues, but there are still a few thoughts that I wanted to share. I did not really want to, but I’ve been asked, on Twitter, and I thought it might be good to write them down, to clarify some ideas I have, but also (probably, hopefully) to get interesting feedbacks.

# Data Science: from Small to Big Data

This Tuesday, I will be in Leuven (in Belgium) at the ACP meeting to give  a talk on Data Science: from Small to Big Data. The talk will take place in the Faculty Club from 6 till 8 pm. Slides could be found online (with animated pictures).


text:
match.data() and get_matches() create a data frame with additional variables for the distance measure, matching weights, and subclasses after matching. This dataset can be used to estimate treatment effects after matching or subclassification. get_matches() is most useful after matching with replacement; otherwise, match.data() is more flexible. See Details below for the difference between them.

match.data(object,
group = "all",
distance = "distance",
weights = "weights",
subclass = "subclass",
data = NULL,
include.s.weights = TRUE,
drop.unmatched = TRUE)

get_matches(object,
distance = "distance",
weights = "weights",
subclass = "subclass",
id = "id",
data = NULL,
include.s.weights = TRUE)

## Arguments

object a matchit object; the output of a call to matchit(). which group should comprise the matched dataset: "all" for all units, "treated" for just treated units, or "control" for just control units. Default is "all". a string containing the name that should be given to the variable containing the distance measure in the data frame output. Default is "distance", but "prop.score" or similar might be a good alternative if propensity scores were used in matching. Ignored if a distance measure was not supplied or estimated in the call to matchit(). a string containing the name that should be given to the variable containing the matching weights in the data frame output. Default is "weights". a string containing the name that should be given to the variable containing the subclasses or matched pair membership in the data frame output. Default is "subclass". a string containing the name that should be given to the variable containing the unit IDs in the data frame output. Default is "id". Only used with get_matches(); for match.data(), the units IDs are stored in the row names of the returned data frame. a data frame containing the original dataset to which the computed output variables (distance, weights, and/or sublcass) should be appended. If empty, match.data() will attempt to find the dataset using the environment of the matchit object, which can make this unreliable if match.data() is used in a fresh R session or environment different from the original calling environment (e.g., inside a function) or if the original dataset changed between calling matchit() and match.data(). It is always safest to supply a data frame, which should have as many rows as and be in the same order as the data in the original call to matchit(). The same goes for get_matches(), which calls match.data() internally. logical; whether to multiply the estimated weights by the sampling weights supplied to matchit(), if any. Default is TRUE. If FALSE, the weights in the match.data() or get_matches() output should be multiplied by the sampling weights before being supplied to the function estimating the treatment effect in the matched data. logical; whether the returned data frame should contain all units (FALSE) or only units that were matched (i.e., have a matching weight greater than zero) (TRUE). Default is TRUE to drop unmatched units.

## Details

match.data() creates a dataset with one row per unit. It will be identical to the dataset supplied except that several new columns will be added containing information related to the matching. When drop.unmatched = TRUE, the default, units with weights of zero, which are those units that were discarded by common support or the caliper or were simply not matched, will be dropped from the dataset, leaving only the subset of matched units. The idea is for the output of match.data() to be used as the dataset input in calls to glm() or similar to estimate treatment effects in the matched sample. It is important to include the weights in the estimation of the effect and its standard error. The subclass column, when created, contains par or subclass membership and should be used to estimate the effect and its standard error. Subclasses will only be included if there is a subclass component in the matchit object, which does not occur with matching with replacement, in which case get_matches() should be used. See vignette("estimating-effects") for information on how to use match.data() output to estimate effects.

get_matches() is similar to match.data(); the primary difference occurs when matching is performed with replacement, i.e., when units do not belong to a single matched pair. In this case, the output of get_matches() will be a dataset that contains one row per unit for each pair they are a part of. For example, if matching was performed with replacement and a control unit was matched to two treated units, that control unit will have two rows in the output dataset, one for each pair it is a part of. Weights are computed for each row, and are equal to the inverse of the number of control units in each control unit's subclass. Unmatched units are dropped. An additional column with unit IDs will be created (named using the id argument) to identify when the same unit is present in multiple rows. This dataset structure allows for the inclusion of both subclass membership and repeated use of units, unlike the output of match.data(), which lacks subclass membership when matching is done with replacement. A match.matrix component of the matchit object must be present to use get_matches(); in some forms of matching, it is absent, in which case match.data() should be used instead. See vignette("estimating-effects") for information on how to use get_matches() output to estimate effects after matching with replacement.

## Value

A data frame containing the data supplied in the data argument or in the original call to matchit() with the computed output variables appended as additional columns, named according the arguments above. For match.data(), the group and drop.unmatched arguments control whether only subsets of the data are returned. See Details above for how match.data() and get_matches() differ. Note that get_matches sorts the data by subclass and treatment status, unlike match.data(), which uses the order of the data.

The returned data frame will contain the variables in the original data set or dataset supplied to data, and the following columns:

distance

The propensity score, if estimated or supplied to the distance argument in matchit() as a vector.

weights

The computed matching weights. These must be used in effect estimation to correctly incorporate the matching.

subclass

Matching strata membership. Units with the same value are in the same stratum.

id

The ID of each unit, corresponding to the row names in the original data or dataset supplied to data. Only included in get_matches output. This column can be used to identify which rows belong to the same unit since the same unit may appear multiple times if reused in matching with replacement.

These columns will take on the name supplied to the corresponding arguments in the call to match.data() or get_matches(). See Examples for an example of rename the distance column to "prop.score". If data or the original dataset supplied to matchit() was a data.table or tbl, the match.data() output will have the same class, but the get_matches() output will always be a base R data.frame. In addition to their base class (e.g., data.frame or tbl), returned objects have the class matchdata or getmatches. This class is important when using rbind() to append matched datasets.

matchit()

rbind.matchdata()

vignette("estimating-effects") for uses of match.data() and get_matches() in estimating treatment effects.

## Examples

data("lalonde")

# 4:1 matching w/replacement
m.out1 <- matchit(treat ~ age + educ + married +
race + nodegree + re74 + re75,
data = lalonde, replace = TRUE,
caliper = .05, ratio = 4)

m.data1 <- match.data(m.out1, data = lalonde,
distance = "prop.score")
dim(m.data1) #one row per matched unit
#>       treat age educ   race married nodegree re74 re75       re78 prop.score
#> NSW1      1  37   11  black       1        1    0    0  9930.0460 0.63876993
#> NSW2      1  22    9 hispan       0        1    0    0  3595.8940 0.22463424
#> NSW3      1  30   12  black       0        0    0    0 24909.4500 0.67824388
#> NSW4      1  27   11  black       0        1    0    0  7506.1460 0.77632408
#> NSW5      1  33    8  black       0        1    0    0   289.7899 0.70163874
#> NSW6      1  22    9  black       0        1    0    0  4056.4940 0.69906990
#> NSW7      1  23   12  black       0        0    0    0     0.0000 0.65368426
#> NSW8      1  32   11  black       0        1    0    0  8472.1580 0.78972311
#> NSW9      1  22   16  black       0        0    0    0  2164.0220 0.77983825
#> NSW10     1  33   12  white       1        0    0    0 12418.0700 0.04292461
#>       weights
#> NSW1        1
#> NSW2        1
#> NSW3        1
#> NSW4        1
#> NSW5        1
#> NSW6        1
#> NSW7        1
#> NSW8        1
#> NSW9        1
#> NSW10       1
g.matches1 <- get_matches(m.out1, data = lalonde,
distance = "prop.score")
dim(g.matches1) #multiple rows per matched unit
#>         id subclass weights treat age educ   race married nodegree        re74
#> 1     NSW1        1    1.00     1  37   11  black       1        1     0.00000
#> 2   PSID69        1    0.25     0  30   17  black       0        0 17827.37000
#> 3  PSID387        1    0.25     0  55    4  black       0        1     0.00000
#> 4  PSID386        1    0.25     0  20   12  black       0        0     0.00000
#> 5  PSID373        1    0.25     0  20   12  black       0        0     0.00000
#> 6     NSW2        2    1.00     1  22    9 hispan       0        1     0.00000
#> 7  PSID111        2    0.25     0  51   11  white       0        1    48.98167
#> 8   PSID66        2    0.25     0  26    8 hispan       0        1  3168.13400
#> 9  PSID339        2    0.25     0  26    9 hispan       0        1  1563.49500
#> 10 PSID150        2    0.25     0  22   11 hispan       0        1  7341.37300
#>        re75      re78 prop.score
#> 1     0.000  9930.046  0.6387699
#> 2  5546.419 14421.130  0.6385542
#> 3     0.000     0.000  0.6357539
#> 4     0.000 11594.240  0.6428929
#> 5     0.000     0.000  0.6428929
#> 6     0.000  3595.894  0.2246342
#> 7  3813.387  1525.014  0.2240794
#> 8  5872.258 11136.150  0.2225951
#> 9     0.000  2862.356  0.2161957
#> 10 2535.097 14187.650  0.2128728

format: plain


text:
# Method to solving this proof with a java app

I'm writing a program to solve this proof, but I don't know how to go about solving it. If anyone has some insight it would be great help. Thanks

• For every odd integer $n$, $3 \leq n \leq 199$, there exists an integer $m \geq 0$, and a prime number $p$, such that $n = 2^{m} + p$.
-
Recursively generate powers of two less than $n$, and check if the difference between $n$ and $2^m$ is prime using an appropriate method? – J. M. Nov 18 '10 at 23:56
You mean a search-routine that's simpler than brute force? – Unreasonable Sin Nov 18 '10 at 23:58
Ingenious. That's basically like solving for p, then checking if it's prime, right? – Trevor Arjeski Nov 19 '10 at 0:00
I planned on doing brute force. – Trevor Arjeski Nov 19 '10 at 0:01
As there are only 46 primes (one of which you don't need) less than 199 and only 7 powers of 2, you could just make the two lists, list the 7*45 sums, sort them, and check. – Ross Millikan Nov 19 '10 at 0:09

I'm answering so this is so it gets off the "Unanswered Questions" list... :)

## Algorithm

Here's a "brute force" (but somewhat fast) approach for a generalized upper limit of $N$. (This case is $N=199$

1. Generate (or obtain) a sorted array of all primes less than $N$. (This is $\pi(N)$, where $\pi$ is the prime counting function)
2. Generate a sorted array of all powers of two less than $N$. (There are $\lfloor\log_2(N)\rfloor$ of these.) This could be done dynamically (not stored in array) using bitshifts.
3. Iterate through all the odd numbers less than $N$ but greater than 3 (there are about $N/2$ of these).
1. For each odd number, iterate through the powers of two
2. Subtract the power of two from the odd number, and do a binary search on the prime array to see if the difference is prime.
3. If you find a prime, the statement has been verified for the odd number. Continue to next odd.

## Runtime:

Approximate $\pi(N) = \text{Li}(N)$, where $\text{Li}(x)$ is the logarithmic integral. Also assume we are given the prime list. Then, we have an upper bound of the runtime: $$\mathcal{O}\left(\frac{N}{2}\lfloor\log_2(N)\rfloor\cdot\log_2(\text{Li}(N))\right)\approx \mathcal{O}\left(N\lg(N)\cdot\lg(\text{Li}(N))\right)$$

(The first factor from iterating through odd numbers, the second from the number of powers of two, and the third factor from the binary search on the prime array.)

## Java Code

public class MathPost {
static final int N = 199;

//Assumed to have all primes under N
static int[] primeArray = {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199};

public static void main(String[] args) {
boolean isTrue = true;
boolean trueForAPowerOfTwo = false;

int powerOfTwo = 1;
int supposedPrime = 0;
int oddNumber;
for (oddNumber = 3; oddNumber <= N; oddNumber += 2) {
//We don't know if the proposition is true for this odd number with any power of two yet
trueForAPowerOfTwo = false;

for (int power = 0; (1 << power) < oddNumber; power++) {
//Find the power of two we are considering
powerOfTwo = 1 << power;

//Determine the difference of the power of two and the odd number.
supposedPrime = oddNumber - powerOfTwo;

//The following sets the "wasTrueForAPowerOfTwo" variable to true if we find a prime difference
trueForAPowerOfTwo |= java.util.Arrays.binarySearch(primeArray, supposedPrime) >= 0;
}

isTrue &= trueForAPowerOfTwo;

if (!isTrue) break;
}

System.out.printf("The proposition is %s\n", isTrue);
if (!isTrue) {
System.out.printf("Counterexample: The odd number %d cannot be expressed as desired.\n", oddNumber);
}

}
}

## The Result

The statement is false. Take the odd (prime) number $127$. It cannot be expressed as the sum of a power of two and another prime: \begin{align} 127-1 &= 126 &\text{ (Not prime)}\\ 127-2 &= 125 &\text{ (Not prime)}\\ 127-4 &= 123 &\text{ (Not prime)}\\ 127-8 &= 119 &\text{ (Not prime)}\\ 127-16 &= 110 &\text{ (Not prime)}\\ 127-32 &= 95 &\text{ (Not prime)}\\ 127-64 &= 63 &\text{ (Not prime)} \end{align}

-
I had solved this on my own 3 years ago, but I really appreciate the time you took to write this answer. – Trevor Arjeski Feb 28 '13 at 15:04
@TrevorArjeski wow. I'm impressed you actually came back to look at this! :) I figured you'd have solved it (3yrs is a while to work on a problem), but thought it was an interesting question and enjoyed writing up a solution... – apnorton Feb 28 '13 at 22:22


text:
# The Groove of Math

## Random Exercise #5: The homotopy fibre of CP^1 ↪ CP^∞ and a Serre Spectral Sequence

[It is really annoying that the title cannot contain any tex code:( ]

Hi all, sorry for not having written much in the blog lately, i was kind of busy. First of all, I’m happy to announce that next I’ll be a graduate student at UT in Austin. I’m so happy and I’m looking forward to start working there. Ok, back to business, on today’s random exercise I will deal with a problem I encountered today. Since I found it funny and interesting I will write it down here. We are asked to study the homotopy fibre of the map

$\mathbb{C}P^1 \hookrightarrow \mathbb{C}P^{\infty}$.

## Random Exercise #4: The nonexistent 5-Manifold in Bredon’s book

Dear all, sorry for not having posted anything so far but my master thesis (stable classification of certain families of four manifolds) and PhD interviews have been keeping me very busy 🙂

This is a quick exercise which I was asked to solve by a friend of mine and since I found it funny I decided to spell some more details about it in order to show the importance of the Bockstein exact sequence. This could be a typical exam exercise in “Topology 2” here in Bonn, so dear fellows students, read and comment if anything is unclear!

The following is exercise 7, page 366 from Bredon’s Topology and Geometry:

Show that there can be no $5$-manifold $M$ with

## Filling the Details #2: Representing 2-homology classes of a 4-manifold.

For the second article under the FtD tag, I’ve chosen a nice proposition which bothered me some times ago. It’s intuitively easy to believe but I found a scarce amount of details spelled out in the proofs one can find on the classical references. So here we have the following result:

Proposition: Every $2$-homology class of a smooth 4-manifold is generated by a surface


text:
LaTeX is a high-quality typesetting system; features are designed for the production of technical and scientific documentation. It’s the de-facto standard for the communication and publication of scientific documents, and available as free software. LaTeX is actually built on the TeX typesetting system created by the legendary Donald Knuth.

LaTeX is nothing more than a series of TeX macros, providing ready made commands for common formatting and layout needs, such as section headings, footnotes, bibliographies and cross references.

The name LaTeX is pronounced lah-teck or lay-teck, and NOT lay-tecks. The characters in TeX, on which LaTeX is based, on the greek word τέχνη (skill, art, technique). Creator Donald Knuth promotes a pronunciation of /tɛx/ (tekh) (using a voiceless velar fricative as in Modern Greek, similar to the ch in loch).

# Install base Tex system

Tex Live is one great (GPLv2 licenced) way to get up and running with a completely standalone TeX system. In addition to the core TeX compiler, can bundle in optional packages and macros. These really are essential if you want to be productive.

On Arch, you can installed the everything with sudo pacman -S textlive-most. For Debian based distros sudo apt install textlive-full

Alternatively cherry pick the individual packages (e.g. I went with these on my Arch system):

• texlive-core base system and libs
• texlive-bin CLI programs and compiler
• texlive-latexextra LaTeX macros
• texlive-publishers
• texlive-fontsextra fonts
• pdflatex for rendering LaTeX to PDF

# Pandoc

Pandoc is an amazing program (written in haskell) that converts document formats from one format to another. Some examples are Markdown to reveal.js. HTML to docx. In the specific case of writing technical reports, I’m most interesting in going from Markdown to LaTeX.

On Arch, install these packages:

• pandoc base Pandoc system
• pandoc-citeproc bibliography support
• pandoc-crossref cross referencing of figures and tables

# Author paper

Write the paper in plain old markdown. Use YAML frontmatter to specify variables for Pandoc to use in its LaTeX template. By default Pandoc will produce its own template pandoc -D latex. Supported front matter options for a LaTeX targetted output format are documented here (along with others). Here’s an example from one of my papers:

---
title: "Semantic Similarity in Binaries with BinHunt"
date: "May 2019"
author: "Benjamin Simmonds, UNSW Canberra"
abstract: Extracting meaningful semantic differences between software binaries without source code is difficult. This is a challenging problem due to the overwhelming amount of syntactic noise that small changes can result in at the assembly level.
---

# Introduction

Software binaries, or architecture specific compiled programs are often

# Create bibliography (BibTeX)

If not familiar with BibTeX, don’t worry the basics are easy.

BibTeX is reference management software for formatting lists of references. The BibTeX tool often coupled with the LaTeX document preparation system.

Create a BibTeX file with all the references you consumed to produce your research. I for example always create a file called paper.bib, and cite every reference I use in this file. An excerpt from one:

@misc{tardos2006algorithm,
title={Algorithm Design},
author={Tardos, Eva and Kleinberg, Jon},
year={2006},
}

@ONLINE {symbolicexecution2019,
author = "Schroeder, Brian and Burget, Joel",
title  = "A gentle introduction to symbolic execution",
month  = "apr",
year   = "2019",
url    = "https://blog.monic.co/a-gentle-introduction-to-symbolic-execution/"
}

@article{eskandari2012graph,
title={A graph mining approach for detecting unknown malwares},
author={Eskandari, Mojtaba and Hashemi, Sattar},
journal={Journal of Visual Languages \& Computing},
volume={23},
number={3},
pages={154--162},
year={2012},
publisher={Elsevier}
}

@inproceedings{jin2012binary,
title={Binary function clustering using semantic hashes},
author={Jin, Wesley and Chaki, Sagar and Cohen, Cory and Gurfinkel, Arie and Havrilla, Jeffrey and Hines, Charles and Narasimhan, Priya},
booktitle={2012 11th International Conference on Machine Learning and Applications},
volume={1},
pages={386--391},
year={2012},
organization={IEEE}
}

Don’t worry, most academic and journal databases support BibTeX formatting out of the box. Using Google Scholar for example locate the relevant paper you want to cite, click the Cite option and then the BibTeX link.

The pandoc-citeproc filter automatically creates a references section at the end of the document, and will substitute all references in the markdown source in academic citation style. To reference a BibTeX entry in the content of the paper, simply prefix the BibTeX tag with an ampersand @, for example:

The novel concept first introduced (@jin2012binary, P42) blah blah blah

References with BibTeX entry:

@inproceedings{jin2012binary,
title={Binary function clustering using semantic hashes},
author={Jin, Wesley and Chaki, Sagar and Cohen, Cory and Gurfinkel, Arie and Havrilla, Jeffrey and Hines, Charles and Narasimhan, Priya},
booktitle={2012 11th International Conference on Machine Learning and Applications},
volume={1},
pages={386--391},
year={2012},
organization={IEEE}
}

# Render the paper as PDF

You could run the following pandoc command on the shell. I prefer to use make, which will happily take on the tedious task of building the PDF version of your paper. Simply run make.

To do this create a new file called Makefile, placing it in the same directory as paper.md (the Markdown file to be converted into LaTeX by pandoc) and paper.bib (a vanilla BibTeX file) files:

all: clean paper

paper:
pandoc -s -F pandoc-crossref -F pandoc-citeproc --bibliography=paper.bib \
--variable papersize=a4paper \
--variable classoption=twocolumn \
-s paper.md -o paper.pdf

tex:
pandoc -s -F pandoc-crossref -F pandoc-citeproc --bibliography=paper.bib \
--variable papersize=a4paper \
-s paper.md -t latex -o paper.tex

clean:
rm -f *.pdf *.log *.tex

By default, running make will kick off the all target, first cleaning up old output files, following up by creating a fresh PDF output. Optionally if you want the raw TeX output, can do so with a make tex.

# Use Git

Authoring documents using Markdown and LaTeX is absolutely life changing. You end up with the most beautifully typeset document output known to man, while maintaining the source for the report in a plain text format. Plain text (unlike nasty proprietary binary formats like docx) lends itself to:

• Diffing (with diff) and patching (with patch)
• Keeping in a version control system (VCS) like Git, to track the when, who, what of all changes.
• Can use your favourite text editor (i.e. Vim)

I personally have a single Git repo for all my papers. In the top level of the repo, create a .gitignore based on the useful TeX.gitignore template from github/gitignore, to prevent generated or temporary files for being commmitted into the repo.

Periodically push the repo to GitHub as a convenient backup and distribution mechanism.

Gist by macogden

TeX

LaTeX


text:
Browse Questions

# According to the properties of lograthmic functions, what is the value of: $\log_a \large \frac{x}{y}$?

Answer: $\log_a \large \frac{x}{y}$$= \log_a x - \log_a b$


text:
# Insert a reference in a caption?

After visiting wiki on labels and cross-referencing, a couple of tex sites (1,2), and a stackexchange solution, the following reference will not appear in the caption. I'm using pdflatex on ubuntu.

The section label that is being referenced (appen) is defined at the end of the doc. All of the following attempts to insert the reference in the caption produce ?? or [?]:

\begin{figure}[!ht]
\caption*{See the ~\ref{appen} for language abbreviations.}
\makebox[\textwidth]{\includegraphics[width=100mm]{{/filepath/hourly.update}.pdf}}
\end{figure}

\begin{figure}[!ht]
\caption*{See the ~\cite{appen} for language abbreviations.}
\makebox[\textwidth]{\includegraphics[width=100mm]{{/filepath/hourly.update}.pdf}}
\end{figure}

\begin{figure}[!ht]
\caption*{See the ~\protect\ref{appen} for language abbreviations.}
\makebox[\textwidth]{\includegraphics[width=100mm]{{/filepath/hourly.update}.pdf}}
\end{figure}

\begin{figure}[!ht]
\caption*{See the ~\protect\cite{appen} for language abbreviations.}
\makebox[\textwidth]{\includegraphics[width=100mm]{{/filepath/hourly.update}.pdf}}
\end{figure}

\section{Appendix}
\label{appen}
\subsection{Language Abbreviations}
...words, figures, tables, and such.

-- edit 2013-02-12

LaTeX Warning: Reference appen' on page 2 undefined on input line 94

Any idea how to insert a reference within a caption?

• As far as I can tell, this question does not share a relationship with the supposed duplicate. The error that I'm receiving is as follows: -- LaTeX Warning: Reference appen' on page 2 undefined on input line 94 – blehman Feb 13 '14 at 6:38
• Did you read Question mark instead of citation number completely? – Reinstate Monica - M. Schröder Feb 13 '14 at 7:24
• Yes. BibTeX is not the problem nor is that package being used. The problem is occurring with a reference within a caption. Your link is a great explaination for another problem dealing with bibliographies. As you can see both in my question and in the answer, I'm labeling a section so ~\cite{} is not useful. Thanks for the help. – blehman Feb 13 '14 at 20:32

The first one (using \ref{appen}) will give you, after two compilations, the expected result. The ones using \cite (which is to be used for bibliographical items) are incorrect. The one using \protect\ref{appen} will also yield the correct result, but \protect is not necessary. An example (notice that leaving a blank space before ~ will produce a superfluous space):

\documentclass{article}
\usepackage[demo]{graphicx}
\usepackage{caption}

\begin{document}

\begin{figure}[!ht]
\caption*{See the Section~\ref{appen} for language abbreviations.}
\makebox[\textwidth]{\includegraphics[width=100mm]{{/filepath/hourly.update}.pdf}}
\end{figure}

\section{Appendix}
\label{appen}

\end{document}

The general mechanism for cross-referencing is to place

\label{<key>}

after an anchor has been generated (typically, after some automatic numbering has been done), and to use

\ref{<key>}

to cross-reference the labelled object.

The demo option for graphicx simply replaces actual figures with black rectangles; do not use that option in your actual document.

• I really like that you warned about the demo option. I remember carrying over a lot of lipsums and similar example producing packages and options from answers on TeX.SX when I first started using the site. – Dennis Feb 13 '14 at 2:26
• @Gonzalo, those suggestions were tried and still produce the same error -- LaTeX Warning: Reference `appen' on page 2 undefined on input line 94. -- Let me know if you have any other ideas. Thanks. – blehman Feb 13 '14 at 6:32
• @user28435 Did you try compiling my example code exactly as it is? What was the result after two compilations? – Gonzalo Medina Feb 13 '14 at 12:13


text:
## Monday, May 30, 2011

### WSOP Utility Analysis revisited, part 1: How much of an edge is necessary to overcome tax and utility effects in the WSOP Main Event?

A few months ago, I looked at the true costs of playing the WSOP Main Event for the typical player after adjusting for the effects of tax and risk aversion. It's worth a read if you missed it. Overall, we found that there are significant tax and utility effects that a player must overcome, and thus that there is significant extra effective rake in playing the event.

Exactly how much of a skill advantage does a prospective WSOP Main Event player need in order to overcome this effective rake? Can we... quantify it?

The old model was built on the assumption that the average player was equally likely to finish in any of the 7,319 places in the tournament. To allow for our player in question to have a level of skill different from that of the rest of the field, we need a way to map player skill to a specific finish probability distribution, a specification for the probabilities of finishing in each of the 7,319 different places.

Connecting skill edge to a finish probability distribution

If we generalize our scope to a tournament with N players, then we know that, when all players are equally skilled, each player will finish in each place with probability 1/N, as we've stuck to in the old model. If a specific player has a different strategy than the other players (whom we will assume to all still be uniformly skilled), then what does his finish probability distribution look like?

We quantify skill through the typical tournament results convention of raw ROI, which we define as the player's pure return on investment in the absence of tax and utility considerations, but including the rake. For example, in the WSOP main event, since $600 of the$10,000 entry fee goes to rake, the average player's raw ROI is -6%. This is ROI as it is commonly used in discussions and in tracking software; we add "raw" to convey that it is in pure dollar terms, without utility or tax effects. Since this is the measure of skill edge that is commonly used, it is desirable for us to define our correspondence between player strategy and finish probability distribution by creating a mapping between raw ROI and a probability distribution.

My first naïve method of constructing such a family of finish probability distributions was to start from the uniform 1/N distribution, but to perturb the terms such that the probability of finishing in each place was equal to a uniform constant plus the probability of finishing in the prior place. In math, this family of distributions looks like:

In addition to the necessary properties of probability distributions (nonnegative probabilities which sum to 1), this naïve distribution has a basic property that would likely be desirable in any probability distribution for tournament finishes: the distribution is increasing in that, when the player's skill edge is positive, he has a higher probability of finishing in each place than he does the next-worst place. (Similarly, when the player's skill edge is negative, he has a higher probability of finishing in each place than he does the next-higher place.)

However, there is no reason that the increasing increments between each place should be equal throughout the entire distribution, as they are here. Our simple guess may not approximate reality very well. There's also another problem with this family of distributions: it turns out that we can't capture every possible ROI (from -100% to a certain 1st place) through this type of distribution.

Note that this model, as well as all of the forthcoming models, produces the desired result of uniform probabilities of 1/N when we set the player's skill to zero.

A second naïve guess for developing a family of finish probability distributions is to reframe the tournament as a series of heads-up matches, and to give our player a specific probability for winning each heads-up match. This generalizes through logarithms to tournament field sizes that are not powers of 2. In this approach, we get a different probability distribution depending on whether or not we start from the probability of getting 1st place and use conditional probabilities downward through the places, or whether we instead start from the bottom and "condition upwards". If we condition downwards, this distribution looks like:

This is again an increasing distribution, and one that should make sense for a bona fide heads-up tournament with 2^k players, but other than that, there is no intuition as to how it might apply to a non-heads-up tournament.

The most elegant and accurate approach that I've come across is one based on the following foundation: we assume that a single player who is superior to his uniform opponents and who is playing in a tournament with equal starting stacks (as is the case in real poker tournaments) might have a similar finish probability distribution to a player who is not superior to his opponents, but instead starts the tournament with a larger chip stack than the rest of the field. (Thanks to Aaron Brown for this idea for modeling finish probability distributions, and if you enjoy both poker and finance and have not read his excellent book The Poker Face of Wall Street, you are missing out.)

To reiterate:

A superior player's skill edge in a tournament could be approximated in a model that removes his strategic advantage but instead gives him a larger starting chip stack.

From this point, we just need a way of mapping one's chip stack to one's finish probability distribution. The most popular method of this is the Independent Chip Model, i.e. the Malmuth-Harville tournament chip valuation algorithm, and that's what I've chosen to use, as it happens to produce an identical result to instead parameterizing the player's skill edge by starting with his probability of finishing in 1st place and then conditioning downwards. We don't get a closed-form formula, but we get an iterative formula that is easy to implement in Excel:

See The Mathematics of Poker for a discussion of other methods of tournament chip stack valuation; it seems that Malmuth-Weitzman is the only reasonable competing theory to ICM out there, and it only differs from ICM in how a busted player's chips are distributed among the remaining players on average. I haven't spent too much time thinking about the differences between these models. My understanding is that each of them diverges from real historical tournament results, though I do not know if anybody has analyzed this rigorously.

Adding this finish probability distribution to the WSOP expected utility analysis

The hard part's done! Now we can just take this parameterized finish probability distribution and plug it into our good old WSOP Main Event expected utility analysis and see what we get.

We see that, under this model, the relationship between the CE payoff and the player's raw ROI is close to linear. Examining the data confirms that this visual intuition is accurate.

In particular, we see that, for our typical player ($80k net worth,$50k income with at least $10k from poker, risk aversion of 0.8), the ROI needed to simply break even after tax and utility considerations is 86%. So that answers the question we were left with when we first looked at expected utility analysis for the WSOP Main Event. After correcting for tax effects and risk aversion, the typical player needs to have a raw ROI of 86% in order to break even in certainty equivalent by playing the WSOP Main Event. Having higher net worth helps alleviate this high threshold. A player with a net worth of$500,000 and a YTD salary of $100,000 needs a 42% ROI to break even. A player with a net worth of$5,000,000 and a YTD salary of $1,000,000 needs only a 7% ROI. Again, we see that the WSOP Main Event is best suited towards the wealthy, despite its "everyman" appeal... though it is precisely the fact that the event appeals to so many everymen that causes the risk to be too great for them. Being less naturally risk-averse also brings down this threshold. The typical player with$80,000 net worth and a YTD salary of $50,000 needs only a 61% ROI if we reduce his ρ from 0.8 to 0.5. Further reducing ρ to 0.2 still leaves him needing a 33% ROI. Reducing it all the way to 0 (i.e. no risk-aversion, and utility is realized on the full amount of dollars remaining after taxes) brings it down to a 11% ROI. Risk aversion is demanding a higher minimum skill level than tax effects are. Coming up next... That's all for today, since building the finish probability distribution took so long. But now that the hard part's out of the way, we can play around with our expanded model in some other practical ways. Note that I realize that this year's WSOP Main Event is likely to have far fewer than 7,319 competitors. For a smaller tournament, both tax and utility effects will be lessened, and lower skill edge will be necessary to overcome these. It would be interesting to generalize this model to an arbitrary N-player tournament, if there were an algorithmic, standardized payout formula that covered all possible field sizes. (Does one exist for any popular tournament series? Drop me a comment if there is one.) For now, knowing that a player with typical risk aversion and tax effects will need to have a pure expectation of$8,600 in the WSOP Main Event is a pretty valuable baseline that can guide decisions. This again illustrates how the large field size and high buyin of the WSOP Main Event effectively crowds out the average player form participating — unless that player is willing to pay a high effective cost.

But what if the player could reduce the size of the buyin for this event?

Next time, I look at the possibility of the player being able to sell off some of his action at 1-to-1 (i.e. no markup), which effectively reduces the buyin from $10,000 to something more manageable for less wealthy players. Indeed, doing this can change an unprofitable opportunity into a profitable opportunity for the risk-averse tournament player. For different levels of raw ROI, we will solve for optimal percentages of one's action to take in the idealized WSOP Main Event. This should be valuable information not only for prospective tournament players, but also those that might buy shares in them. Stay tuned. ## Saturday, May 7, 2011 ### To play or not to play: Optimal game selection with risky operators There are few certainties in the post-Black Friday online poker world. One fact that we are all still sure of is that nothing has changed with respect to the fact that playing online poker is legal on the federal level in the U.S. If an American player finds a poker site willing to take his business, and if he can manage to get his money to and from that site, he is not violating any federal law by doing so. So while the remaining U.S.-facing online poker sites are unlikely to be anywhere near as large, reputable, or financially liquid as what players are used to, as long as the risk of losing one's funds is small enough, there may still be some (severely-reduced) profit opportunities for American poker players, or, at the very least, a means of continuing to stay competitive at poker by practicing one's skills at one's leisure. The remaining U.S. sites should be expected to be different than PokerStars and Full Tilt Poker in at least a few incredibly important ways: • Since the subset of former American players which chooses to move to these sites will be highly skewed towards serious, professional players, all players should expect their winrates to drop significantly. • Deposits and withdrawals will be more costly and more unstable, and players should account for some nonzero probability of never receiving a cashout. • In the event of either voluntary or government-induced site closure, due to liquidity issues and the fact that these sites are less reputable, there is some probability that U.S. players would never get their account balances returned to them. Not too comforting. To be sure, even if playing on the remaining smaller sites is otherwise acceptable and a reasonable substitute for the experience of pre-Black Friday online poker (which will not be the case for all poker players), these issues are serious and will annihilate a lot of utility. These added risks cause players to have to make tradeoffs involving their profit potential and the amount of money they risk in their account balances at these sites. Putting it all together, how big are these costs? How bad do these risks have to be before it's not worth even trying to continue to play online? We can quantify this with a model which can estimate the effects of these risks on expected utility through Monte Carlo simulation. Model inputs We start by making assumptions about various aspects of the risks and costs of playing on these sites. Some of these are deterministic, and others will be uncertain. Once we have settled upon reasonable estimates of these values for all of these parameters, we can vary the most critical parameters to find break-even thresholds for expected utility, which can guide player decisions. • Fees associated with deposits and withdrawals — Deposits and withdrawals at smaller sites are costly, in both inconvenience and fees, so managing account balances carefully will be important to avoid incurring too many of these expenses. Our model will consider fixed fees for depositing or withdrawing, which we can expand to include intangible costs that reflect the inconveniences of these money transfers. • Deposit/withdrawal strategies — The simplest way to describe a rule that guides when to cash out or deposit is a set of four numbers. When the player's account balance is below some critical level (such as when it is too low to play his chosen stakes and number of tables), he should redeposit to bring his account balance up to some higher threshold level. Similarly, since he doesn't want to keep an unnecessary amount of money in a risky account, when his account balance hits some upper critical level, he should withdraw to bring his account balance down to some lower threshold level (which might be close to, if not exactly the same as, the threshold level for deposits). • Winrate, standard deviation, and play volume — The model will approximate poker results by increments of the appropriate normal distribution. To match the rest of the model, rather than looking at winrate and standard deviation per hand, we can look at these on a per-day basis by factoring in the amount of hands the player plans to play each day. The number of days remaining in the year will be one of the inputs. • Per-day probability of site closure — In our model, at the end of each day, we will assume that there is a fixed probability that the site will close forever. Each day will be independent of the last. • Probability of getting paid if the site closes — When the site closes, there is a chance that all player balances will be lost. • Probability of getting paid on each cashout — For each cashout prior to closure, there is a chance that the player will never see his money (or never be able to cash out in the first place). We can assume that this probability is less than that of the probability of getting paid when the site closes. • Tax and utility functions That's a lot of inputs, and a lot of inputs that we can't necessarily get great estimates of. Since the non-standard parts of this model are those pertaining to the risks of having money with the online site, we can ignore the lesser (or at least more standard) uncertainty on personal game performance parameters and instead assume that we know the player's exact winrate. Parameter choices We'll again work with our typical poker player:$80k net worth, $40k non-poker income,$10k year-to-date poker winnings, and isoelsatic utility with ρ=0.8, paying both federal and NJ income taxes. For simplicity, we'll assume that any money lost due to site-specific risks is treated the same as a poker loss and directly deductible against poker winnings, though I believe this is not exactly true.

UPDATE 06/18/2011: The above assumption may only be valid for those who file their taxes as professional poker players. Amateur poker players likely cannot deduct these losses at all in most cases, and this has serious implications for the model. See Update/correction to risky site model: Casualty losses and tax effects.

I'm not entirely familiar with costs of deposit and withdrawal at all of the remaining U.S.-facing poker sites, but for Carbon Poker, it looks like roughly $17 in costs to make a deposit, and about$5 per withdrawal for a medium-volume player.

We'll look at two possible levels of stakes for our player: multitabling $0.25/$0.50 NL at an hourly winrate of $10 with a standard deviation of about$70, or multitabling $0.50/$1.00 NL with an hourly winrate of $15 and a standard deviation of about$140. His $0.50/$1 winrate is higher in dollars, but lower relative to the stakes, which introduces not only the usual higher relative risk to his results, but also a need to keep a higher balance with the site. Keep in mind that he will have a much tougher time achieving these winrates on a remaining U.S.-facing website, since the player pools will be significantly tougher and it may not always be possible to play as many tables as he would normally be accustomed to. In the post-Black Friday market, I would expect winrates to be significantly lower than the good old days, and most players should probably drop down at least a level or two.

We assume that, for $0.25/$0.50 NL, he will redeposit when his balance falls below $800, and he will redeposit up to$1,200. He will withdraw down to $1,200 when his balance hits$2,000. For $0.50/$1 NL, we will double each of these deposit/withdrawal thresholds. Once we settle upon the other parameters, we can play with this to see which withdrawal strategy is optimal.

As a default, based on no science at all and simply my gut assessment of the risks of the current poker market, we'll assume that the per-day probability of site closure is 0.002 (i.e. the site lasts, on average, 500 days), that there is a 50% probability of players getting paid if the site closes, and a 75% probability of getting paid on any given cashout prior to closure.

We'll have him start playing in early May, so there's about N=240 days left in the year.

Algorithm

To calculate the expected utility of playing with these risks, these winrate parameters, and this utility function, we will simulate the system through the following steps:
1. Initialize the player's starting bankroll by him making his first deposit, tracking the costs of doing so in a running net total starting from his year-to-date winnings. This running net will NOT include his day-to-day poker results; money on the site is not counted as a gain or loss until it is successfully withdrawn.
2. At the start of each day, the player plays his daily poker session, and we adjust his site balance by a normal random variable with appropriate mean and variance. Since he can't lose more than he has on the site, we truncate this normal distribution on both tails by preventing the magnitude of the swing from exceeding his existing balance.
3. After his session, if his bankroll is below his deposit threshold, he deposits according to his strategy. We track both the deposit itself and the costs of doing so in the running net.
4. Then, if his bankroll is above his withdrawal threshold, he withdraws according to his strategy. He ends up receiving the cashout according to the chosen probability, which will be added to the running net, less the costs. If he gets unlucky and misses on the cashout, we assume it is lost forever.
5. At the end of each day, the site closes down according to the chosen per-day probability. If the site shuts down, the player's balance is returned to him according to the probability of getting paid if the site closes. Either way, if the site closes, we exit the loop and go directly to step 7, as there's no longer any poker to be played (we neglect the possibility of choosing another remaining site for the rest of the year).
6. Repeat steps 2-5 for each of the remaining N days.
7. At the end of the year, if the site is still open, for the purposes of evaluating his year-end utility, the player withdraws his balance (and receives it according to the chosen probability). This gives his final net winnings for the year, and we evaluate the after-tax utility of this amount.
8. Repeat steps 1-7 in a Monte Carlo simulation to simulate the average expected utility.

Some results

For the parameters we chose, we can look at how the results change as we perturb the most uncertain parameters, those related to site-specific risks.

First, if we disregard the assumption of getting paid 50% of the time if the site closes, and instead vary that, we can see what the effects of additional liquidity are and find the break-even liquidity level. Since this liquidity and the probability of successful cashouts are linked, it is desirable to adjust both at the same time, so as we vary the liquidity level, we will assume that the risk of losing a cashout is always half the risk of losing a bankroll in the event of closure.

The rightmost point is when the site is fully secure, which is what decision-makers have been roughly used to in the pre-Black Friday environment. In this case, the player prefers to play his more pure-EV-profitable stake of 100NL. As the site gets less and less secure, the need to keep a higher bankroll at the site for 100NL will create more and more risk, and 50NL becomes better when the site liquidity falls below about 38%. For less than 5% liquidity, the increase in the player's utility by playing either stake is negative, so the player should refrain from playing at all, unless the entertainment value or the value of being able to practice and work on his game is worth the cost.

If we return the liquidity probability to its original default rate of 50%, we can instead look at how expected utility varies with the per-day probability of site closure.

When changing the average lifetime of the poker site, we see that the effect is less linear. Again, the rightmost point corresponds to perfect safety, a site that has no risk of closing. As the probability of closing is increased (and we assume a 50% chance of not being able to get money back after closure), the expected utility drops off significantly. 50NL overtakes 100NL around a per-day closure probability of 0.004 (i.e. a mean lifetime of 250 days). Not playing at all becomes the best choice at a per-day closure probability of around 0.01 (i.e. a mean lifetime of 100 days).

If we return to the original, fixed best estimates of these parameters (50% liquidity, average site lifetime of 500 days), we can tweak the deposit and withdrawal strategies to see which one works best. We'll stick with just 50NL here, since it seems to be the better choice. We need to keep the lower deposit threshold at $800 in order for the player to have enough bankroll to be able to play a bunch of tables of 50NL during his sessions, but we can vary the other three account balance strategy parameters. We see that we can actually do better than the deposit and withdrawal strategy that we initially chose. It turns out to be better to wait to withdraw until at a higher balance, but still to deposit up to just$1,200. I suppose that this asymmetry comes about due to our utility function; we have a higher relative risk tolerance when we have earned more money on the year. We see that the "tightest" deposit and withdrawal thresholds do not do well here. This is all due to the interplay between the site-specific risks and the costs of moving money, and this behavior will probably change when any of these are modified.

Conclusions

If we compare the middle parts of each of the graphs to the rightmost points (which correspond to perfect safety), we can really see how much our utility suffers under these new risks. Players will need to have significant edges over their competition to overcome this.

In these results, the "shape" of the results is no surprise. We're not seeing anything we wouldn't have been able to guess without doing the math, but quantifying it provides a useful framework for guiding our play decisions under the uncertainty of the site-specific risks.

For our example player here, we see that, under the default assumptions with 50% liquidity and an average site lifetime of 500 days, playing 100NL will be a better risk-adjusted value than 50NL, though it's close enough that the player should move down if he feels the liquidity risk is greater, or that the sites have a significantly shorter expected lifetime. To dissuade him from playing entirely, there would have to be either a very low probability of getting paid back if the site closed, or withdrawals while the site is still operating would have to be less than 50% likely to go through. So, for this winning player, it looks like playing on a remaining U.S.-facing website should be better than quitting entirely as long as he has at least some faith in the current market. Many less skilled players with thinner winrates will be forced out of the market entirely from these additional risks.

This model is well-suited to being dynamically updated as the year goes on. The time horizon and year-to-date winnings will change constantly, which will shape bankroll decisions. Also, it should be useful to update the site-specific risks with better estimates as time goes on and as more information develops in the wake of Black Friday, or as the market positioning of the remaining U.S.-facing sites change.

While it's easy to see results by varying the site-specific parameters in this model when the player-specific variables (utility, income, winrate, etc.) are fixed, it's hard to draw broader conclusions over a more general player base. This is the sort of model that is best applied on a case-by-case basis, with each user's particular play variables fixed and known. If enough people are interested in using the spreadsheet I wrote for this, I might consider cleaning it up and hosting it somewhere.


text:
#### PHP Conditional Statements Publishing Model Question Paper 1

12th Standard

Reg.No. :
•
•
•
•
•
•

Computer Applications

Time : 01:00:00 Hrs
Total Marks : 50
5 x 1 = 5
1. What will be the output of the following PHP code?
< ? php
 x;
if (x)
print “hi” ;
else
print “how are u”;
? >

(a)

how are u

(b)

hi

(c)

error

(d)

no output

2. What will be the output of the following PHP code ?
< ?php
x;
if (x == 0)
print “hi” ;
else
print “how are u”;
print “hello”
? >

(a)

how are uhello

(b)

hihello

(c)

hi

(d)

no output

3. Statement which is used to make choice between two options and only option is to be performed is written as

(a)

if statement

(b)

if else statement

(c)

then else statement

(d)

else one statement

4. What will be the output of the following PHP code ?
< ?php
a = “”;
if (a)
print “all”;
if
else
print “some”;
? >

(a)

all

(b)

some

(c)

error

(d)

no output

5. What will be the output of the following PHP code ?
< ?php
if (-100)
print “hi” ;
else
print “how are u”;
? >

(a)

how are u

(b)

hi

(c)

error

(d)

no output

6. 5 x 2 = 10
7. Define Conditional Statements in PHP

8. Define if statement in PHP.

9. What is if else statement in PHP?

10. Write the Syntax of Switch statement.

11. Compare if and if else statement.

12. 5 x 3 = 15
13. Write the features Conditional Statements in PHP.

14. Differentiate Switch and if else statement

15. Differentiate if statement and if else if else stament.

16. Write a php program to check whether the given number is positive or negative.

17. Write, a PHP program to display independenc day and republic day the greetings using If elseif else statement.

18. 4 x 5 = 20
19. Explain the if elseif else statement. With suitable example.

20. Explain if else statement in PHP. With an example.

21. Write a PHP code to display the days of weak using switch statement

22. Write a PHP program to display week days using switch case statement.


text:
• Customer Solutions ▾
• Competitions
• Community ▾
with —

# U.S. Census Return Rate Challenge

Finished
Friday, August 31, 2012
Sunday, November 11, 2012
\$1,000 • 244 teams

# External Data (deadline for new data sources is passed)

» Next
Topic
 Rank 7th Posts 13 Joined 20 Dec '10 Email user maternaj wrote: David, back to the never ending story regarding the external datasets.. ;) APPROVED--http://www.census.gov/geo/www/2010census/tractrel/trftxt/us2010trf.txt (only area fields, not POP and HU fields). I was one of those who pointed out about the statement that POP and HU fields are from 2010 census even the field is called POP00, HU00. However, I have just done some search and there are also fields called POPPCT00 and HUPCT00 that look clean and relating really to 2000 data. On the other hand, the brackets clause from the approved datasets page I would read that these are forbidden as well. The deadline is coming really quickly so the definitive yes/no for this will be appreciated.   maternaj, The definition of POPPCT00 is "Calculated Percentage of the POP00 this record (POP10PT) contains (to 2 decimal points)" and the definition of HUPCT00 is "Calculated Percentage of the HU00 this record (HU10PT) contains (to 2 decimal points)". From this, it seems like they are derived from fields POP00 and HU00. In that case they should not be allowed. Can someone please confirm if my understanding is correct? #196 / Posted 6 months ago
 Rank 5th Posts 10 Thanks 3 Joined 7 Jul '11 Email user Shashi Godbole wrote: maternaj wrote: David, back to the never ending story regarding the external datasets.. ;) APPROVED--http://www.census.gov/geo/www/2010census/tractrel/trftxt/us2010trf.txt (only area fields, not POP and HU fields). I was one of those who pointed out about the statement that POP and HU fields are from 2010 census even the field is called POP00, HU00. However, I have just done some search and there are also fields called POPPCT00 and HUPCT00 that look clean and relating really to 2000 data. On the other hand, the brackets clause from the approved datasets page I would read that these are forbidden as well. The deadline is coming really quickly so the definitive yes/no for this will be appreciated.   maternaj, The definition of POPPCT00 is "Calculated Percentage of the POP00 this record (POP10PT) contains (to 2 decimal points)" and the definition of HUPCT00 is "Calculated Percentage of the HU00 this record (HU10PT) contains (to 2 decimal points)". From this, it seems like they are derived from fields POP00 and HU00. In that case they should not be allowed. Can someone please confirm if my understanding is correct? Shashi, hmm, this is strange, when I was looking at some sample data it looked like the columns (POPPCT00, HUPCT00) have nothing to do with *10PT columns. At least for the merged tracts (where I think the info is more useful). Have just now checked for split tracts and it looks that the actual 2010 data are really used.. Sorry for the confusion, I have not noticed the description before and was doing the research based on the actual data only... :( Anyway, due to the time constraints (came with this "theory" too late) we decided not to use these two fields anyway but I think it would be many people's opinion that the whole "external data thread" was a bit of a nightmare for this competition.. :) #197 / Posted 6 months ago
 Rank 12th Posts 65 Thanks 34 Joined 14 May '10 Email user Hi all! maternaj wrote: Anyway, due to the time constraints (came with this "theory" too late) we decided not to use these two fields anyway but I think it would be many people's opinion that the whole "external data thread" was a bit of a nightmare for this competition.. :)   Especially for people who know the bad English! All the best, Alex. #198 / Posted 6 months ago
 Rank 7th Posts 13 Joined 20 Dec '10 Email user maternaj wrote: Shashi Godbole wrote: maternaj wrote: David, back to the never ending story regarding the external datasets.. ;) APPROVED--http://www.census.gov/geo/www/2010census/tractrel/trftxt/us2010trf.txt (only area fields, not POP and HU fields). I was one of those who pointed out about the statement that POP and HU fields are from 2010 census even the field is called POP00, HU00. However, I have just done some search and there are also fields called POPPCT00 and HUPCT00 that look clean and relating really to 2000 data. On the other hand, the brackets clause from the approved datasets page I would read that these are forbidden as well. The deadline is coming really quickly so the definitive yes/no for this will be appreciated.   maternaj, The definition of POPPCT00 is "Calculated Percentage of the POP00 this record (POP10PT) contains (to 2 decimal points)" and the definition of HUPCT00 is "Calculated Percentage of the HU00 this record (HU10PT) contains (to 2 decimal points)". From this, it seems like they are derived from fields POP00 and HU00. In that case they should not be allowed. Can someone please confirm if my understanding is correct? Shashi, hmm, this is strange, when I was looking at some sample data it looked like the columns (POPPCT00, HUPCT00) have nothing to do with *10PT columns. At least for the merged tracts (where I think the info is more useful). Have just now checked for split tracts and it looks that the actual 2010 data are really used.. Sorry for the confusion, I have not noticed the description before and was doing the research based on the actual data only... :( Anyway, due to the time constraints (came with this "theory" too late) we decided not to use these two fields anyway but I think it would be many people's opinion that the whole "external data thread" was a bit of a nightmare for this competition.. :) Indeed it was. #199 / Posted 6 months ago
 DavidChudzicki Competition Admin Kaggle Admin Posts 424 Thanks 106 Joined 21 Nov '10 Email user I agree this has been a bit of a pain. My apologies... we'll try to learn from it for the future. #200 / Posted 6 months ago