R/gh.R

In FamEvent: Family Age-at-Onset Data Simulation and Penetrance Estimation

gh <- function(g,d,p){
xk<- c(-5.38748089001,-4.60368244955,-3.94476404012,-3.34785456738,-2.78880605843,-2.25497400209,-1.73853771212, -1.2340762154,-0.737473728545, -0.245340708301,0.245340708301, 0.737473728545, 1.2340762154,1.73853771212,2.25497400209, 2.78880605843,  3.34785456738, 3.94476404012,4.60368244955, 5.38748089001)
wk <- c(2.22939364553415129252E-13, 4.3993409922731805536E-10, 1.086069370769281694E-7, 7.80255647853206369415E-6, 2.28338636016353967257E-4, 0.00324377334223786183218, 0.0248105208874636108822, 0.1090172060200233200138, 0.2866755053628341297197, 0.46224366960061008965, 0.46224366960061008965, 0.28667550536283412972, 0.1090172060200233200138, 0.0248105208874636108822, 0.00324377334223786183218, 2.283386360163539672571E-4, 7.8025564785320636941E-6, 1.086069370769281694E-7, 4.39934099227318055363E-10, 2.22939364553415129252E-13)
yk <- exp(sqrt(2)*sqrt(p)*xk)
apply(cbind(d, g), 1, function(x) sum(wk*yk^x[1]*exp(-x[2]*yk))/sqrt(pi) )

#wt <- c(0.898591961453,0.704332961176, 0.62227869619, 0.575262442852,0.544851742366, 0.524080350949,0.509679027117,0.499920871336,0.493843385272,0.490921500667, 0.490921500667,0.493843385272,0.499920871336,0.509679027117,0.524080350949,0.544851742366, 0.575262442852, 0.62227869619, 0.704332961176,0.898591961453)
#apply(wt*exp(-t(as.matrix(g)%*%exp(as.matrix(t(sqrt(p)*xk)))) +d*sqrt(p)*xk - xk^2/2),2, sum)

}