This script requires the MixKin definitions. You have to evaluate MixKin prior to running these examples.

Examples of applications

Example : LRs for a random mixture, genotype prediction, LRs for relatives of the PoI, ROC curves, random man, familial search

Generate and inspect the mixture

DropoutProbs={0.2,0.3,0.7};
DropIn=0.05;

mix=RandomMixtureReplicates[
(* number of persons *)3,
(* replicates *)3,
(* dropout probs *)DropoutProbs,
(* dropins per locus expected *)DropIn];

(* display the mixture *)
a=Table[Alleles[[L,#]]&/@mix[[2,i,L]],{L,NGM},{i,Length[mix[[2]]]}];
Do[a[[i]]=PrependTo[a[[i]],Loci[[NGM[[i]]]]],{i,15}];
a//MatrixForm

(* display the donors *)
donors=Table[Alleles[[L,#]]&/@mix[[1,i,L]]//Sort,{L,NGM},{i,Length[mix[[1]]]}];
Do[donors[[i]]=PrependTo[donors[[i]],Loci[[NGM[[i]]]]],{i,15}];
donors//MatrixForm

Calculate standard LRs

(* calculate the LR for PoI=donor 1 *)
LR1=Monitor[LR[
(* the replicates *) mix[[2]],
(* Hp: donor 1 *){mix[[1,1]]},
(* Hd: no specified contributors *){},
(* number of donors *)3,
(* their respective dropout probabilities *)DropoutProbs,
(* drop in parameter *)DropIn,
(* kit *)NGM,
(* theta *)0.01,
(* print table after evaluation *)"Y"],res//MatrixForm]

(* Log 10 of the LR for donor 1 *)
Total[Log[10,LR1]]

(* equal dropout is faster but usually gives lower weight of evidence *)
(* LR for donor 1 with equal dropout for all donors, all the same as the real one for the first donor *)
LReq=Monitor[LR[mix[[2]],{mix[[1,1]]},{},3,{DropoutProbs[[1]],DropoutProbs[[1]],DropoutProbs[[1]]},DropIn,NGM,0.01,"Y"],res//MatrixForm];

(* LR for donor 1 with equal dropout for all donors, all the same as the real one for the first donor *)
Total[Log[10,LReq]]

(* H1: 3p mix with dropout (0.3, 0.3, 0.3) and c =0.5 *)
(* H2: 2p mix with dropout (0.1, 0.1) and c=0.5 *)
LRForNumberOfContributors[mix[[2]],{},{},3,{0.3,0.3,0.3},0.5,2,{0.1,0.1},0.5,NGM,"N"]

Total[Log[10,%]]

(* we calculate the LR as a function of the dropout probabilities, for fixed c=0.05 *)
AlgebraicLR=LR[mix[[2]],{mix[[1,1]]},{},3,{d1,d2,d3},0.05,NGM,0.01,"N"];//Timing

(* check: if an unknown donor has dropout one, then we get the LR for two individuals *)
AlgebraicLR/.{d1->0.1,d2->0.3,d3->1}
AlgebraicLR/.{d1->0.1,d2->1,d3->0.3}
LR[mix[[2]],{mix[[1,1]]},{},2,{0.1,0.3},0.05,NGM,0.01,"N"]

(* check: if everyone has dropout one, then all alleles must be contaminant and LR=1 on each locus *)
AlgebraicLR/.{d1->1,d2->1,d3->1}//Timing

(* sensitivity analysis: LR for donor 1 as a function of d1 if donors 2 and 3 have fixed dropout  *)
ListPlot[Table[{x,Total[Log[10,AlgebraicLR/.{d1->x,d2->DropoutProbs[[2]],d3->DropoutProbs[[3]]}]]},{x,Points}],PlotMarkers->Automatic,Joined->True,GridLines->Automatic,LabelStyle->Larger,AxesLabel->{"d1","Log(LR)"}]

(* set the dropout for third donor to 0.5 *)
AlgebraicLR2donors=AlgebraicLR/.{d3->0.5};

ListPlot[Table[{x,Total[Log[10,AlgebraicLR2donors/.{d1->x,d2->0.3}]]},{x,Points}],PlotMarkers->Automatic,Joined->True,GridLines->Automatic,LabelStyle->Larger,AxesLabel->{"d1","Log(LR)"}]//Timing

(* the same expression, but now simplified. This takes longer to compute, but... *)
AlgebraicLR2donors=AlgebraicLR/.{d3->DropoutProbs[[3]]}//Simplify;

(* the plots are now much easier to obtain *)
ListPlot[Table[{x,Total[Log[10,AlgebraicLR2donors/.{d1->x,d2->DropoutProbs[[2]]}]]},{x,Points}],PlotMarkers->Automatic,Joined->True,GridLines->Automatic,LabelStyle->Larger,AxesLabel->{"d1","Log(LR)"}]//Timing

AlgebraicLR=Simplify/@AlgebraicLR;

ListPlot[Table[{x,Total[Log[10,AlgebraicLR/.{d1->x,d2->0.3,d3->0.5}]]},{x,Points}],PlotMarkers->Automatic,Joined->True,GridLines->Automatic,LabelStyle->Larger,AxesLabel->{"d1","Log(LR)"}]//Timing

Integrating out the probability of dropout

DropoutProbs={0.2,0.3,0.7};
DropIn=0;

mix=RandomMixtureReplicates[
(* number of persons *)3,
(* replicates *)3,
(* dropout probs *)DropoutProbs,
(* dropins per locus expected *)DropIn];

donor 1

Liks=LRIncludingLikelihoods[
(* the replicates *)mix[[2]],
(* Hp *){mix[[1,1]]},
(* Hd *){},
(* number of donors *)3,
(* dropout probabilities *){d1,d2,d3},
(* Drop-in *)0,
(* kit *) NGM,
(* theta *)0,
(* do not show the expressions *)"N"];//Timing

(* integrate over the three dropout probabilities *)
(* may take a few up to 30 minutes *)
(* much slower if dropin >  0 since the algebraic LR is then a much larger expression *)
NIntegrate[Product[Liks[[1,L]],{L,15}],{d1,0,1},{d2,0,1},{d3,0,1}]/NIntegrate[Product[Liks[[2,L]],{L,15}],{d1,0,1},{d2,0,1},{d3,0,1}]//Quiet//Timing

(* for comparison: the LR obtained with the generating parameters *)
Monitor[LR[mix[[2]],{mix[[1,1]]},{},3,DropoutProbs,DropIn,NGM,0,"Y"],res//MatrixForm]

(* one free dropout parameter and the other two fixed makes the computation obviously much faster *)
Liks1U=Liks/.{d2->DropoutProbs[[2]],d3->DropoutProbs[[3]]}//Simplify;

NIntegrate[Product[Liks1U[[1,L]],{L,15}],{d1,0,1}]/NIntegrate[Product[Liks1U[[2,L]],{L,15}],{d1,0,1}]//Timing

(* two free parameters *)
Liks2U=Liks/.{d3->DropoutProbs[[3]]}//Simplify;

NIntegrate[Product[Liks2U[[1,L]],{L,15}],{d1,0,1},{d2,0,1}]/NIntegrate[Product[Liks2U[[2,L]],{L,15}],{d1,0,1},{d2,0,1}]//Timing

donor 2

Liks=LRIncludingLikelihoods[mix[[2]],{mix[[1,2]]},{},3,{d1,d2,d3},0,NGM,0,"N"]//Simplify;//Timing

(* one free dropout parameter. use the actual dropout rates for the two other donors *)
Liks1U=Liks/.{d1->DropoutProbs[[2]],d3->DropoutProbs[[3]]}//Simplify;

NIntegrate[Product[Liks1U[[1,L]],{L,15}],{d2,0,1}]/NIntegrate[Product[Liks1U[[2,L]],{L,15}],{d2,0,1}]

(* two free parameters *)
Liks2U=Liks/.{d3->DropoutProbs[[3]]}//Simplify;

NIntegrate[Product[Liks2U[[1,L]],{L,15}],{d1,0,1},{d2,0,1}]/NIntegrate[Product[Liks2U[[2,L]],{L,15}],{d1,0,1},{d2,0,1}]

(* three free parameters *)
NIntegrate[Product[Liks[[1,L]],{L,15}],{d1,0,1},{d2,0,1},{d3,0,1}]/NIntegrate[Product[Liks[[2,L]],{L,15}],{d1,0,1},{d2,0,1},{d3,0,1}]//Quiet

(* WoE obtained by integrating out all dropout parameters *)
Log[10,%]

(* WoE with the dropout parameters used to generate the mixture *)
Monitor[LR[mix[[2]],{mix[[1,2]]},{},3,{DropoutProbs[[2]],DropoutProbs[[1]],DropoutProbs[[3]]},DropIn,NGM,0,"Y"],res//MatrixForm]//Timing

donor 3

Liks=LRIncludingLikelihoods[mix[[2]],{mix[[1,3]]},{},3,{d1,d2,d3},0,NGM,0,"N"]//Simplify;//Timing

(* one free dropout parameter. use the actual dropout rates for the two other donors *)
Liks1U=Liks/.{d1->DropoutProbs[[3]],d2->DropoutProbs[[2]]}//Simplify;

NIntegrate[Product[Liks1U[[1,L]],{L,15}],{d3,0,1}]/NIntegrate[Product[Liks1U[[2,L]],{L,15}],{d3,0,1}]

(* two free parameters *)
Liks2U=Liks/.{d1->DropoutProbs[[3]]}//Simplify;

NIntegrate[Product[Liks2U[[1,L]],{L,15}],{d2,0,1},{d3,0,1}]/NIntegrate[Product[Liks2U[[2,L]],{L,15}],{d2,0,1},{d3,0,1}]

(* three free parameters *)
NIntegrate[Product[Liks[[1,L]],{L,15}],{d1,0,1},{d2,0,1},{d3,0,1}]/NIntegrate[Product[Liks[[2,L]],{L,15}],{d1,0,1},{d2,0,1},{d3,0,1}]//Quiet

(* WoE obtained by integrating out all dropout parameters *)
Log[10,%]

(* WoE with the dropout parameters used to generate the mixture *)
Monitor[LR[mix[[2]],{mix[[1,3]]},{},3,{DropoutProbs[[3]],DropoutProbs[[1]],DropoutProbs[[2]]},DropIn,NGM,0,"Y"],res//MatrixForm]//Timing

Calculate LRs opposing the PoI to unknown relatives instead of unrelated unknowns

(* this takes a little longer, up to 4 minutes for 3p *)
(* now in Hd we replace donor 1 by an unknown relative of his, of various types *)
(* so we get hypotheses Hp: donor 1 + 2 unknowns versus Hd: brother of donor 1 + 2 unknowns *)
LRrel1={{"sib","parent/child","halfsib","first cousin"},(Total[Log[10,LRSuspectOpposedToKin[mix[[2]],{mix[[1,1]]},3,DropoutProbs,DropIn,NGM,#]]])&/@{{0.25,0.5,0.25},{0,1,0},{0.5,0.5,0},{0.75,0.25,0}}};//Timing
LRrel1=Transpose[LRrel1];
AppendTo[LRrel1,{"unrelated",Total[Log[10,LR1]]}];Transpose[LRrel1]//MatrixForm

(* validation instrument *)
(* the sum of the sib-LR above and the Log[10,Mean[10^rm]] below should be equal to the LR for donor 1 *)
(* and indeed it almost exactly is *)
rm=LRSampleFromHdForRelatives[mix[[2]],{mix[[1,1]]},3,DropoutProbs,0.05,NGM,{0.25,0.5,0.25},1000000];
rm//Histogram
Log[10,Mean[10^rm]]+LRrel1[[1,2]]
Total[Log[10,LR[mix[[2]],{mix[[1,1]]},{},3,DropoutProbs,0.05,NGM,0,"N"]]]

(* sample 100,000 brothers of donor 1 and calculate the LR that they get as PoI with unequal parameters *)
(* sample 100,000 brothers of donor 1 and calculate the LR that they get as PoI with equal dropout 0.2*)
(* pair the histograms, Left: correct; Right: equal dropout *)
(* observe that brothers are excluded much better with the correct parameters if the real dropout values are quite different *)
s1=LRSampleFromHdForRelatives[mix[[2]],{mix[[1,1]]},3,DropoutProbs,0.05,NGM,{0.25,0.5,0.25},100000];
s2=LRSampleFromHdForRelatives[mix[[2]],{mix[[1,1]]},3,{DropoutProbs[[1]],DropoutProbs[[1]],DropoutProbs[[1]]},0.05,NGM,{0.25,0.5,0.25},100000]; 
PairedHistogram[s1,s2]

Length[Select[s1,#>0&]]
Length[Select[s2,#>0&]]

(* the mean LR should be one for a non - contributor *)
(* but very many simulations are needed to get close to that average *)
rm=LRSampleFromHdForRelatives[mix[[2]],{mix[[1,1]]},3,DropoutProbs,DropIn,NGM,{1,0,0},1000000];
rm//Histogram
Log[10,Mean[10^rm/.Indeterminate->0]]

LR distributions for all three donors

DropoutProbs={0.2,0.3,0.7};
DropIn=0.05;

mix=RandomMixtureReplicates[
(* number of persons *)3,
(* replicates *)3,
(* dropout probs *)DropoutProbs,
(* dropins per locus expected *)DropIn];

(* LRs to be expected when we look for donor 1 with the correct dropout probabilities *)
LRSampleFromHp[mix[[2]],{},3,DropoutProbs,0.05,NGM,100000]//Histogram

(* with wrong dropout probabilities *)
LRSampleFromHp[mix[[2]],{},3,{0.3,0.3,0.3},0.05,NGM,100000]//Histogram

(* now we suppose that we know donor 3 *)
LRSampleFromHp[mix[[2]],{mix[[1,3]]},3,{DropoutProbs[[3]],DropoutProbs[[1]],DropoutProbs[[2]]},0.05,NGM,100000]//Histogram

(* now we suppose that we know donors 2 and 3 *)
(* we now sometimes sample the same! Indeed we can better derive the donor's genotype now *)
LRSampleFromHp[mix[[2]],{mix[[1,2]],mix[[1,3]]},3,{DropoutProbs[[2]],DropoutProbs[[3]],DropoutProbs[[1]]},0.05,NGM,100000]//Histogram

(* LRs to be expected when we look for donor 2 *)
LRSampleFromHp[mix[[2]],{},3,{DropoutProbs[[2]],DropoutProbs[[1]],DropoutProbs[[2]]},0.05,NGM,100000]//Histogram

(* LRs to be expected when we look for donor 3 *)
LRSampleFromHp[mix[[2]],{},3,{DropoutProbs[[3]],DropoutProbs[[1]],DropoutProbs[[2]]},0.05,NGM,100000]//Histogram

Genotype probabilities

DropoutProbs={0.1,0.4,0.8};
DropIn=0;

mix=RandomMixtureReplicates[
(* number of persons *)3,
(* replicates *)3,
(* dropout probs *)DropoutProbs,
(* dropins per locus expected *)DropIn];

(* unconditional on any known donors *)
GenosDonor1=GenotypeDistribution[mix[[2]],{},3,DropoutProbs,DropIn,NGM];

GenosDonor1[[7]]//MatrixForm

(* genotypes of the second donor, conditional on the major *)
GenosDonor2=GenotypeDistribution[mix[[2]],{mix[[1,1]]},3,DropoutProbs,DropIn,NGM];

MatrixForm/@GenosDonor2

GenosMajorDonor=GenotypeDistribution[mix[[2]],{mix[[1,3]],mix[[1,2]]},3,{DropoutProbs[[3]],DropoutProbs[[2]],DropoutProbs[[1]]},DropIn,NGM];

MatrixForm/@GenosMajorDonor

donors=Table[Alleles[[L,#]]&/@mix[[1,i,L]]//Sort,{L,NGM},{i,Length[mix[[1]]]}];
Do[donors[[i]]=PrependTo[donors[[i]],Loci[[NGM[[i]]]]],{i,15}];
donors//MatrixForm

Using wrong models: random men and LR calculations

LR6=Monitor[LREqualDropout[mix[[2]],{mix[[1,1]]},{},2,0.2,0.05,NGM,"Y"],res//MatrixForm];
res//MatrixForm

Total[Log[10,LR6]]

Total[Log[10,LR1]]

rm=LRSampleFromHdForRelatives[mix[[2]],{mix[[1,1]]},2,{0.2,0.2},0.05,NGM,{1,0,0},100000];
rm//Histogram
Max[rm]
Quantile[rm,0.99]

rmhp=LRSampleFromHp[mix[[2]],{},2,{0.2,0.2},0.05,NGM,100000];rmhp//Histogram

Sum[If[a>Total[Log[10,LR6]],1/10^a,0],{a,rmhp}]/Length[rmhp]

LR7=Monitor[LREqualDropout[mix[[2]],{mix[[1,2]]},{},2,0.5,0.05,NGM,"Y"],res//MatrixForm];
res//MatrixForm

Sum[If[a>Total[Log[10,LR7]],1/10^a,0],{a,rmhp}]/Length[rmhp]

LR distributions for relatives of the donors

(* LRs for 10^5 brothers of donor 1 in favour of being the first donor *)
LRSampleFromHdForRelatives[mix[[2]],{mix[[1,1]]},3,DropoutProbs,DropIn,NGM,{0.25,0.5,0.25},100000]//Histogram

(* same but now evaluated with equal dropout: usually more false positives *)
AvDropout=Table[Mean[DropoutProbs],{k,Length[DropoutProbs]}];
LRSampleFromHdForRelatives[mix[[2]],{mix[[1,1]]},3,AvDropout,DropIn,NGM,{0.25,0.5,0.25},100000]//Histogram

(* cousins of donor 3, equal dropout *)
LRSampleFromHdForRelatives[mix[[2]],{mix[[1,3]]},3,AvDropout,DropIn,NGM,{0.75,0.25,0},100000]//Histogram

(* double first cousins of donor 2 *)
(* DFC: children when two brothers have children with two sisters *)
LRSampleFromHdForRelatives[mix[[2]],{mix[[1,2]]},3,AvDropout,DropIn,NGM,{9/16,6/16,1/16},100000]//Histogram

ROC curves for all three donors

(* ROC curve for donor 1 *)
MixtureROC[mix[[2]],{},3,DropoutProbs,DropIn,NGM,10000,Black]

(* ROC curve for donor 2 *)
MixtureROC[mix[[2]],{},3,{DropoutProbs[[2]],DropoutProbs[[1]],DropoutProbs[[3]]},DropIn,NGM,10000,Black]

(* ROC curve for donor 3 *)
MixtureROC[mix[[2]],{},3,{DropoutProbs[[3]],DropoutProbs[[1]],DropoutProbs[[2]]},DropIn,NGM,10000,Black]

LR comparison between correct and wrong dropout choices

(* LRs we get for the first donor, for two choices of evaluation: correct dropout and equal dropout *)
(* we take a sample from the Hp-distribution induced by the first specified model, i.e., the correct values for the probabilities of dropout *)
(* then we calculate the LRs that we would have obtained instead with the second set of parameters *)
s=LRSampleFromHpForDifferentParameters[
(* the replicates *) mix[[2]],
(* known contributors *) {},
(* first number of contributors *)3,
(* first set of dropout probabilities *)DropoutProbs,
(* first drop-in parameter *) DropIn,
(* second number of contributors *)3,
(* second set of dropout probabilities *)AvDropout,
(* second drop in parameter *)DropIn,
(* kit *)NGM,
(* sample size *)10000];

Show[ListPlot[s,GridLines->Automatic,AxesLabel->{"True LR","Equal Dropout LR"},LabelStyle->Larger],Plot[x,{x,0,Max[Transpose[s][[1]]]}]]

s=LRSampleFromHpForDifferentParameters[mix[[2]],{},3,DropoutProbs,0.05,3,{1,1,1},0.05,NGM,10000];

(* validation instrument: if the second set of parameters has dropout probabilities equal to one then all LR's are equal to one *)
Show[ListPlot[s,GridLines->Automatic,AxesLabel->{"True LR","Equal Dropout LR"},LabelStyle->Larger],Plot[x,{x,0,Max[Transpose[s][[1]]]}]]

s=LRSampleFromHpForDifferentParameters[mix[[2]],{},3,{1,1,1},0.05,3,DropoutProbs,0.05,NGM,10000];

(* validation instrument: if the first set of parameters has dropout probabilities equal to one then all LR's with the second set correspond to LRs from random persons *)
Show[ListPlot[s,GridLines->Automatic,AxesLabel->{"True LR","Equal Dropout LR"},LabelStyle->Larger],Plot[x,{x,0,Max[Transpose[s][[1]]]}]]

(* LRs we get for the second donor, for two choices of evaluation: correct dropout and equal dropout *)
s2=LRSampleFromHpForDifferentParameters[mix[[2]],{},3,{DropoutProbs[[2]],DropoutProbs[[1]],DropoutProbs[[3]]},DropIn,3,AvDropout,DropIn,NGM,10000];

Show[ListPlot[s2,GridLines->Automatic,AxesLabel->{"True LR","Equal Dropout LR"},LabelStyle->Larger],Plot[x,{x,0,Max[Transpose[s2][[1]]]}]]

(* LRs we get for the third donor, for two choices of evaluation: correct dropout and equal dropout *)
s3=LRSampleFromHpForDifferentParameters[mix[[2]],{},3,{DropoutProbs[[3]],DropoutProbs[[1]],DropoutProbs[[2]]},DropIn,3,AvDropout,DropIn,NGM,10000];

Show[ListPlot[s3,GridLines->Automatic,AxesLabel->{"True LR","Equal Dropout LR"},LabelStyle->Larger],Plot[x,{x,0,Max[Transpose[s3][[1]]]}]]

Familial Search

DropoutProbs={0.2,0.3,0.7};
DropIn=0.05;

mix=RandomMixtureReplicates[
(* number of persons *)3,
(* replicates *)3,
(* dropout probs *)DropoutProbs,
(* dropins per locus expected *)DropIn];

(* ROC curve for a familial search for parents/children of the first donor, without any donors assumed as known *)
FamilialSearchingROC[mix[[2]],{},3,DropoutProbs,DropIn,NGM,{0,1,0},10000,Green]

Hpsample=LRSampleFromFamilialHp[mix[[2]],{},3,DropoutProbs,DropIn,NGM,{0,1,0},10000];
Hpsample//Histogram

Length[Select[Hpsample,#>3&]]/Length[Hpsample]//N

HpsampleSibs=LRSampleFromFamilialHp[mix[[2]],{},3,DropoutProbs,DropIn,NGM,{0.25,0.5,0.25},10000];
HpsampleSibs//Histogram

Length[Select[HpsampleSibs,#>3&]]/Length[HpsampleSibs]//N

(* ROC curve for a familial search for parents/children of the second donor, assuming the first donor to be known *)
FamilialSearchingROC[mix[[2]],{mix[[1,1]]},3,DropoutProbs,DropIn,NGM,{0,1,0},10000,Green]

We assume here that the two minor donors are known and investigate the TPR and FPR for a familial search targeted at parents and children of donor 1 (major)

(* ROC curve for parents and children *)
FamilialSearchingROC[mix[[2]],{mix[[1,2]],mix[[1,3]]},3,{0.3,0.5,0.2},0.05,NGM,{0,1,0},10000,Black]

Example: KI - distribution with a given person

(* we regard a reference sample as a one-person mixture without dropin or dropout *)
mix=RandomMixtureReplicates[
(* number of persons *)1,
(* aantal replicas *)1,
(* dropoutkansen *){0},
(* dropin parameter *)0];

hp=LRSampleFromFamilialHp[mix[[2]],{},1,{0},0,NGM,{0.25,0.5,0.25},100000];//Timing

(* Histogram of SI with brothers, NGM kit, for the person generated in mix[[1,1]]*)
Histogram[hp]

hd=LRSampleFromFamilialHd[mix[[2]],{},1,{0},0,NGM,{0.25,0.5,0.25},100000];//Timing

(* Histogram of SI with urelated individuals, NGM kit, for the person generated in mix[[1,1]]*)
Histogram[hd]

(* ROC curve based on using Hp sample as importance sample for Hd *)
ListPlot[Table[{Log[10,Sum[If[hp[[i]]>x,10^-hp[[i]],0],{i,Length[hp]}]/Length[hp]],Length[Select[hp,#>x&]]/Length[hp]},{x,-5,10,0.2}],PlotRange->All,GridLines->Automatic,Joined->True,PlotMarkers->Automatic,AxesLabel->{"Log(10,FPR)","TPR"},LabelStyle->Larger]

(* ROC curves computed directly, using importance sampling with sample size 10^5, for sibs, parent/child and half-sibs *)
lp1=FamilialSearchingROC[mix[[2]],{},1,{0},0,NGM,{0.25,0.5,0.25},100000,Red]
lp2=FamilialSearchingROC[mix[[2]],{},1,{0},0,NGM,{0,1,0},100000,Green]
lp3=FamilialSearchingROC[mix[[2]],{},1,{0},0,NGM,{0.5,0.5,0},100000,Blue]
Show[lp1,lp2,lp3]

Distiguish half sibling from first cousins

(* create 100,000 first cousins and look at H1: first cousins; H2: half siblings, NGM *)
t1=Table[fc1=RandomProfile[32];fc2=CreateFirstCousin[fc1];Log[10,LRfamilial[{0.75,0.25,0},fc1,fc2,NGM]/LRfamilial[{0.5,0.5,0},fc1,fc2,NGM]],{i,100000}];
Mean[t1]
StandardDeviation[t1]
t1//Histogram

(* create 100,000 half-siblings and look at H1: half siblings; H2: first cousins, NGM *)
t2=Table[hs1=RandomProfile[32];hs2=CreateHalfSib[hs1];Log[10,LRfamilial[{0.5,0.5,0},hs1,hs2,NGM]/LRfamilial[{0.75,0.25,0},hs1,hs2,NGM]],{i,100000}];Mean[t2]
StandardDeviation[t2]
t2//Histogram

(* create 100,000 first cousins and look at H1: first cousins; H2: half siblings, 32 loci  *)
t3=Table[fc1=RandomProfile[32];fc2=CreateFirstCousin[fc1];Log[10,LRfamilial[{0.75,0.25,0},fc1,fc2,Range[32]]/LRfamilial[{0.5,0.5,0},fc1,fc2,Range[32]]],{i,100000}];Mean[t3]
StandardDeviation[t3]
t3//Histogram

(* create 100,000 half-siblings and look at H1: half siblings; H2: first cousins, 32 loci *)
t4=Table[hs1=RandomProfile[32];hs2=CreateHalfSib[hs1];Log[10,LRfamilial[{0.5,0.5,0},hs1,hs2,Range[32]]/LRfamilial[{0.75,0.25,0},hs1,hs2,Range[32]]],{i,100000}];Mean[t4]
StandardDeviation[t4]
t4//Histogram