ES_recomb<-function(G,child) {
  nParent=nrow(G$b);nChild=nrow(child$b);nVar=ncol(G$b)
  for (i_child in 1:nChild) {
    parentID=sample(1:nParent,2)
    coID=sample(c(0,1),nVar,replace=T)
    child$b[i_child,]=G$b[parentID[1],]
    child$b[i_child,which(coID==1)]=G$b[parentID[2],which(coID==1)]
    child$sigma[i_child,]=0.5*(G$sigma[parentID[1],]+G$sigma[parentID[2],])
  }
  return(child)
}
ES_mutate<-function(child,tau){
  nChild=nrow(child$b);nVar=ncol(child$b)
  child$sigma<-child$sigma*exp(matrix(rnorm(nChild*nVar)*tau,nrow=nChild))
  child$b=child$b+child$sigma*matrix(rnorm(nChild*nVar),nrow=nChild,ncol=nVar)
  return(child)
}
ES_er<-function(b,x,y){
  yhat<-x%*%b
  return(sum((y-yhat)^2))
}

x=matrix(rnorm(4*50,mean=10,sd=2),ncol=4);x=cbind(rep(1,50),x)
y=x%*%c(10,2,5,-3,-5)+rnorm(50,mean=0,sd=2);
    
G=list();child=list();
nGen=1000;nParent=30;nChild=60;tau=1;nVar=5
G$b=matrix(rnorm(nVar*nParent,mean=0,sd=1),ncol=nVar)
G$sigma=matrix(runif(nVar*nParent),ncol=nVar)+0.5
child$b=matrix(0,nrow=nChild,ncol=nVar)
child$sigma=matrix(0,nrow=nChild,ncol=nVar)
optHist=matrix(0,nrow=nGen,ncol=1)
for (i_gen in 1:nGen) {
  child<-ES_recomb(G,child)
  child<-ES_mutate(child,tau)
  fitG=apply(G$b,1,ES_er,x,y);fitC=apply(child$b,1,ES_er,x,y)
  fitT=c(fitG,fitC);fitMin=sort(fitT,index.return=T)
  tempB=rbind(G$b,child$b);tempS=rbind(G$sigma,child$sigma)
  G$b=tempB[fitMin$ix[1:nParent],]
  G$sigma=tempS[fitMin$ix[1:nParent],]
  optHist[i_gen]=fitMin$x[1]
}

>  head(G$b)
         [,1]     [,2]     [,3]      [,4]      [,5]
[1,] 8.597338 2.118279 5.081328 -2.982018 -5.062286
[2,] 8.597347 2.118279 5.081328 -2.982018 -5.062286
[3,] 8.597336 2.118280 5.081328 -2.982018 -5.062286
[4,] 8.597333 2.118280 5.081328 -2.982018 -5.062286
[5,] 8.597343 2.118280 5.081328 -2.982018 -5.062286
[6,] 8.597336 2.118280 5.081328 -2.982018 -5.062287
> summary(lm(y~x[,2:5]))

Call:
lm(formula = y ~ x[, 2:5])

Residuals:
    Min      1Q  Median      3Q     Max 
-4.1866 -1.7043  0.3254  1.3890  4.1851 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.5973     3.9641   2.169   0.0354 *  
x[, 2:5]1     2.1183     0.2235   9.476 2.73e-12 ***
x[, 2:5]2     5.0813     0.2077  24.464  < 2e-16 ***
x[, 2:5]3    -2.9820     0.1606 -18.566  < 2e-16 ***
x[, 2:5]4    -5.0623     0.2120 -23.874  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.342 on 45 degrees of freedom
Multiple R-squared:  0.9737,	Adjusted R-squared:  0.9713 
F-statistic: 415.7 on 4 and 45 DF,  p-value: < 2.2e-16