Data Simulation Of Heritability Estimation Project

## [1] "Document was last updated at 2021-11-26."

---

Table of contents

1 yTKy_yTy

Heritability estimation expression

\[h^2=\frac{y^TKy-y^Ty}{tr(K^2)-n}\]

---

Estimation of yTKy_yTy

import numpy as np
import pandas as pd
n=1000
m=10000
#B=10
for w in range(3):
    m=m*(w+1)
    frq=[np.random.uniform(0.2,0.8) for i in range(m)]

    dat=np.zeros((n,m))
    for i in range(n):
        dat[i,:]=np.random.binomial(2, frq,m)
  
    x=pd.DataFrame(dat).apply(lambda x: (x-x.mean())/x.std())
    k=x.values @ x.values.T/x.shape[1]
    beta=np.random.normal(0,0.1,m)
    vg=0.5
    y1=x@beta
    ve=np.var(y1)*(1-vg)/vg
    e=np.random.normal(0,np.sqrt(ve),n)
    y=[y1[i]+e[i] for i in range(len(e))]
    y=pd.DataFrame(y).apply(lambda x: (x-x.mean())/x.std())
    p1=y.values.T @ k @ y.values
    p2=y.values.T@y.values
    print('n:{},m:{},yTKy:{},yTy:{}'.format(n,m,p1[0,0],p2[0,0]))

## n:1000,m:10000,yTKy:1066.8029172737927,yTy:999.0000000000001
## n:1000,m:20000,yTKy:1024.783732848335,yTy:999.0000000000001
## n:1000,m:60000,yTKy:1020.2180197567307,yTy:999.0000000000001

2 LB_tr(K2)

Basis theory

\[E[z^TAz]=tr[A]\]

Explaination

leads to the following unbiased estimator of the trace of \(K_2\) given \(B\) random vectors, \(z_1\);…;\(z_B\), drawn independently from a distribution with zero mean and identity covariance matrix \(I_N\):

\[L_B=tr[K^2]=\frac{1}{B}\sum_{b}z^T_bKKz_b\\ =\frac{1}{B}\frac{1}{M^2}\sum_{b}z^T_bXX^TXX^Tz_b\\ =\frac{1}{B}\frac{1}{M^2}\sum_{b}||XX^Tz_b||2^2\] ————————————–

Estimation of LB_tr(K2)

import numpy as np
import pandas as pd

def sum_list(ls):
    su=0
    for i in ls:
        su=su+i**2
    return su

def estimate_LB(x,B):
  n=x.shape[0]
  m=x.shape[1]
  LB=[]
  for i in range(B):
      zb=np.random.randn(n)
      L=x@(x.T@zb)
      LB.append(sum_list(L)/m**2)
  return(np.mean(LB))

n=1000
m=10000
B=10
for w in range(3):
    m=m*(w+1)
    frq=[np.random.uniform(0.2,0.8) for i in range(m)]

    dat=np.zeros((n,m))
    for i in range(n):
        dat[i,:]=np.random.binomial(2, frq,m)
  
    x=pd.DataFrame(dat).apply(lambda x: (x-x.mean())/x.std())
    k=x.values @ x.values.T/x.shape[1]
    tril=[k[i,j] for i in range(n) for j in range(i)]
    me=1/np.var(tril)
    lb=estimate_LB(x,B)
    print('n:{},m:{},LB:{},tr(k2):{}'.format(n,m,lb, n+n**2/me))

## n:1000,m:10000,LB:1076.2366503734575,tr(k2):1099.6898862687221
## n:1000,m:20000,LB:1047.0020229660327,tr(k2):1049.670293740906
## n:1000,m:60000,LB:1004.4503395967738,tr(k2):1016.6185149904431

3 Simulate_matrix trace

Simulate_matrix trace

import numpy as np
import pandas as pd
n=1000
m=10000
for w in range(2):
    m=m*(w+1)
    frq=[np.random.uniform(0.2,0.8) for i in range(m)]

    dat=np.zeros((n,m))
    for i in range(n):
        dat[i,:]=np.random.binomial(2, frq,m)
  
    x=pd.DataFrame(dat).apply(lambda x: (x-x.mean())/x.std())
    k=x.values @ x.values.T/x.shape[1]
    tril=[k[i,j] for i in range(n) for j in range(i)]
    me=1/np.var(tril)
    me
    def simulate_trace(n,me):
        k2_trace=n+n**2/me
        k3_trace=n+2+3*n**2/me
        k4_trace=n+3+(6*n**2+4*n+4)/me+(2*n**3+2*n**2-n)/me**2
        return(k2_trace, k3_trace,k4_trace)

    def ture_trace(k):
        k2=k@k.T
        k3=k2@k
        k4=k3@k.T
        k2_trace=np.trace(k2)
        k3_trace=np.trace(k3)
        k4_trace=np.trace(k4)
        return(k2_trace, k3_trace,k4_trace)
    simulate=simulate_trace(n,me)
    true=ture_trace(k)
    print('n:{},m:{},w:{} ,simulate:{},true:{}'.format(n,m,w,simulate,true))

## 10052.44802312114
## n:1000,m:10000,w:0 ,simulate:(1099.4782562118153, 1300.4347686354456, 1620.079477077645),true:(1098.5072383097033, 1307.5404234656746, 1656.9730908345628)
## 20071.324081566778
## n:1000,m:20000,w:1 ,simulate:(1049.8223234270024, 1151.466970281007, 1307.1029190139045),true:(1048.8350283147008, 1150.965521816546, 1312.9347453162452)

simulate_matrix trace-times

import numpy as np
import pandas as pd
n=300
#m=2400
for w in range(15):
    m=n*w
    frq=[np.random.uniform(0.2,0.8) for i in range(m)]

    dat=np.zeros((n,m))
    for i in range(n):
        dat[i,:]=np.random.binomial(2, frq,m)
  
    x=pd.DataFrame(dat).apply(lambda x: (x-x.mean())/x.std())
    k=x.values @ x.values.T/x.shape[1]
    tril=[k[i,j] for i in range(n) for j in range(i)]
    me=1/np.var(tril)
    me
    def simulate_trace(n,me):
        k2_trace=n+n**2/me
        k3_trace=n+2+3*n**2/me
        k4_trace=n+3+(6*n**2+4*n+4)/me+(2*n**3+2*n**2-n)/me**2
        return(k2_trace, k3_trace,k4_trace)

    def ture_trace(k):
        k2=k@k.T
        k3=k2@k
        k4=k3@k.T
        k2_trace=np.trace(k2)
        k3_trace=np.trace(k3)
        k4_trace=np.trace(k4)
        return(k2_trace, k3_trace,k4_trace)
    simulate=simulate_trace(n,me)
    true=ture_trace(k)
    if abs(simulate[2]-true[2])/true[2]<0.01:
        print('w:{} ,simulate:{},true:{}'.format(w,simulate,true))
        break
    else:
        continue

## nan
## 302.39124166699753
## 602.7878960258561
## 912.6446477779571
## 1208.5548415457224
## 1512.1571125116254
## 1824.787391314342
## 2122.429990235825
## 2438.993325901222
## 2724.52113120342
## 3058.0184652164703
## 3364.437291320357
## 3685.613312081603
## w:12 ,simulate:(324.4192736402856, 375.25782092085683, 453.8308860052123),true:(323.4281900916129, 374.22017366640284, 457.4581183477777)
## 
## <string>:10: RuntimeWarning: invalid value encountered in true_divide

4 me

Subexample_me

import numpy as np
import pandas as pd

def sum_list(ls):
    su=0
    for i in ls:
        su=su+i**2
    return su

def estimate_LB(x,B):
  n=x.shape[0]
  m=x.shape[1]
  LB=[]
  for i in range(B):
      zb=np.random.randn(n)
      L=x@(x.T@zb)
      LB.append(sum_list(L)/m**2)
  return(np.mean(LB),np.var(LB)/B)

n=2000
m=2000
#m=n*13
frq=[np.random.uniform(0.2,0.8) for i in range(m)]
dat=np.zeros((n,m))
for i in range(n):
    dat[i,:]=np.random.binomial(2, frq,m)
for j in range(3):
   w=round(n/(j+1))
   dat2=dat[0:w,:]  
   x=pd.DataFrame(dat).apply(lambda x: (x-x.mean())/x.std())
   x2=pd.DataFrame(dat2).apply(lambda x: (x-x.mean())/x.std())
   k=x.values @ x.values.T/x.shape[1]
   k2=x2.values @ x2.values.T/x2.shape[1]
   tril=[k[i,j] for i in range(n) for j in range(i)]
   me=1/np.var(tril)
   tril2=[k2[i,j] for i in range(w) for j in range(i)]
   me2=1/np.var(tril2)
   print('w:{},n:{},m:{},me:{},me2:{}'.format(w,n,m,me,me2))

## w:2000,n:2000,m:2000,me:2000.748000427016,me2:2000.748000427016
## w:1000,n:2000,m:2000,me:2000.748000427016,me2:2007.8253664192634
## w:667,n:2000,m:2000,me:2000.748000427016,me2:2019.3233134649427

5 Simulation of tr(K4)

**Tr(k4)_estimate12-true**

import numpy as np
import pandas as pd


def sum_list(ls):
    su=0
    for i in ls:
        su=su+i**2
    return su

def estimate_LB(x,B):
  n=x.shape[0]
  m=x.shape[1]
  LB=[]
  for i in range(B):
      zb=np.random.randn(n)
      L=x@(x.T@zb)
      LB.append(sum_list(L)/m**2)
  return(np.mean(LB),np.var(LB)/B)


n=1000
m1=1000
B=5
#a=[]
for w in range(15):
    m=m1*(w+1)
    frq=[np.random.uniform(0.2,0.8) for i in range(m)]
    dat=np.zeros((n,m))
    for i in range(n):
        dat[i,:]=np.random.binomial(2, frq,m)  
    x=pd.DataFrame(dat).apply(lambda x: (x-x.mean())/x.std())
    k=x.values @ x.values.T/x.shape[1]
    tril=[k[i,j] for i in range(n) for j in range(i)]
    me=1/np.var(tril)
    lb=estimate_LB(x,B)
    k4_trace1=np.trace(k@k.T@k@k.T)
    k4_trace3=lb[1]*B/2
    k4_trace2=n+(3*n**3-3*n**2+2*n)/(n-1)**3+6*n**2/me+(4*n**2+4*n)/((n-1)*me)+(2*n**3+2*n**2-n)/me**2
    print('w:{},n:{},m:{},k4_true:{},k4_estimate1:{},k4_estimate2:{}'.format(w+1,n,m,k4_trace1,k4_trace2,k4_trace3))

## w:1,n:1000,m:1000,k4_true:13913.74502006159,k4_estimate1:8960.767102184476,k4_estimate2:6161.254679801518
## w:2,n:1000,m:2000,k4_true:5588.130278590672,k4_estimate1:4487.375718743659,k4_estimate2:6205.260960592643
## w:3,n:1000,m:3000,k4_true:3689.634377018763,k4_estimate1:3220.0822687420527,k4_estimate2:7350.240502690656
## w:4,n:1000,m:4000,k4_true:2883.707843183356,k4_estimate1:2623.5903152261403,k4_estimate2:646.4649551349632
## w:5,n:1000,m:5000,k4_true:2445.912959656189,k4_estimate1:2282.391072293752,k4_estimate2:3578.7791455794613
## w:6,n:1000,m:6000,k4_true:2167.3961452550116,k4_estimate1:2056.8126219782953,k4_estimate2:2482.734544187215
## w:7,n:1000,m:7000,k4_true:1977.6861237339424,k4_estimate1:1898.2440132196093,k4_estimate2:170.70438591902789
## w:8,n:1000,m:8000,k4_true:1839.7833202315314,k4_estimate1:1780.7568085382436,k4_estimate2:3337.133420797455
## w:9,n:1000,m:9000,k4_true:1737.4736445964877,k4_estimate1:1691.6650731251611,k4_estimate2:1235.636431701536
## w:10,n:1000,m:10000,k4_true:1657.8292294384673,k4_estimate1:1621.4045208040413,k4_estimate2:908.4630726209447
## w:11,n:1000,m:11000,k4_true:1592.5048054636677,k4_estimate1:1563.360867930893,k4_estimate2:1144.2743666206588
## w:12,n:1000,m:12000,k4_true:1537.8089073924493,k4_estimate1:1513.9630342212229,k4_estimate2:1958.482319581203
## w:13,n:1000,m:13000,k4_true:1493.8834234708079,k4_estimate1:1474.177555852842,k4_estimate2:1143.5865622379004
## w:14,n:1000,m:14000,k4_true:1457.1963508671538,k4_estimate1:1440.8324628297094,k4_estimate2:1205.6516861700052
## w:15,n:1000,m:15000,k4_true:1424.7649580890852,k4_estimate1:1410.8661098095859,k4_estimate2:690.17906015712

**Tr(k4)_steps**

import numpy as np
import pandas as pd
n=1000
#m=2400
me1=[]
for w in range(15):
    m=n*(w+1)
    frq=[np.random.uniform(0.2,0.8) for i in range(m)]

    dat=np.zeros((n,m))
    for i in range(n):
        dat[i,:]=np.random.binomial(2, frq,m)
  
    x=pd.DataFrame(dat).apply(lambda x: (x-x.mean())/x.std())
    k=x.values @ x.values.T/x.shape[1]
    tril=[k[i,j] for i in range(n) for j in range(i)]
    me=1/np.var(tril)
    me1.append(me)

def k4_steps(n,me):
    k4_step11=1+6/me+3/me**2
    k4_step12=(1/(n-1)**2+1/me)*(1+1/me)
    k4_step13=1/(n-1)**4+6/(me*(n-1)**2)+3/me**2
    k4_step14=(1/(n-1)**2+1/me)**2
    k4_step21=(1+1/me)*(1/(n-1)**2+1/me)
    k4_step22=1/(n-1)**2+1/me
    k4_step23=(1/(n-1)**2+1/me)**2
    k4_step24=1/(n-1)**4
    k4_step25=-1/(n-1)**3
    df = pd.DataFrame([[k4_step11,k4_step12,k4_step13,k4_step14,k4_step21,k4_step22,k4_step23,k4_step24, k4_step25]],columns=['k4_step11','k4_step12','k4_step13','k4_step14','k4_step21','k4_step22','k4_step23','k4_step24', 'k4_step25'])
    return(df)

for me in me1:
   n=300
   df=k4_steps(n,me)
   print(df)

##    k4_step11  k4_step12  k4_step13  ...  k4_step23     k4_step24     k4_step25
## 0   1.005971   0.001007   0.000003  ...   0.000001  1.251167e-10 -3.740989e-08
## 
## [1 rows x 9 columns]
##    k4_step11  k4_step12     k4_step13  ...     k4_step23     k4_step24     k4_step25
## 0   1.002996   0.000511  7.811856e-07  ...  2.604786e-07  1.251167e-10 -3.740989e-08
## 
## [1 rows x 9 columns]
##    k4_step11  k4_step12     k4_step13  ...     k4_step23     k4_step24     k4_step25
## 0   1.001995   0.000344  3.541122e-07  ...  1.181208e-07  1.251167e-10 -3.740989e-08
## 
## [1 rows x 9 columns]
##    k4_step11  k4_step12     k4_step13  ...     k4_step23     k4_step24     k4_step25
## 0   1.001496   0.000261  2.032865e-07  ...  6.784558e-08  1.251167e-10 -3.740989e-08
## 
## [1 rows x 9 columns]
##    k4_step11  k4_step12     k4_step13  ...     k4_step23     k4_step24     k4_step25
## 0    1.00119    0.00021  1.314612e-07  ...  4.390380e-08  1.251167e-10 -3.740989e-08
## 
## [1 rows x 9 columns]
##    k4_step11  k4_step12     k4_step13  ...     k4_step23     k4_step24     k4_step25
## 0   1.000998   0.000178  9.435438e-08  ...  3.153487e-08  1.251167e-10 -3.740989e-08
## 
## [1 rows x 9 columns]
##    k4_step11  k4_step12     k4_step13  ...     k4_step23     k4_step24     k4_step25
## 0   1.000854   0.000154  7.043838e-08  ...  2.356287e-08  1.251167e-10 -3.740989e-08
## 
## [1 rows x 9 columns]
##    k4_step11  k4_step12     k4_step13  ...     k4_step23     k4_step24     k4_step25
## 0   1.000747   0.000136  5.502562e-08  ...  1.842528e-08  1.251167e-10 -3.740989e-08
## 
## [1 rows x 9 columns]
##    k4_step11  k4_step12     k4_step13  ...     k4_step23     k4_step24     k4_step25
## 0   1.000662   0.000122  4.407403e-08  ...  1.477475e-08  1.251167e-10 -3.740989e-08
## 
## [1 rows x 9 columns]
##    k4_step11  k4_step12     k4_step13  ...     k4_step23     k4_step24     k4_step25
## 0   1.000596    0.00011  3.635251e-08  ...  1.220092e-08  1.251167e-10 -3.740989e-08
## 
## [1 rows x 9 columns]
##    k4_step11  k4_step12     k4_step13  ...     k4_step23     k4_step24     k4_step25
## 0   1.000544   0.000102  3.084762e-08  ...  1.036595e-08  1.251167e-10 -3.740989e-08
## 
## [1 rows x 9 columns]
##    k4_step11  k4_step12     k4_step13  ...     k4_step23     k4_step24     k4_step25
## 0   1.000499   0.000094  2.642446e-08  ...  8.891565e-09  1.251167e-10 -3.740989e-08
## 
## [1 rows x 9 columns]
##    k4_step11  k4_step12     k4_step13  ...     k4_step23     k4_step24     k4_step25
## 0   1.000458   0.000088  2.272153e-08  ...  7.657255e-09  1.251167e-10 -3.740989e-08
## 
## [1 rows x 9 columns]
##    k4_step11  k4_step12     k4_step13  ...     k4_step23     k4_step24     k4_step25
## 0   1.000427   0.000082  2.012246e-08  ...  6.790897e-09  1.251167e-10 -3.740989e-08
## 
## [1 rows x 9 columns]
##    k4_step11  k4_step12     k4_step13  ...     k4_step23     k4_step24     k4_step25
## 0   1.000398   0.000078  1.779152e-08  ...  6.013917e-09  1.251167e-10 -3.740989e-08
## 
## [1 rows x 9 columns]

6 Estimate of the iteration of var(LB)

**var(LB)_estimate_var(LB)_times**

import numpy as np
import pandas as pd
def sum_list(ls):
    su=0
    for i in ls:
        su=su+i**2
    return su

def estimate_LB(x,B):
  n=x.shape[0]
  m=x.shape[1]
  LB=[]
  for i in range(B):
      zb=np.random.randn(n)
      L=x@(x.T@zb)
      LB.append(sum_list(L)/m**2)
  return(np.mean(LB),np.var(LB)/B)


n=1000
#m=10000
B=10
for w in range(20):
    m=n*(w+1)
    frq=[np.random.uniform(0.2,0.8) for i in range(m)]

    dat=np.zeros((n,m))
    for i in range(n):
        dat[i,:]=np.random.binomial(2, frq,m)
  
    x=pd.DataFrame(dat).apply(lambda x: (x-x.mean())/x.std())
    k=x.values @ x.values.T/x.shape[1]
    k4_trace=np.trace(k@k.T@k@k.T)
    lb=estimate_LB(x,B)
    #print('n:{},m:{},LB:{},var(LB):{},estimate_var(LB):{}'.format(n,m,lb[0],lb[1],(2*k4_trace)/B))
    if abs(lb[1]-(2*k4_trace)/B)/lb[0]<0.01:
        print('w:{},n:{},m:{},LB:{},var(LB):{},estimate_var(LB):{}'.format(w+1,n,m,lb[0],lb[1],(2*k4_trace)/B))
        break
    else:
        continue

## w:10,n:1000,m:10000,LB:1101.422945472793,var(LB):326.9165406421916,estimate_var(LB):331.60733856969966

7 Simulations of tr(K3) and tr(K4) by R

Simulation of tr(K3)

k3_solve=function(n,m){
  set.seed(1)
  frq=runif(m, 0.2, 0.8)
  dat=matrix(0, n, m)
  for(i in 1:n) {
    dat[i,]=rbinom(m, 2, frq)
  }
  x=apply(dat, 2, scale)
  G=x%*%t(x)/m
  
  F=G%*%t(G)
  K3=F%*%G
  target=sum(diag(K3))
  #estimate
  u11=1
  u12=-1/(n-1)
  
  #s11=var(G[col(G)==row(G)])
  #s12=var(G[col(G)<row(G)])
  s11=1/m
  s12=1/m
  #大样本时采用样本方差代替总体方差
  estimate1=u11^3+3*u11*s11
  estimate2=u11*(u12^2+s12)
  estimate3=u12^3
  estimate=n*estimate1+3*n*(n-1)*estimate2+n*(n-1)*(n-2)*estimate3
  target=round(target,4)
  estimate=round(estimate,4)
  error=round(abs(target-estimate),4)
  error_rate=round(error/target,4)
  return(c(target,estimate,error,error_rate))
}
n=c(1000,2000,3000)
m=c(10000,20000,30000)
#n=c(4000)
#m=c(60000)
True=c()
estimate=c()
error=c()
error_rate=c()
value=c()
for (i in seq(1,3,1)){
  re=k3_solve(n[i],m[i])
  True[i]=re[1] 
  estimate[i]=re[2]
  error[i]=re[3]
  error_rate[i]=re[4]
}
result=data.frame(n,m,True,estimate,error,error_rate)
result

##      n     m     True estimate   error error_rate
## 1 1000 10000 1308.070 1302.003  6.0665     0.0046
## 2 2000 20000 2617.973 2602.001 15.9716     0.0061
## 3 3000 30000 3927.763 3902.001 25.7618     0.0066

---

Simulation of tr(K4)

k4_solve=function(n,m,k){
  frq=runif(m, 0.2, 0.8)
  dat=matrix(0, n, m)
  for(i in 1:n) {
    dat[i,]=rbinom(m, 2, frq)
  }
  x=apply(dat, 2, scale)
  if (k==1)
  {
    G=x%*%t(x)/m
    F=G%*%t(G)
    K4=F%*%t(F)
    target=sum(diag(K4))
    return(target)
  }
  else
  {
    #estimate
    u11=1
    u12=-1/(n-1)
    #s11=var(G[col(G)==row(G)])
    #s12=var(G[col(G)<row(G)])
    s11=1/m
    s12=1/m
    diag1=u11^4+6*u11^2*s11+3*s11^2
    diag2=(u11^2+s11)*(u12^2+s12)
    diag3=u12^4+6*u12^2*s12+3*s12^2
    diag4=(u12^2+s12)^2
    estimate1=diag1+2*(n-1)*diag2+(n-1)*diag3+(n-1)*(n-2)*diag4
    not_diag1 = (u11^2+s11)*(u12^2+s12)
    not_diag2 = u12^2+s12
    not_diag3 = (u12^2+s12)^2
    not_diag4 = u12^4
    not_diag5 = u12^3
    estimate2=2*not_diag1+2*not_diag2+(n-2)*not_diag3+(n-2)*(n-3)*not_diag4+4*(n-2)*not_diag5
    estimate=n*estimate1+n*(n-1)*estimate2
    return(estimate)
  }
}
n=c(1000,2000,3000)
m=c(10000,20000,30000)
True=c()
estimate=c()
error=c()
error_rate=c()
for (i in seq(1,3,1)){
  True[i]=k4_solve(n[i],m[i],1) 
  estimate[i]=k4_solve(n[i],m[i],0)
  error[i]=abs(True[i]-estimate[i])
  error_rate[i]=error[i]/True[i]
}
result=data.frame(n,m,True,estimate,error,error_rate)
result

##      n     m     True estimate     error error_rate
## 1 1000 10000 1657.556 1623.417  34.13949 0.02059628
## 2 2000 20000 3319.216 3243.413  75.80306 0.02283764
## 3 3000 30000 4979.548 4863.412 116.13558 0.02332251

8 Simulation of heritablity

Simulate_h^2

solve_varg1=function(n,m,vg=0.5){
  frq=runif(m, 0.2, 0.8)
  dat=matrix(0, n, m)
  for(i in 1:n){
    dat[i,]=rbinom(m, 2, frq)
  }
  x=apply(dat, 2, scale)
  k=x%*%t(x)/m
  beta=rnorm(m,0,0.1)
  y1=x%*%beta
#  vg=0.5
  ve=var(y1)[1]*(1-vg)/vg
  e=rnorm(n,0,sqrt(ve))
  y=y1+e
  y=scale(y)
  estimate=(t(y)%*%k%*%y-t(y)%*%y)*m/n^2
  return(estimate)
}

#N=500
#M=2000
N=c(1000,2000,3000)
M=c(10000,20000,30000)
result=c()
for (k in seq(1,3,1)){
  for (i in seq(1,5,1)) {
    result[i]=solve_varg1(N[k],M[k])
  }
  print(c(N[k],M[k],mean(result)))
  result=c()
}

## [1] 1.000000e+03 1.000000e+04 4.810543e-01
## [1] 2.000000e+03 2.000000e+04 5.225054e-01
## [1] 3.000000e+03 3.000000e+04 5.054843e-01

9 Mailman case

import numpy as np
from numpy import random
#import itertools
from itertools import product


m=b=8
s=[0,1,2]
n=len(s)**m
A=random.randint(0,len(s),size=(m,n))
w=[x for x in product(s, repeat=b)]
u=np.array(w).T
#print(u.shape)
#print(A.shape)

p=np.zeros((n,n))

for i in range(A.shape[1]):
    for j in range(u.shape[1]):
        if (u[:,i]==A[:,j]).all():
            p[i][j]=1
#print(p)
x=np.random.randn(n,1)
z=p@x

un=np.hsplit(u,len(s))
zn=np.vsplit(z,len(s))

l=[]
q=[]
uns=[]
r=[]
for i in range(len(s)):
    l.append(un[i][0,:])
    uns.append(un[i][1:,:])
    q.append(l[i]@zn[i])
    r.append(uns[i]@zn[i])
wq=np.sum(q,axis=0)
w1=wq[:,np.newaxis]
w2=np.sum(r,axis=0)
w=np.concatenate((w1,w2))
print(w)

## [[ -93.5471264 ]
##  [ -29.52865911]
##  [  59.92585975]
##  [   5.93752644]
##  [ -55.37070067]
##  [-197.46114539]
##  [ -21.44760053]
##  [  19.93881951]]

def split_function(u,z):
    un=np.hsplit(u,len(s))
    zn=np.vsplit(z,len(s))
    zn_sum=np.sum(zn,axis=0)
    l=[]
    q=[]
    uns=[]
    #r=[]
    for i in range(len(s)):
        l.append(un[i][0,:])
        uns.append(un[i][1:,:])
        q.append(l[i]@zn[i])
       # r.append(uns[i]@zn[i])
    wq=np.sum(q,axis=0)
    w1=wq[:,np.newaxis]
 #   w2=uns[0]@zn_sum
 #   w=np.concatenate((w1,w2))
    return(w1,un[0][1:,:],zn_sum)
v=split_function(u,z)
print(v[0])

## [[-93.5471264]]

while v[1].shape[0]>1:
    b=split_function(v[1],v[2])
    print(b[0])
    v=b

## [[-29.52865911]]
## [[59.92585975]]
## [[5.93752644]]
## [[-55.37070067]]
## [[-197.46114539]]
## [[-21.44760053]]

final=v[1]@v[2]
print(final) 

## [[19.93881951]]

10 Conclusion

Reslut table

import sys
from prettytable import PrettyTable
import importlib
importlib.reload(sys)

## <module 'sys' (built-in)>

def solve_trace(n,me):
    k2=n+n**2/me
    k3=n+2+3*n**2/me
    k4=n+3+(6*n**2+4*n+4)/me+(2*n**3+2*n**2-n)/me**2
    return(k2,k3,k4)

def mse_B_solve(n,me,var_g,B,k):
    R=k[2]-4*k[1]+6*k[0]-3*n
    P=2*k[1]-6*k[0]+4*n
    Q=k[0]-n
    S=k[2]
    var=[]
    bias=[]
    mse=[]
    b=(2*R*var_g**2+2*P*var_g+2*Q)/Q**2
    for i in range(1,1+B):
        var.append(2*S*var_g**2/(i*Q**2)+(2*R*var_g**2+2*P*var_g+2*Q)/Q**2)
        bias.append(4*S**2*var_g**2/(i**2*Q**4))
        mse.append(2*S*var_g**2/(i*Q**2)+(2*R*var_g**2+2*P*var_g+2*Q)/Q**2+4*S**2*var_g**2/(i**2*Q**4))
        if (mse[i-1]-b<0.1*b):
            return(var[i-1],bias[i-1],mse[i-1],i,b)
            break
B=1000
##heritability
var_g=0.2

table = PrettyTable(['n','m','me','LB','tr(k3)','tr(k4)','var(h2)','bias','mse','B','b'])
h=0.2
n1=[10000,50000,100000,300000,500000]
m1=[3000000,5000000,7000000,10000000,15000000,20000000]
for j in range(len(n1)):    
    for i in range(len(m1)):
        n=n1[j]
        m=m1[i]
        me=0.1*m
        df=[n,m,me]       
        k=solve_trace(n, me)
        res=mse_B_solve(n,me,var_g,B,k)
        for p in range(len(k)):
             df.append(k[p])       
        for q in range(len(res)):
             df.append(res[q])
        table.add_row(df)   
print(table)

## +--------+----------+-----------+--------------------+--------------------+--------------------+------------------------+------------------------+------------------------+----+------------------------+
## |   n    |    m     |     me    |         LB         |       tr(k3)       |       tr(k4)       |        var(h2)         |          bias          |          mse           | B  |           b            |
## +--------+----------+-----------+--------------------+--------------------+--------------------+------------------------+------------------------+------------------------+----+------------------------+
## | 10000  | 3000000  |  300000.0 | 10333.333333333334 |      11002.0       |    12025.357791    |  0.006604114783487988  | 8.329491648090052e-06  |  0.006612444275136078  | 15 |  0.006026897609519991  |
## | 10000  | 5000000  |  500000.0 |      10200.0       |      10602.0       | 11211.080807960001 |  0.011021038207916528  | 2.375960924056883e-05  |  0.011044797817157097  | 23 |  0.010046161615920006  |
## | 10000  | 7000000  |  700000.0 | 10142.857142857143 | 10430.57142857143  | 10864.282046510203 |  0.015405900173017453  | 4.4280676742356486e-05 |  0.01545018084975981   | 32 |  0.014075025622319957  |
## | 10000  | 10000000 | 1000000.0 |      10100.0       |      10302.0       |   10605.04020399   |  0.022021662112629334  | 8.886271820503252e-05  |  0.022110524830834366  | 45 |    0.02013632163192    |
## | 10000  | 15000000 | 1500000.0 | 10066.666666666666 |      10202.0       | 10403.915647106667 |  0.03308156346356094   | 0.00019531205890312557 |  0.03327687552246407   | 67 |  0.03028648164792029   |
## | 10000  | 20000000 | 2000000.0 |      10050.0       |      10152.0       | 10303.520051997499 |  0.04420127808711001   | 0.0003431082757006707  |  0.04454438636281068   | 89 |   0.0404966416639199   |
## | 50000  | 3000000  |  300000.0 | 58333.333333333336 |      75002.0       | 102781.50001277777 | 0.00026689896961766394 | 1.4019575420272755e-08 | 0.0002669129891930842  | 5  | 0.00024321811201471994 |
## | 50000  | 5000000  |  500000.0 |      55000.0       |      65002.0       |   81003.4200078    | 0.0004402794788856685  | 3.4280772195056614e-08 | 0.00044031375965786353 | 7  | 0.00040324934402495994 |
## | 50000  | 7000000  |  700000.0 | 53571.42857142857  | 60716.28571428571  | 71942.07143418369  | 0.0006134315662613333  | 6.283953545428932e-08  | 0.0006134944057967876  | 9  |   0.0005632959360352   |
## | 50000  | 10000000 | 1000000.0 |      52500.0       |      57502.0       | 65253.205003949995 | 0.0008793256262369744  | 1.4413792732583228e-07 | 0.0008794697641643002  | 11 | 0.0008033946240505598  |
## | 50000  | 15000000 | 1500000.0 | 51666.666666666664 |      55002.0       | 60114.24666931111  | 0.0013190552576812394  |  3.33040679665349e-07  | 0.0013193882983609047  | 15 | 0.0012036359040761616  |
## | 50000  | 20000000 | 2000000.0 |      51250.0       |      53752.0       |  57565.6012519875  | 0.0017590973306334317  |  6.01587520929488e-07  | 0.0017596989181543612  | 19 | 0.0016039731841017602  |
## | 100000 | 3000000  |  300000.0 | 133333.33333333334 |      200002.0      | 322226.77778999996 | 6.740127400109998e-05  | 8.410237802318999e-10  | 6.740211502488021e-05  | 4  | 6.160119200087998e-05  |
## | 100000 | 5000000  |  500000.0 |      120000.0      |      160002.0      |   228003.8800076   | 0.00011072333120182401 | 2.0794307719408027e-09 | 0.00011072541063259595 | 5  | 0.00010160317600152001 |
## | 100000 | 7000000  |  700000.0 | 114285.71428571429 | 142859.14285714284 | 189799.5306177551  |  0.00015400635600252   | 3.8441463216156016e-09 | 0.0001540102001488416  | 6  |  0.00014160612000216   |
## | 100000 | 10000000 | 1000000.0 |      110000.0      |      130002.0      | 162003.42000389998 | 0.00022012701257499424 | 8.569831214027765e-09  | 0.00022013558240620827 | 7  | 0.00020161233600311996 |
## | 100000 | 15000000 | 1500000.0 | 106666.66666666667 |      120002.0      | 140892.16444706667 | 0.00032980592889413307 | 1.985060200257922e-08  | 0.00032982577949613567 | 9  | 0.0003016274960047198  |
## | 100000 | 20000000 | 2000000.0 |      105000.0      |      115002.0      |  130503.205001975  | 0.0004396132247341673  | 3.603271196728599e-08  | 0.0004396492574461346  | 11 |  0.00040164865600632   |
## | 300000 | 3000000  |  300000.0 |      600000.0      |     1200002.0      |   2700009.00001    | 7.800020666677779e-06  | 9.000060000166669e-12  | 7.800029666737779e-06  | 4  | 7.2000186666755565e-06 |
## | 300000 | 5000000  |  500000.0 |      480000.0      |      840002.0      |  1596006.1200068   | 1.2629678148169136e-05 | 2.4264932318807728e-11 | 1.2629702413101454e-05 | 4  | 1.1644489185201975e-05 |
## | 300000 | 7000000  |  700000.0 | 428571.4285714286  | 685716.2857142857  | 1181637.7346989796 |  1.75186073333642e-05  | 5.109646142799804e-11  | 1.751865842982563e-05  | 4  | 1.608897155558025e-05  |
## | 300000 | 10000000 | 1000000.0 |      390000.0      |      570002.0      |   894004.3800037   | 2.4963135555601234e-05 | 1.2181738019597626e-10 | 2.496325737298143e-05  | 4  | 2.2755717333369877e-05 |
## | 300000 | 15000000 | 1500000.0 |      360000.0      |      480002.0      | 684003.8800025333  | 3.6907036800067555e-05 | 2.310426211647013e-10  | 3.690726784268872e-05  | 5  | 3.386701955561185e-05  |
## | 300000 | 20000000 | 2000000.0 |      345000.0      |      435002.0      |  583503.645001925  | 4.882039516058256e-05  | 3.6902396719257034e-10 | 4.882076418454975e-05  | 6  | 4.4978395851927906e-05 |
## | 500000 | 3000000  |  300000.0 | 1333333.3333333335 |     3000002.0      | 8277793.0000077775 | 2.958403574401119e-06  | 1.420869225735474e-12  | 2.958404995270345e-06  | 4  | 2.720003136000895e-06  |
## | 500000 | 5000000  |  500000.0 |     1000000.0      |     2000002.0      |   4500009.000006   | 4.6800074400024005e-06 | 3.2400129600216002e-12 | 4.680010680015361e-06  | 4  |  4.32000672000192e-06  |
## | 500000 | 7000000  |  700000.0 | 857142.8571428572  | 1571430.5714285714 |  3153068.10204551  | 6.414412918403679e-06  | 6.110810658095267e-12  | 6.414419029214338e-06  | 4  | 5.920011840002943e-06  |
## | 500000 | 10000000 | 1000000.0 |      750000.0      |     1250002.0      |  2250005.5000035   | 9.040024160005602e-06  | 1.2960063360117762e-11 | 9.040037120068961e-06  | 4  |  8.32002240000448e-06  |
## | 500000 | 15000000 | 1500000.0 | 666666.6666666666  |     1000002.0      | 1611115.666669111  | 1.3480050960008802e-05 | 3.364019024037107e-11  | 1.3480084600199043e-05 | 4  | 1.2320047680007041e-05 |
## | 500000 | 20000000 | 2000000.0 |      625000.0      |      875002.0      | 1312504.125001875  | 1.7664086784011518e-05 | 4.515868385337509e-11  | 1.7664131942695372e-05 | 5  |  1.63200825600096e-05  |
## +--------+----------+-----------+--------------------+--------------------+--------------------+------------------------+------------------------+------------------------+----+------------------------+