This is MixKin 1.6.
No guarantees. No systematic support.
Written by Klaas Slooten (Netherlands Forensic Institute/ VU  University Amsterdam)
for information or questions.
A licence of Wolfram Mathematica is required to run this script.

The easiest is simply to select “Evaluation>Evaluate Notebook”. When prompted if you want to evaluate the Initialization Cells, best choose No. Otherwise some cells are evaluated twice, which is not problematic but unnecessary.
There are only definitions in this file, no examples. These are in separate files.

Please cite mixkin by citing the FSI:G article “Distinguishing between donors and their relatives in complex DNA mixtures with binary models”.

Change Log

Changes in 1.5.4. Added definition HypLikelihoodEqualDropoutTheta and HypLikelihoodTheta for a single likelihood

Changes in 1.5.3. Added definition LRinclLikelihoodsNoThetaFixedOrder

Changes in 1.5.2. w.r.t 1.5.1.
	- kit PPF6C added

Changes in 1.5.1. w.r.t 1.5
	- extra definitions for LR’s for familial searching

Changes in 1.5. w.r.t 1.4.3.1
- removed some duplicated definitions
- it was time for a round version number.

Changes in 1.4.3.1 w.r.t 1.4.3
 - for equal dropout LR's the resulting matrix-form LR table was output if the dropout probs were numerical, but without checking what the value of the monitor variable was. Fixed that.

Changes in 1.4.3. w.r.t 1.4.2
- LR definition now always with theta.  The computation automatically chooses the most efficient one, ie. when theta=0 the LR computation is done and when all dropout probabilities are the same, the LR for equal dropout is chosen.

Changes in 1.4.2 w.r.t 1.4.1
- theta added for standard LR calculations

Change in 1.4.1 w.r.t. 1.4.0.1
- user can change the allele frequency database, by choosing between the built in frequencies and supplying their own file in the same format
- kits NGM and SGMPlus are automatically defined from the data

Change in 1.4.0.1. wrt 1.4
- in definition of LRSameDonor we now only sum over genotypes with positive probability. Without this, there may be division by zero if there are alleles with frequency zero (which is not the case with the default frequencies).

Changes in 1.4. w.r.t. 1.3.2
- some tidying up removing unused cells
- slightly more comment
- probability of homozygous dropout can be chosen to be as described in "Estimating the probability of dropout of STR alleles in forensic genetics" (Tvedebrink et.al., FSI:G 2009) for LR computations (not for mixture simulations).

Changes in 1.3.2. w.r.t. 1.3.1.
- added definition LRSuspectOpposedToKinLikelihoodsSortedKit which computes the kinship-LR in the same way as for the standard LR, i.e., with the loci ordered according to increasing number of recorded alleles in the mixture. This makes comparisons of results between standard LR and kinship LR easier.

Changes in 1.3.1 w.r.t. 1.3:
- drop-in parameter c>1 allowed for Poisson model.

Changes in 1.3 w.r.t. 1.2.:
- Dropin model must be chosen: linear corresponds to (1-cp_a) and Poisson to Exp(-cp_a) as the term in the probability not to detect an allele. Dialog box opens.
- loci expanded with HdPlex, now 32 loci included. 
- mixture simulation function RandomMixtureReplicatesWithSomeDonorsFixed where some of the donors of the mixture can be fixed.
- function LRFirstUnknownThirdSib: computes LR for an unknown in a mixture to be a sibling of two persons who are known siblings.

1. Get allele frequencies, define kits

StartingTime=DateString[];

AlleleFreqChoice=ChoiceDialog["Do you want to use the built-in allele frequencies (NL-Caucasian 2085 men)?",{"Yes please"->"FLDO","No thanks, I import my own"->"Own"}]

(* Dutch allele frequencies based on a sample of 2085 individuals. The HdPlex loci are typed for only 335 individuals. *)
If[AlleleFreqChoice=="FLDO",RawDataAlleleFrequencies={{"Allele frequencies","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"D1S1656","","TPOX","","D2S441","","D2S1338","","D3S1358","","FGA","","D5S818","","CSF1PO","","SE33","","D7S820","","D8S1179","","D10S1248","","TH01","","VWA","","D12S391","","D13S317","","PENTA E","","D16S539","","D18S51","","D19S433","","PENTA D","","D21S11","","D22S1045","","D7S1517","","D3S1744","","D2S1360","","D6S474","","D4S2366","","D8S1132","","D5S2500","","D21S2055","","D10S2325",""},{10.`,0.0016786570743405277`,7.`,0.00047961630695443646`,8.`,0.00047961630695443646`,12.`,0.00023980815347721823`,10.`,0.00023980815347721823`,17.`,0.0011990407673860913`,7.`,0.00047961630695443646`,6.`,0.00023980815347721823`,10.2`,0.00023980815347721823`,7.`,0.01918465227817746`,8.`,0.01918465227817746`,7.`,0.00023980815347721823`,5.`,0.007194244604316546`,10.`,0.00023980815347721823`,14.`,0.00023980815347721823`,7.`,0.0014388489208633094`,5.`,0.07146282973621103`,7.`,0.00023980815347721823`,9.`,0.0009592326139088729`,10.`,0.00047961630695443646`,2.2`,0.00023980815347721823`,24.2`,0.00023980815347721823`,10.`,0.0014388489208633094`,16.`,0.0029850746268656717`,13.`,0.0029850746268656717`,18.`,0.0029850746268656717`,13.`,0.26865671641791045`,9.`,0.3492537313432836`,14.`,0.0014925373134328358`,9.`,0.0029850746268656717`,16.1`,0.056716417910447764`,7.`,0.14626865671641792`},{11.`,0.07482014388489208`,8.`,0.5414868105515588`,9.`,0.00023980815347721823`,14.`,0.00023980815347721823`,11.`,0.0019184652278177458`,18.`,0.017985611510791366`,8.`,0.003597122302158273`,7.`,0.0011990407673860913`,11.`,0.00047961630695443646`,8.`,0.16546762589928057`,9.`,0.016786570743405275`,10.`,0.0007194244604316547`,6.`,0.20983213429256595`,12.`,0.00047961630695443646`,15.`,0.041486810551558756`,8.`,0.10983213429256594`,6.`,0.0023980815347721825`,8.`,0.016786570743405275`,10.`,0.007913669064748202`,11.`,0.005035971223021582`,5.`,0.0007194244604316547`,25.2`,0.00023980815347721823`,11.`,0.1405275779376499`,17.`,0.004477611940298508`,14.`,0.10149253731343283`,19.`,0.013432835820895522`,14.`,0.19402985074626866`,10.`,0.18059701492537314`,16.`,0.011940298507462687`,10.`,0.09253731343283582`,17.1`,0.008955223880597015`,8.`,0.05373134328358209`},{12.`,0.12398081534772182`,9.`,0.09472422062350119`,10.`,0.18872901678657075`,15.`,0.00047961630695443646`,13.`,0.0059952038369304565`,19.`,0.060671462829736206`,9.`,0.030455635491606715`,8.`,0.0023980815347721825`,11.2`,0.00047961630695443646`,9.`,0.17985611510791366`,10.`,0.08896882494004796`,11.`,0.0026378896882494006`,7.`,0.18081534772182256`,13.`,0.0019184652278177458`,16.`,0.03741007194244604`,9.`,0.07410071942446043`,7.`,0.1707434052757794`,9.`,0.12949640287769784`,11.`,0.013908872901678658`,12.`,0.06690647482014389`,6.`,0.00047961630695443646`,26.`,0.0016786570743405277`,12.`,0.013429256594724222`,18.`,0.041791044776119404`,15.`,0.0626865671641791`,20.`,0.10746268656716418`,15.`,0.14626865671641792`,11.`,0.05970149253731343`,17.`,0.13134328358208955`,11.`,0.26119402985074625`,18.1`,0.010447761194029851`,9.`,0.0880597014925373`},{13.`,0.05731414868105515`,10.`,0.05731414868105515`,11.`,0.3431654676258993`,16.`,0.038129496402877695`,14.`,0.12278177458033573`,19.1`,0.00047961630695443646`,10.`,0.050599520383693045`,9.`,0.022541966426858514`,12.`,0.004556354916067146`,10.`,0.2592326139088729`,11.`,0.07745803357314149`,12.`,0.035491606714628296`,8.`,0.11103117505995204`,14.`,0.09880095923261392`,16.3`,0.00023980815347721823`,10.`,0.065947242206235`,8.`,0.010311750599520384`,10.`,0.06522781774580336`,12.`,0.15707434052757793`,12.1`,0.0009592326139088729`,7.`,0.006714628297362111`,27.`,0.042925659472422065`,13.`,0.0059952038369304565`,19.`,0.1253731343283582`,16.`,0.12835820895522387`,21.`,0.09253731343283582`,16.`,0.282089552238806`,12.`,0.15223880597014924`,18.`,0.2`,12.`,0.18358208955223881`,19.1`,0.2582089552238806`,10.`,0.10149253731343283`},{14.`,0.08752997601918466`,11.`,0.258273381294964`,11.3`,0.04244604316546763`,17.`,0.19400479616306954`,15.`,0.24652278177458034`,19.2`,0.0009592326139088729`,11.`,0.35203836930455634`,10.`,0.25491606714628295`,12.2`,0.00023980815347721823`,10.3`,0.00023980815347721823`,12.`,0.15947242206235013`,13.`,0.3160671462829736`,8.3`,0.00023980815347721823`,15.`,0.09760191846522781`,17.`,0.10047961630695443`,11.`,0.29640287769784174`,9.`,0.010551558752997603`,11.`,0.3215827338129496`,13.`,0.13788968824940048`,12.2`,0.00023980815347721823`,8.`,0.02038369304556355`,28.`,0.1693045563549161`,14.`,0.05083932853717027`,20.`,0.12686567164179105`,17.`,0.3194029850746269`,22.`,0.29850746268656714`,17.`,0.1044776119402985`,13.`,0.15671641791044777`,19.`,0.13134328358208955`,13.`,0.07014925373134329`,20.1`,0.05223880597014925`,11.`,0.16417910447761194`},{14.3`,0.0016786570743405277`,12.`,0.04628297362110311`,12.`,0.04748201438848921`,18.`,0.08369304556354916`,16.`,0.24172661870503595`,20.`,0.134052757793765`,12.`,0.37673860911270984`,10.3`,0.00023980815347721823`,13.`,0.007673860911270983`,11.`,0.19256594724220621`,12.1`,0.00023980815347721823`,14.`,0.30719424460431655`,9.`,0.13525179856115108`,16.`,0.20599520383693048`,17.3`,0.022781774580335732`,12.`,0.29376498800959233`,10.`,0.08800959232613909`,12.`,0.2642685851318945`,14.`,0.1647482014388489`,13.`,0.2134292565947242`,9.`,0.20863309352517986`,28.2`,0.00023980815347721823`,15.`,0.32517985611510786`,21.`,0.10298507462686567`,18.`,0.16119402985074627`,23.`,0.13582089552238805`,18.`,0.0029850746268656717`,14.`,0.09104477611940298`,20.`,0.13880597014925372`,14.`,0.08208955223880597`,21.`,0.0014925373134328358`,12.`,0.19552238805970149`},{15.`,0.13117505995203838`,13.`,0.0014388489208633094`,12.3`,0.0016786570743405277`,19.`,0.11654676258992806`,17.`,0.21630695443645084`,20.2`,0.00023980815347721823`,13.`,0.1748201438848921`,11.`,0.30695443645083936`,13.2`,0.0016786570743405277`,11.3`,0.00023980815347721823`,13.`,0.31966426858513186`,15.`,0.19640287769784173`,9.3`,0.34652278177458035`,17.`,0.2741007194244604`,18.`,0.17553956834532375`,13.`,0.11223021582733812`,11.`,0.09544364508393284`,13.`,0.1793764988009592`,15.`,0.14292565947242206`,13.2`,0.018225419664268584`,10.`,0.09904076738609113`,29.`,0.2011990407673861`,16.`,0.37074340527577937`,22.`,0.11641791044776119`,19.`,0.15223880597014924`,24.`,0.06119402985074627`,19.`,0.0014925373134328358`,15.`,0.010447761194029851`,21.`,0.12985074626865672`,15.`,0.21940298507462686`,21.1`,0.005970149253731343`,13.`,0.14626865671641792`},{15.3`,0.07482014388489208`,"","",13.`,0.02565947242206235`,20.`,0.1515587529976019`,18.`,0.15443645083932853`,21.`,0.16690647482014387`,14.`,0.00959232613908873`,12.`,0.3330935251798561`,13.3`,0.00047961630695443646`,12.`,0.1278177458033573`,14.`,0.19496402877697844`,16.`,0.11726618705035971`,10.`,0.009112709832134292`,18.`,0.19976019184652277`,18.3`,0.022062350119904078`,14.`,0.045563549160671464`,12.`,0.19784172661870503`,14.`,0.022541966426858514`,16.`,0.12326139088729017`,14.`,0.3661870503597122`,11.`,0.13621103117505995`,29.2`,0.0007194244604316547`,17.`,0.0841726618705036`,23.`,0.08059701492537313`,20.`,0.058208955223880594`,25.`,0.08507462686567165`,"","","","",22.`,0.11492537313432835`,16.`,0.06417910447761194`,22.`,0.0014925373134328358`,14.`,0.08059701492537313`},{16.`,0.11870503597122303`,"","",13.3`,0.00047961630695443646`,21.`,0.030455635491606715`,19.`,0.00935251798561151`,21.2`,0.003117505995203837`,15.`,0.0014388489208633094`,13.`,0.06402877697841727`,14.`,0.026858513189448444`,13.`,0.045563549160671464`,15.`,0.09184652278177458`,17.`,0.02038369304556355`,"","",19.`,0.10599520383693045`,19.`,0.10311750599520383`,15.`,0.0007194244604316547`,13.`,0.09376498800959233`,15.`,0.00023980815347721823`,17.`,0.10743405275779377`,14.2`,0.02326139088729017`,12.`,0.23213429256594725`,29.3`,0.00047961630695443646`,18.`,0.006714628297362111`,24.`,0.14477611940298507`,21.`,0.013432835820895522`,26.`,0.09850746268656717`,"","","","",23.`,0.1044776119402985`,17.`,0.019402985074626865`,22.1`,0.0014925373134328358`,15.`,0.022388059701492536`},{16.1`,0.00023980815347721823`,"","",14.`,0.29376498800959233`,22.`,0.03141486810551559`,20.`,0.0007194244604316547`,22.`,0.17218225419664268`,16.`,0.00023980815347721823`,14.`,0.012470023980815348`,14.2`,0.0014388489208633094`,14.`,0.009112709832134292`,16.`,0.026378896882494004`,18.`,0.0033573141486810554`,"","",20.`,0.014388489208633093`,19.3`,0.009832134292565946`,"","",14.`,0.06139088729016787`,16.`,0.00023980815347721823`,18.`,0.07170263788968825`,15.`,0.18297362110311752`,13.`,0.21223021582733814`,30.`,0.2587529976019185`,19.`,0.0009592326139088729`,25.`,0.21044776119402986`,"","",27.`,0.04776119402985075`,"","","","",24.`,0.03283582089552239`,18.`,0.004477611940298508`,23.`,0.011940298507462687`,17.`,0.0014925373134328358`},{16.3`,0.05491606714628297`,"","",14.1`,0.00023980815347721823`,23.`,0.09376498800959233`,"","",22.1`,0.00023980815347721823`,"","",15.`,0.0016786570743405277`,14.3`,0.00023980815347721823`,15.`,0.0007194244604316547`,17.`,0.005035971223021582`,19.`,0.00023980815347721823`,"","",21.`,0.0007194244604316547`,20.`,0.12062350119904076`,"","",15.`,0.050359712230215826`,"","",19.`,0.03669064748201439`,15.2`,0.04220623501199041`,13.4`,0.0011990407673860913`,30.2`,0.0328537170263789`,"","",26.`,0.03283582089552239`,"","",28.`,0.028358208955223882`,"","","","",25.`,0.0029850746268656717`,"","",24.`,0.008955223880597015`,"",""},{17.`,0.05779376498800959`,"","",15.`,0.04964028776978417`,24.`,0.12278177458033573`,"","",22.2`,0.007434052757793765`,"","",16.`,0.00023980815347721823`,15.`,0.03525179856115108`,"","","","","","","","","","",20.2`,0.00023980815347721823`,"","",16.`,0.05227817745803358`,"","",20.`,0.019664268585131893`,16.`,0.053237410071942444`,14.`,0.06786570743405276`,31.`,0.09160671462829736`,"","",27.`,0.008955223880597015`,"","",29.`,0.016417910447761194`,"","","","","","","","",25.`,0.13134328358208955`,"",""},{17.1`,0.0007194244604316547`,"","",16.`,0.0059952038369304565`,25.`,0.11247002398081535`,"","",23.`,0.147242206235012`,"","","","",15.2`,0.00047961630695443646`,"","","","","","","","","","",20.3`,0.00047961630695443646`,"","",17.`,0.049160671462829736`,"","",21.`,0.009112709832134292`,16.2`,0.016786570743405275`,14.1`,0.0007194244604316547`,31.2`,0.07793764988009592`,"","",28.`,0.0014925373134328358`,"","",30.`,0.008955223880597015`,"","","","","","","","",26.`,0.14776119402985075`,"",""},{17.3`,0.14196642685851318`,"","","","",26.`,0.021582733812949638`,"","",23.2`,0.0033573141486810554`,"","","","",15.3`,0.0011990407673860913`,"","","","","","","","","","",21.`,0.13093525179856116`,"","",18.`,0.02422062350119904`,"","",22.`,0.004316546762589928`,17.`,0.005275779376498801`,15.`,0.009832134292565946`,32.`,0.01606714628297362`,"","","","","","",31.`,0.0014925373134328358`,"","","","","","","","",27.`,0.019402985074626865`,"",""},{18.`,0.0038369304556354917`,"","","","",27.`,0.0016786570743405277`,"","",24.`,0.13764988009592327`,"","","","",16.`,0.04364508393285372`,"","","","","","","","","","",22.`,0.10143884892086331`,"","",19.`,0.010551558752997603`,"","",23.`,0.0009592326139088729`,17.2`,0.0033573141486810554`,16.`,0.003117505995203837`,32.2`,0.07338129496402877`,"","","","","","",32.`,0.0014925373134328358`,"","","","","","","","",28.`,0.013432835820895522`,"",""},{18.3`,0.05371702637889689`,"","","","",28.`,0.0009592326139088729`,"","",24.2`,0.0011990407673860913`,"","","","",16.2`,0.00023980815347721823`,"","","","","","","","","","",23.`,0.07505995203836931`,"","",20.`,0.007434052757793765`,"","",24.`,0.0014388489208633094`,18.`,0.00023980815347721823`,17.`,0.00023980815347721823`,33.`,0.0028776978417266188`,"","","","","","","","","","","","","","","","",29.`,0.016417910447761194`,"",""},{19.`,0.00047961630695443646`,"","","","","","","","",25.`,0.1026378896882494`,"","","","",16.3`,0.0009592326139088729`,"","","","","","","","","","",24.`,0.036930455635491605`,"","",21.`,0.0028776978417266188`,"","","","",18.2`,0.0009592326139088729`,18.`,0.00023980815347721823`,33.2`,0.02613908872901679`,"","","","","","","","","","","","","","","","",30.`,0.020895522388059702`,"",""},{19.3`,0.013669064748201438`,"","","","","","","","",26.`,0.03741007194244604`,"","","","",17.`,0.06354916067146282`,"","","","","","","","","","",25.`,0.015587529976019185`,"","",22.`,0.0007194244604316547`,"","","","",19.2`,0.00023980815347721823`,"","",34.`,0.00023980815347721823`,"","","","","","","","","","","","","","","","",31.`,0.026865671641791045`,"",""},{20.3`,0.0009592326139088729`,"","","","","","","","",27.`,0.004796163069544365`,"","","","",17.2`,0.00023980815347721823`,"","","","","","","","","","",26.`,0.003597122302158273`,"","",23.`,0.00047961630695443646`,"","","","","","","","",34.1`,0.00047961630695443646`,"","","","","","","","","","","","","","","","",32.`,0.029850746268656716`,"",""},{"","","","","","","","","","",28.`,0.00023980815347721823`,"","","","",17.3`,0.0009592326139088729`,"","","","","","","","","","",27.`,0.0014388489208633094`,"","","","","","","","","","","","",34.2`,0.0019184652278177458`,"","","","","","","","","","","","","","","","",33.`,0.058208955223880594`,"",""},{"","","","","","","","","","","","","","","","",18.`,0.0671462829736211`,"","","","","","","","","","",28.`,0.00023980815347721823`,"","","","","","","","","","","","",35.`,0.00047961630695443646`,"","","","","","","","","","","","","","","","",34.`,0.06567164179104477`,"",""},{"","","","","","","","","","","","","","","","",18.2`,0.0007194244604316547`,"","","","","","","","","","",29.`,0.00023980815347721823`,"","","","","","","","","","","","",35.2`,0.00023980815347721823`,"","","","","","","","","","","","","","","","",35.`,0.04029850746268657`,"",""},{"","","","","","","","","","","","","","","","",18.3`,0.0009592326139088729`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",36.`,0.007462686567164179`,"",""},{"","","","","","","","","","","","","","","","",19.`,0.06258992805755396`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",37.`,0.004477611940298508`,"",""},{"","","","","","","","","","","","","","","","",19.1`,0.00047961630695443646`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",19.2`,0.005275779376498801`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",20.`,0.04844124700239808`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",20.1`,0.00047961630695443646`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",20.2`,0.015587529976019185`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",20.3`,0.00023980815347721823`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",21.`,0.03213429256594724`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",21.2`,0.017745803357314148`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",22.`,0.008633093525179856`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",22.2`,0.02853717026378897`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",23.`,0.0026378896882494006`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",23.2`,0.027098321342925658`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",24.`,0.0007194244604316547`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",24.2`,0.03213429256594724`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",25.2`,0.035731414868105514`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",25.3`,0.0009592326139088729`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",26.2`,0.06354916067146282`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",27.`,0.00047961630695443646`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",27.2`,0.08513189448441247`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",27.3`,0.00023980815347721823`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",28.2`,0.08321342925659472`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",29.2`,0.07218225419664269`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",30.2`,0.056594724220623505`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",31.`,0.00023980815347721823`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",31.2`,0.026618705035971222`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",32.`,0.0019184652278177458`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",32.2`,0.011031175059952039`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",33.`,0.00959232613908873`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",33.2`,0.004556354916067146`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",34.`,0.0023980815347721825`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",34.2`,0.0007194244604316547`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",35.`,0.0011990407673860913`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",35.2`,0.00047961630695443646`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",36.`,0.00047961630695443646`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""},{"","","","","","","","","","","","","","","","",37.`,0.00023980815347721823`,"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",""}}];

If[AlleleFreqChoice=="Own",RawDataAlleleFrequenciesFile=SystemDialogInput["FileOpen"];RawDataAlleleFrequencies=Import[RawDataAlleleFrequenciesFile][[1]];]

PopulationStatistics=If[AlleleFreqChoice=="Own",RawDataAlleleFrequenciesFile,"FLDO 2085"]

NumberOfLoci=Length[RawDataAlleleFrequencies[[2]]]/2

Loci=Table[RawDataAlleleFrequencies[[2,2*i-1]],{i,NumberOfLoci}]

(* change spelling of vWA *)
Loci=Loci/.{"VWA"->"vWA"}

(* define kits *)
NGM=Position[Loci,#]&/@{"D10S1248","vWA","D16S539","D2S1338","D8S1179","D21S11","D18S51","D22S1045","D19S433","TH01","FGA","D2S441","D3S1358","D1S1656","D12S391"}//Flatten
SGMPlus=Position[Loci,#]&/@{"D3S1358","vWA","D16S539","D2S1338","D8S1179","D21S11","D18S51","D19S433","TH01","FGA"}//Flatten
Identifiler=Position[Loci,#]&/@{"D8S1179","D21S11","D7S820","CSF1PO","D3S1358","TH01","D13S317","D16S539","D2S1338","D19S433","vWA","TPOX","D18S51","D5S818","FGA"}//Flatten
MiniFiler=Position[Loci,#]&/@{"D13S317","D7S820","D2S1338","D21S11","D16S539","D18S51","CSF1PO","FGA"}//Flatten
PPF6C=Position[Loci,#]&/@{"D3S1358","D1S1656","D2S441","D10S1248","D13S317","PENTA E","D16S539","D18S51","D2S1338","CSF1PO","PENTA D","TH01","vWA","D21S11","D7S820","D5S818","TPOX","D8S1179","D12S391","D19S433","SE33","D22S1045","FGA"}//Flatten

TypedKits={};FailedKits={};
If[Length[NGM]==15,AppendTo[TypedKits,"NGM"],AppendTo[FailedKits,"NGM"]];
If[Length[SGMPlus]==10,AppendTo[TypedKits,"SGMPlus"],AppendTo[FailedKits,"SGMPlus"]];
If[Length[Identifiler]==15,AppendTo[TypedKits,"Identifiler"],AppendTo[FailedKits,"Identifiler"]];
If[Length[MiniFiler]==8,AppendTo[TypedKits,"MiniFiler"],AppendTo[FailedKits,"MiniFiler"]];
If[Length[PPF6C]==23,AppendTo[TypedKits,"PPF6C"],AppendTo[FailedKits,"PPF6C"]];

TypedKits

If[Length[FailedKits]>0,CreateDialog[Column[{"The kits ", FailedKits," have not (completely) been found.",DefaultButton[]}]]];

(* Ladder lengths *)
Quiet[LL=Table[0,{i,1,NumberOfLoci}];For[i=1,i≤NumberOfLoci, i++,L=2;While[RawDataAlleleFrequencies[[L+2,2*i-1]]≠"" ,L++];LL[[i]]=L-1]];LL

(* tables with alleles, denoted by repeat number *)
Alleles=Table[0,{i,NumberOfLoci},{j,1,1}];For[i=1,i≤NumberOfLoci,i++,Alleles[[i]]=Table[0,{j,LL[[i]]}];For[k=1,k≤LL[[i]],k++,Alleles[[i,k]]=RawDataAlleleFrequencies[[k+2,2*i-1]]]]

(* the numbers 1,..,k for a ladder with k alleles, on all loci *)
AlleleNumbers=Table[Range[LL[[i]]],{i,1,NumberOfLoci}];

Fr=Table[RawDataAlleleFrequencies[[k+2,2*i]],{i,NumberOfLoci},{k,LL[[i]]}];

GiveAlleleNumber[L_,a_]:=Sum[KroneckerDelta[Alleles[[L,i]],a]*i,{i,LL[[L]]}]

(* some points between zero and one for which we will calculate LRs when they represent dropout probabilities *)
Points=Join[Range[9]/100,Range[9]/10,Range[9]/100+0.9]//N;

(* same, but fewer points *)
FewerPoints=Join[{0.01},Range[9]/10,{0.99}]//N;

Kinship

Sibling index

(* SI on one locus, no mutation *)OneLocusSI[{a1_,a2_,b1_,b2_},i_]:=1/4+(KroneckerDelta[a1,b1]*(KroneckerDelta[a2,b2]+Fr[[i,b2]]) +KroneckerDelta[a1,b2]*(KroneckerDelta[a2,b1]+Fr[[i,b1]])+KroneckerDelta[a2,b1]*Fr[[i,b2]]+KroneckerDelta[a2,b2]*Fr[[i,b1]])/(8Fr[[i,b1]]*Fr[[i,b2]])

SI[list1_,list2_,kit_]:=Product[OneLocusSI[{list1[[2*i-1]],list1[[2*i]],list2[[2*i-1]],list2[[2*i]]},i],{i,kit}]

(* create a sibling's genotype on locus i for a person with genotype (a,b) *)
CreateSibOneLocus[{a_,b_},i_]:=Module[{ibd=RandomChoice[{none,left,right,both}], result}, If[a==0,result={0,0},If[ibd==none,result={RandomChoice[Fr[[i]]-> AlleleNumbers[[i]]],RandomChoice[Fr[[i]]-> AlleleNumbers[[i]]]}]; If[ibd==left,result={a,RandomChoice[Fr[[i]]-> AlleleNumbers[[i]]]}]; If[ibd==right,result={b,RandomChoice[Fr[[i]]-> AlleleNumbers[[i]]]}]; If[ibd==both,result={a,b}]];result];

(* create a sib on all NumberOfLoci loci *)CreateSib[list_]:=Module[{nieuwesib={}},For[i=1,i≤NumberOfLoci,i++,nieuwesib=Join[nieuwesib,CreateSibOneLocus[{list[[2*i-1]],list[[2*i]]},i]]];nieuwesib];

Parent - Child Index

(* PI on locus i between persons with genotypes (a1,a2) and (b1,b2), no mutation *)OneLocusPI[{a1_,a2_,b1_,b2_},i_]:=(KroneckerDelta[a1,b1]+KroneckerDelta[a2,b1])/(4Fr[[i,b1]])+(KroneckerDelta[a1,b2]+KroneckerDelta[a2,b2])/(4Fr[[i,b2]])

(* PI on kit, a set of markers *)
PI[list1_,list2_,kit_]:=Product[OneLocusPI[{list1[[2*i-1]],list1[[2*i]],list2[[2*i-1]],list2[[2*i]]},i],{i,kit}]

(* create genotype of a parent or child on locus i for a person with genotype (a,b) heeft, disregarding mutation *)
CreateChildOneLocus[{a_,b_},i_]:=Module[{ibd=RandomChoice[{left,right}], result}, If[a==0,result={0,0}, If[ibd==left,result={a,RandomChoice[Fr[[i]]-> AlleleNumbers[[i]]]}]; If[ibd==right,result={b,RandomChoice[Fr[[i]]-> AlleleNumbers[[i]]]}]];result];

(* create a child on all NumberOfLoci loci *)CreateChild[list_]:=Module[{nieuwkind={}},For[i=1,i≤NumberOfLoci,i++,nieuwkind=Join[nieuwkind,CreateChildOneLocus[{list[[2*i-1]],list[[2*i]]},i]]];nieuwkind];

Half - sibling index

OneLocusHSI[{a1_,a2_,b1_,b2_},i_]:=1/2+(KroneckerDelta[a1,b1]+KroneckerDelta[a2,b1])/(8Fr[[i,b1]])+(KroneckerDelta[a1,b2]+KroneckerDelta[a2,b2])/(8Fr[[i,b2]])

HSI[list1_,list2_,kit_]:=Product[OneLocusHSI[{list1[[2*i-1]],list1[[2*i]],list2[[2*i-1]],list2[[2*i]]},i],{i,kit}]

(* create a halfsibling's genotype on locus i for a person with genotype (a,b) *)
CreateHalfSibOneLocus[{a_,b_},i_]:=Module[{ibd=RandomChoice[{0.5,0.25,0.25}->{none,left,right}], result}, If[a==0,result={0,0},If[ibd==none,result={RandomChoice[Fr[[i]]-> AlleleNumbers[[i]]],RandomChoice[Fr[[i]]-> AlleleNumbers[[i]]]}]; If[ibd==left,result={a,RandomChoice[Fr[[i]]-> AlleleNumbers[[i]]]}]; If[ibd==right,result={b,RandomChoice[Fr[[i]]-> AlleleNumbers[[i]]]}]];result];

(* create a halfsib on all NumberOfLoci loci *)CreateHalfSib[list_]:=Module[{nieuwehalfsib={}},For[i=1,i≤NumberOfLoci,i++,nieuwehalfsib=Join[nieuwehalfsib,CreateHalfSibOneLocus[{list[[2*i-1]],list[[2*i]]},i]]];nieuwehalfsib];

First cousins

(* create a halfsibling's genotype on locus i for a person with genotype (a,b) *)
CreateFirstCousinOneLocus[{a_,b_},i_]:=Module[{ibd=RandomChoice[{0.75,0.125,0.125}->{none,left,right}], result}, If[a==0,result={0,0},If[ibd==none,result={RandomChoice[Fr[[i]]-> AlleleNumbers[[i]]],RandomChoice[Fr[[i]]-> AlleleNumbers[[i]]]}]; If[ibd==left,result={a,RandomChoice[Fr[[i]]-> AlleleNumbers[[i]]]}]; If[ibd==right,result={b,RandomChoice[Fr[[i]]-> AlleleNumbers[[i]]]}]];result];

(* create a halfsib on all NumberOfLoci loci *)CreateFirstCousin[list_]:=Module[{NewCousin={}},For[i=1,i≤NumberOfLoci,i++,NewCousin=Join[NewCousin,CreateFirstCousinOneLocus[{list[[2*i-1]],list[[2*i]]},i]]];NewCousin];

LR for kinship using IBD coefficients

(* given two persons with genotypes (a1, a2) and (b1, b2) on locus i;LR for H1: they are related according to IBD coefficents (k0,k1,k2) versus H2: unrelated *) LRfamilialOneLocus[{k0_,k1_,k2_},{a1_,a2_,b1_,b2_},i_]:=k0+k1*((KroneckerDelta[a1,b1]+KroneckerDelta[a2,b1])/(4Fr[[i,b1]])+(KroneckerDelta[a1,b2]+KroneckerDelta[a2,b2])/(4Fr[[i,b2]]))+k2*((KroneckerDelta[a1,b1]KroneckerDelta[a2,b2]+KroneckerDelta[a1,b2]KroneckerDelta[a2,b1])/(2(Fr[[i,b1]]Fr[[i,b2]])))
(* same on a kit *)
LRfamilial[vector_,list1_,list2_,kit_]:=Product[LRfamilialOneLocus[vector,{list1[[2*i-1]],list1[[2*i]],list2[[2*i-1]],list2[[2*i]]},i],{i,kit}]

Child of two parents

Child[p1_,p2_]:=Module[{res={}},Do[AppendTo[res,{RandomChoice[{p1[[2L-1]],p1[[2L]]}],RandomChoice[{p2[[2L-1]],p2[[2L]]}]}],{L,NumberOfLoci}];Flatten[res]]

Mixture simulations and LR definitions

Mixture simulation

(* dialog to choose the homozygous dropout probability *)
HomDropout=ChoiceDialog["Choose relation between homozygous and heterozygous dropout probability",{"Square"->"Sq","Logistic regression (Tvedebrink 2009 FSI:G)"->"Tved"}]

(* choose the parameter for tvedebrink model if it has been selected *)
If[HomDropout=="Tved",beta=Input["Select the value of beta to use in the formula D=2^beta*d/(1+(2^beta-1)d)"]]

(* check value: it must be a negative number *)
If[HomDropout=="Tved"&&(Not[NumberQ[beta]]||beta≥0),MessageDialog["beta has to be a negative number. Script aborts "];Abort[]];

(* definition of homozygous dropout probability in terms of heterozygous dropout probability *)
(* from Cowell et. al. Analysis of DNA mixtures with artefacts *)
Tved[d_,beta_]:=2^beta*d/(1+(2^beta-1)d)

(* dialog to choose the drop-in model in the LR evaluation. Linear is the first order approximation of Poisson *)
(* Poisson will enable all c≥0 and linear only small c *)
DropinModel=ChoiceDialog["Pick drop-in model",{Linear->"Linear",Poisson->"Poisson"}]

KD[{x_,y_}]:=KroneckerDelta[x,y]

(* random match probability *)
RMP[profile_,kit_]:=Product[(2-KD[{profile[[L,1]],profile[[L,2]]}])Fr[[L,profile[[L,1]]]]Fr[[L,profile[[L,2]]]],{L,kit}]

(* generate a random profile according to the allele frequencies *)RandomProfile[FirstSoManyLoci_]:=RandomChoice[Fr[[#]]->Range[LL[[#]]]]&/@Table[Ceiling[i/2],{i,2*FirstSoManyLoci}]

(* create a mixture with n donors, in rep replicates given that their heterozygous dropout probabilities are in the vector het *)
RandomMixtureNoContamination[n_,rep_,het_]:=Module[{MixtureDonors=Table[{},{t,n}],profiles=Table[Partition[RandomProfile[NumberOfLoci],2],{t,n}],mix=Table[{},{L,NumberOfLoci},{r,rep}]},MixtureDonors=profiles;
Do[Do[Do[If[profiles[[k,L,1]]==profiles[[k,L,2]],profiles[[k,L]]=RandomChoice[{het[[k]]^2,1-het[[k]]^2}->{{},profiles[[k,L]]}],profiles[[k,L]]=RandomChoice[{het[[k]]^2,het[[k]](1-het[[k]]),het[[k]](1-het[[k]]),(1-het[[k]])^2}->{{},{profiles[[k,L,1]]},{profiles[[k,L,2]]},profiles[[k,L]]}]],{L,NumberOfLoci}],{k,n}];Do[Do[AppendTo[mix[[L,r]],profiles[[i,L]]],{i,n}],{L,NumberOfLoci}];Do[mix[[L,r]]=Sort[DeleteDuplicates[Flatten[mix[[L,r]]]]],{L,NumberOfLoci}];profiles=MixtureDonors;,{r,rep}];profiles=MixtureDonors;{MixtureDonors,Transpose[mix]}]

(*  n : number of donors;
	rep: number of replicates;
	het: vector with heterozygous dropout probabilities;
	c: contamination parameters; 
 output is a table with first entry the table with donor profiles and second entry the replicates (alleles sorted) *)
RandomMixtureReplicates[n_,rep_,het_,c_]:=Module[{start=RandomMixtureNoContamination[n,rep,het], extra=Table[{},{L,NumberOfLoci}],mix=Table[{},{L,NumberOfLoci}]},
mix=start;
Do[Do[extra[[L]]=Table[RandomChoice[{c*Fr[[L,i]],1-c*Fr[[L,i]]}->{{i},{}}],{i,LL[[L]]}],{L,NumberOfLoci}];
Do[mix[[2,k,L]]=Sort[DeleteDuplicates[Flatten[AppendTo[mix[[2,k,L]],extra[[L]]]]]],{L,NumberOfLoci}],{k,rep}];mix]

(* create a mixture with n donors, in rep replicates given that their heterozygous dropout probabilities are in the vector het *)
(* the first of the n donors are in the second argument *)
(* no drop-in *)
RandomMixtureNoContaminationWithSomeDonorsFixed[n_,Donors_,rep_,het_]:=Module[{MixtureDonors=Table[{},{t,n}],profiles=Join[Donors,Table[Partition[RandomProfile[NumberOfLoci],2],{t,n-Length[Donors]}]],mix=Table[{},{L,NumberOfLoci},{r,rep}]},MixtureDonors=profiles;
Do[Do[Do[If[profiles[[k,L,1]]==profiles[[k,L,2]],profiles[[k,L]]=RandomChoice[{het[[k]]^2,1-het[[k]]^2}->{{},profiles[[k,L]]}],profiles[[k,L]]=RandomChoice[{het[[k]]^2,het[[k]](1-het[[k]]),het[[k]](1-het[[k]]),(1-het[[k]])^2}->{{},{profiles[[k,L,1]]},{profiles[[k,L,2]]},profiles[[k,L]]}]],{L,NumberOfLoci}],{k,n}];Do[Do[AppendTo[mix[[L,r]],profiles[[i,L]]],{i,n}],{L,NumberOfLoci}];Do[mix[[L,r]]=Sort[DeleteDuplicates[Flatten[mix[[L,r]]]]],{L,NumberOfLoci}];profiles=MixtureDonors;,{r,rep}];profiles=MixtureDonors;{MixtureDonors,Transpose[mix]}]

(* create a mixture with n donors, in rep replicates given that their heterozygous dropout probabilities are in the vector het *)
(* the first of the n donors are in the second argument *)
(* with drop-in *)
RandomMixtureReplicatesWithSomeDonorsFixed[n_,Donors_,rep_,het_,c_]:=Module[{start=RandomMixtureNoContaminationWithSomeDonorsFixed[n,Donors,rep,het], extra=Table[{},{L,NumberOfLoci}],mix=Table[{},{L,NumberOfLoci}]},
mix=start;
Do[Do[extra[[L]]=Table[RandomChoice[{c*Fr[[L,i]],1-c*Fr[[L,i]]}->{{i},{}}],{i,LL[[L]]}],{L,NumberOfLoci}];
Do[mix[[2,k,L]]=Sort[DeleteDuplicates[Flatten[AppendTo[mix[[2,k,L]],extra[[L]]]]]],{L,NumberOfLoci}],{k,rep}];mix]

Likelihood ratios, no relatedness

UseStutterModel=ChoiceDialog["Do you want to use the stutter model?",{"Yes"->"Y","No"->"N"}]

ParentOfMixAllele[L_,a_,MixtureAlleles_]:=Module[{StAllel,StAllelPos},
StAllel=Alleles[[L,MixtureAlleles[[L,a]]]]+1;
StAllelPos=If[MemberQ[Alleles[[L]],StAllel],Position[Alleles[[L]],StAllel][[1,1]],0];If[StAllelPos>0&& MemberQ[MixtureAlleles[[L]],StAllelPos],Position[MixtureAlleles[[L]],StAllelPos][[1,1]],0]]

(* the probability to see allele x from a person with genotype (a,b), no contamination *)AlleleProbOnePerson[x_,a_,b_,d_]:=Which[HomDropout=="Sq",If[a==b,If[a==x ,1-d^2,0],If[a==x||b==x,1-d,0]],HomDropout=="Tved",If[a==b,If[a==x ,1-Tved[d,beta],0],If[a==x||b==x,1-d,0]]]

(* the probability that an allele is visible in a mixture, given
L: the locus;
a: the allele;
GenotypesoOnLocus: the genotypes of all donors on locus L;
d: the heterozygous dropout probabilities of the donors;
c: dropin parameter *)
Which[UseStutterModel=="N",If[DropinModel=="Linear",
AlleleProb[L_,a_,GenotypesOnLocus_,d_,c_]:=1-(1-c*MixtureFr[[L,a]])*Product[1-AlleleProbOnePerson[a,GenotypesOnLocus[[i,1]],GenotypesOnLocus[[i,2]],d[[i]]],{i,Length[d]}],AlleleProb[L_,a_,GenotypesOnLocus_,d_,c_]:=1-Exp[-c*MixtureFr[[L,a]]]*Product[1-AlleleProbOnePerson[a,GenotypesOnLocus[[i,1]],GenotypesOnLocus[[i,2]],d[[i]]],{i,Length[d]}]],
UseStutterModel=="Y",
If[DropinModel=="Linear",
AlleleProb[L_,a_,GenotypesOnLocus_,d_,c_]:=1-(1-c*MixtureFr[[L,a]])*(1-StutterProb*(1-Product[1-AlleleProbOnePerson[StutterPositions[[L,a]],GenotypesOnLocus[[i,1]],GenotypesOnLocus[[i,2]],d[[i]]],{i,Length[d]}]))*Product[1-AlleleProbOnePerson[a,GenotypesOnLocus[[i,1]],GenotypesOnLocus[[i,2]],d[[i]]],{i,Length[d]}],AlleleProb[L_,a_,GenotypesOnLocus_,d_,c_]:=1-Exp[-c*MixtureFr[[L,a]]]*(1-StutterProb*(1-Product[1-AlleleProbOnePerson[StutterPositions[[L,a]],GenotypesOnLocus[[i,1]],GenotypesOnLocus[[i,2]],d[[i]]],{i,Length[d]}]))*Product[1-AlleleProbOnePerson[a,GenotypesOnLocus[[i,1]],GenotypesOnLocus[[i,2]],d[[i]]],{i,Length[d]}]]]

(* 1 replicate
MixtureOnLocus: the mixture on locus L;
L: the locus;
GenotypesoOnLocus: the genotypes of all donors on locus L;
d: the heterozygous dropout probabilities of the donors;
c: dropin parameter *)
SingleMixtureProbDonorsKnown[MixtureOnLocus_,L_,GenotypesOnLocus_,d_,c_]:=Product[AlleleProb[L,x,GenotypesOnLocus,d,c],{x,DeleteDuplicates[MixtureOnLocus]}]*Product[1-AlleleProb[L,x,GenotypesOnLocus,d,c],{x,Complement[Range[MixtureLL[[L]]],MixtureOnLocus]}]

(* likelihood of all replicates for
 MixturesOnAllLoci: the mixtures (all loci specified);
L: the locus;
Kgeno: genotypes of known donors;
n: total number of donors;
d: the heterozygous dropout probabilities of the donors;
c: dropin parameter *)
MixtureProb[MixturesOnAllLoci_,L_,Kgeno_,n_,d_,c_]:=Module[{u,x1,x2,x3,x4,x5,x6,x,p,mult,g,indices={},Ugeno={},AllGeno={}},
u=n-Length[Kgeno];
If[u>0,indices=Range[MixtureLL[[L]]*(MixtureLL[[L]]+1)/2]];
If[u>1,Do[indices[[t]]=Table[{indices[[t]],x},{x,indices[[t]],MixtureLL[[L]]*(MixtureLL[[L]]+1)/2}],{t,Length[indices]}];indices=Flatten[indices,1]];
If[u>2,Do[Do[indices[[t]]=Table[Join[indices[[t]],{x}],{x,Last[indices[[t]]],MixtureLL[[L]]*(MixtureLL[[L]]+1)/2}],{t,Length[indices]}];indices=Flatten[indices,1],{i,u-2}]];
g=DeleteDuplicates[Sort/@Tuples[Range[MixtureLL[[L]]],2]];
x=Take[{x1,x2,x3,x4,x5,x6},n];
x=Take[x,-u];
p=Quiet[Permutations[x]/.{x1->d[[1]],x2->d[[2]],x3->d[[3]],x4->d[[4]],x5->d[[5]],x6->d[[6]]}];
p=Join[Take[d,Length[Kgeno]],#]&/@p;
If[u>0,Ugeno=Flatten/@Table[g[[x]],{x,indices}]];
mult=If[u==0,{1},If[u>1,Table[Factorial[u]/Product[Factorial[x],{x,Tally[index][[All,2]]}],{index,indices}],Table[1,{Length[indices]}]]];
Do[AppendTo[AllGeno,Kgeno[[i,L]]],{i,Length[Kgeno]}];AllGeno=Flatten[AllGeno];
Do[If[Max[Length/@Kgeno[[i]]]>2||Min[Length/@Kgeno[[i]]]<2,Print["impossible genotypes"];Abort[]],{i,Length[Kgeno]}];
If[u>0,AllGeno=Join[AllGeno,#]&/@Ugeno];
If[u==0,AllGeno={Partition[AllGeno,2]},AllGeno=Partition[#,2]&/@AllGeno];
If[Length[d]≠n,Print["wrong length dropout vector! calculation aborted"];Abort[]];
Tabel=Table[Product[SingleMixtureProbDonorsKnown[MixturesOnAllLoci[[t,L]],L,AllGeno[[i]],Take[{x1,x2,x3,x4,x5,x6},n],c],{t,Length[MixturesOnAllLoci]}],{i,Length[AllGeno]}];
Sum[Sum[Quiet[Tabel[[i]]/.{x1->t[[1]],x2->t[[2]],x3->t[[3]],x4->t[[4]],x5->t[[5]],x6->t[[6]]}],{t,p}]*10^Total[Log[10,MixtureFr[[L,Take[Flatten[AllGeno[[i]]],-2*u]]]]]*mult[[i]]/Factorial[u]*10^Total[Log[10,2-KD/@Take[AllGeno[[i]],-u]]],{i,Length[AllGeno]}]]

(* likelihood of all replicates for
 MixturesOnAllLoci: the mixtures (all loci specified);
L: the locus;
Kgeno: genotypes of known donors;
n: total number of donors;
d: the heterozygous dropout probability, same for all of the donors;
c: dropin parameter *)
MixtureProbEqualDropout[MixturesOnAllLoci_,L_,Kgeno_,n_,d_,c_]:=Module[{u,x1,x2,x3,x4,x5,x6,x,p,mult,g,indices={},Ugeno={},AllGeno={}},
u=n-Length[Kgeno];
If[u>0,indices=Range[MixtureLL[[L]]*(MixtureLL[[L]]+1)/2]];
If[u>1,Do[indices[[t]]=Table[{indices[[t]],x},{x,indices[[t]],MixtureLL[[L]]*(MixtureLL[[L]]+1)/2}],{t,Length[indices]}];indices=Flatten[indices,1]];
If[u>2,Do[Do[indices[[t]]=Table[Join[indices[[t]],{x}],{x,Last[indices[[t]]],MixtureLL[[L]]*(MixtureLL[[L]]+1)/2}],{t,Length[indices]}];indices=Flatten[indices,1],{i,u-2}]];
g=DeleteDuplicates[Sort/@Tuples[Range[MixtureLL[[L]]],2]];
x=Take[{x1,x2,x3,x4,x5,x6},n];
x=Take[x,-u];
If[u>0,Ugeno=Flatten/@Table[g[[x]],{x,indices}]];
mult=If[u==0,{1},If[u>1,Table[Factorial[u]/Product[Factorial[x],{x,Tally[index][[All,2]]}],{index,indices}],Table[1,{Length[indices]}]]];
Do[AppendTo[AllGeno,Kgeno[[i,L]]],{i,Length[Kgeno]}];AllGeno=Flatten[AllGeno];
Do[If[Max[Length/@Kgeno[[i]]]>2||Min[Length/@Kgeno[[i]]]<2,Print["impossible genotypes"];Abort[]],{i,Length[Kgeno]}];
If[u>0,AllGeno=Join[AllGeno,#]&/@Ugeno];
If[u==0,AllGeno={Partition[AllGeno,2]},AllGeno=Partition[#,2]&/@AllGeno];
Tabel=Table[Product[SingleMixtureProbDonorsKnown[MixturesOnAllLoci[[t,L]],L,AllGeno[[i]],Take[{d,d,d,d,d,d},n],c],{t,Length[MixturesOnAllLoci]}],{i,Length[AllGeno]}];
Sum[
Tabel[[i]]*10^Total[Log[10,MixtureFr[[L,Take[Flatten[AllGeno[[i]]],-2*u]]]]]*mult[[i]]*10^Total[Log[10,2-KD/@Take[AllGeno[[i]],-u]]],{i,Length[AllGeno]}]]

(* LR assuming equal dropout; slightly faster implementation of the same model.*)
(* d is now a single number taken as probability of dropout for all donors *)
LREqualDropout[MixturesOnAllLoci_,HpDonors_,HdDonors_,n_,d_,c_,kit_,monitor_]:=Module[{MixtureAlleles,start=SessionTime[],sortkit},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
If[monitor≠"Y" && monitor ≠ "N",Print["Incorrect choice for monitor value, aborted"];Abort[]];
If[Length[HpDonors]>n && Length[HdDonors]>n,Print["There are more hypothesized contributors than donors. Computation cancelled"]];
If[c<0 || (c>1 && DropinModel=="Linear"), Print["inadmissible choice of dropin parameter. Computation aborted"];Abort[]];
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];MixtureLL=Length/@MixtureFr ;
sortkit=Sort[kit,MixtureLL[[#1]]<MixtureLL[[#2]]&];
res=Table[{Loci[[sortkit[[L]]]],MixtureLL[[sortkit[[L]]]]-1,1,1,0,"Queued"},{L,Length[kit]}];PrependTo[res,{"Locus","distinct alleles", "LR locus","cumulative LR","Log(10,cumu-LR)", "time for locus (s)"}];
AllelesHpDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HpDonors[[i,L,1]],HpDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HpDonors]}];
Do[Do[If[Length[AllelesHpDonors[[i,L]]]==0,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHpDonors[[i,L]]]==1,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HpDonors]}];
AllelesHdDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HdDonors[[i,L,1]],HdDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HdDonors]}];
Do[Do[If[Length[AllelesHdDonors[[i,L]]]==0,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHdDonors[[i,L]]]==1,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HdDonors]}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
Do[start=SessionTime[];res[[L+1,6]]="Processing";res[[L+1,3]]=MixtureProbEqualDropout[
(* de replica's, alle loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes bekende donoren Hp*)AllelesHpDonors,
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropin parameter *)c]/MixtureProbEqualDropout[
(* de replica's, alle loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes bekende donoren Hd *)AllelesHdDonors,
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropin parameter *)c];
If[L>1,res[[L+1,4]]=res[[L,4]]*res[[L+1,3]];res[[L+1,5]]=Log[10,res[[L+1,4]]],res[[L+1,4]]=res[[L+1,3]];res[[L+1,5]]=Log[10,res[[L+1,4]]]];
res[[L+1,6]]=SessionTime[]-start,{L,Length[kit]}];
res[[Length[kit]+1,4]]=res[[Length[kit]+1,3]]res[[Length[kit],4]];If[(NumberQ[d]=={True}||NumberQ[d]==True)&&monitor=="Y",Print[res//MatrixForm]];Table[res[[k,3]],{k,2,Length[res]}]
]

(* the output contains three vectors *)
(* output[[1]] contains the likelihoods under Hp *)
(* output[[2]] contains the likelihoods under Hd *)
(* output[[3]] contains the likelihood ratios *)
(* the order of the loci is not changed *)
(* see also LR definition *)
LRinclLikelihoodsNoThetaFixedOrder[MixturesOnAllLoci_,HpDonors_,HdDonors_,n_,d_,c_,kit_,monitor_]:=Module[{MixtureAlleles,start=SessionTime[]},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
If[monitor≠"Y" && monitor ≠ "N",Print["Incorrect choice for monitor value, aborted"];Abort[]];
If[Length[HpDonors]>n && Length[HdDonors]>n,Print["There are more hypothesized contributors than donors. Computation cancelled"]];
If[c<0 || (c>1 && DropinModel=="Linear"), Print["inadmissible choice of dropin parameter. Computation aborted"];Abort[]];
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];MixtureLL=Length/@MixtureFr ;
res=Table[{Loci[[L]],MixtureLL[[L]]-1,1,1,1,1,1,"Queued"},{L,Length[kit]}];PrependTo[res,{"Locus","distinct alleles", "Lh Hp","Lh Hd","LR locus","cumulative LR","Log(10,cumu-LR)", "time for locus (s)"}];
AllelesHpDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HpDonors[[i,L,1]],HpDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HpDonors]}];
Do[Do[If[Length[AllelesHpDonors[[i,L]]]==0,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHpDonors[[i,L]]]==1,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HpDonors]}];
AllelesHdDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HdDonors[[i,L,1]],HdDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HdDonors]}];
Do[Do[If[Length[AllelesHdDonors[[i,L]]]==0,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHdDonors[[i,L]]]==1,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HdDonors]}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
Do[start=SessionTime[];res[[L+1,3]]=MixtureProb[
(* de replica's, alle loci *)AllReplicates,
(* locus *)kit[[L]],
(* genotypes bekende donoren Hp*)AllelesHpDonors,
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropin parameter *)c];
res[[L+1,4]]=MixtureProb[
(* de replica's, alle loci *)AllReplicates,
(* locus *)kit[[L]],
(* genotypes bekende donoren Hd *)AllelesHdDonors,
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropin parameter *)c];
res[[L+1,5]]=res[[L+1,3]]/res[[L+1,4]];
If[L>1,res[[L+1,6]]=res[[L,6]]*res[[L+1,5]];res[[L+1,7]]=Log[10,res[[L+1,6]]],res[[L+1,6]]=res[[L+1,5]];res[[L+1,7]]=Log[10,res[[L+1,6]]]];
res[[L+1,8]]=SessionTime[]-start,{L,Length[kit]}];
res[[Length[kit]+1,6]]=res[[Length[kit]+1,5]]res[[Length[kit],6]];If[DeleteDuplicates[NumberQ/@d]=={True}&&monitor=="Y",Print[res//MatrixForm]];{Table[res[[k,3]],{k,2,Length[res]}],Table[res[[k,4]],{k,2,Length[res]}],Table[res[[k,5]],{k,2,Length[res]}]}
]

(* the output contains three vectors *)
(* output[[1]] contains the likelihoods under Hp *)
(* output[[2]] contains the likelihoods under Hd *)
(* output[[3]] contains the likelihood ratios *)
(* the order of the loci is obtained by sorting from fewer to more recorded alleles *)
(* see also LR definition *)
LRinclLikelihoodsNoTheta[MixturesOnAllLoci_,HpDonors_,HdDonors_,n_,d_,c_,kit_,monitor_]:=Module[{MixtureAlleles,start=SessionTime[],sortkit},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
If[monitor≠"Y" && monitor ≠ "N",Print["Incorrect choice for monitor value, aborted"];Abort[]];
If[Length[HpDonors]>n && Length[HdDonors]>n,Print["There are more hypothesized contributors than donors. Computation cancelled"]];
If[c<0 || (c>1 && DropinModel=="Linear"), Print["inadmissible choice of dropin parameter. Computation aborted"];Abort[]];
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];MixtureLL=Length/@MixtureFr ;
sortkit=Sort[kit,MixtureLL[[#1]]<MixtureLL[[#2]]&];
res=Table[{Loci[[sortkit[[L]]]],MixtureLL[[sortkit[[L]]]]-1,1,1,1,1,1,"Queued"},{L,Length[kit]}];PrependTo[res,{"Locus","distinct alleles", "Lh Hp","Lh Hd","LR locus","cumulative LR","Log(10,cumu-LR)", "time for locus (s)"}];
AllelesHpDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HpDonors[[i,L,1]],HpDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HpDonors]}];
Do[Do[If[Length[AllelesHpDonors[[i,L]]]==0,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHpDonors[[i,L]]]==1,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HpDonors]}];
AllelesHdDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HdDonors[[i,L,1]],HdDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HdDonors]}];
Do[Do[If[Length[AllelesHdDonors[[i,L]]]==0,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHdDonors[[i,L]]]==1,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HdDonors]}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
Do[start=SessionTime[];res[[L+1,3]]=MixtureProb[
(* de replica's, alle loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes bekende donoren Hp*)AllelesHpDonors,
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropin parameter *)c];
res[[L+1,4]]=MixtureProb[
(* de replica's, alle loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes bekende donoren Hd *)AllelesHdDonors,
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropin parameter *)c];
res[[L+1,5]]=res[[L+1,3]]/res[[L+1,4]];
If[L>1,res[[L+1,6]]=res[[L,6]]*res[[L+1,5]];res[[L+1,7]]=Log[10,res[[L+1,6]]],res[[L+1,6]]=res[[L+1,5]];res[[L+1,7]]=Log[10,res[[L+1,6]]]];
res[[L+1,8]]=SessionTime[]-start,{L,Length[kit]}];
res[[Length[kit]+1,6]]=res[[Length[kit]+1,5]]res[[Length[kit],6]];If[DeleteDuplicates[NumberQ/@d]=={True}&&monitor=="Y",Print[res//MatrixForm]];{Table[res[[k,3]],{k,2,Length[res]}],Table[res[[k,4]],{k,2,Length[res]}],Table[res[[k,5]],{k,2,Length[res]}]}
]

(* gives the LR for the *)
(*MixturesOnAllLoci: all replicates *)
(* HpDonors: the donors present under Hp *)
(* HdDonors: the donors present under Hd *)
(* n: total number of contributors *)
(* d: vector of heterozygous dropout probabilities *)
(* c: dropin parameter *)
(* monitor: if this is filled in as "Y" then the output contains a table in matrixform with the LRs per locus *)
(* intermediate results can still be monitored otherwise *)
LRArbitraryDropoutNoTheta[MixturesOnAllLoci_,HpDonors_,HdDonors_,n_,d_,c_,kit_,monitor_]:=Module[{MixtureAlleles,start=SessionTime[],sortkit},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
If[monitor≠"Y" && monitor ≠ "N",Print["Incorrect choice for monitor value, aborted"];Abort[]];
If[Length[HpDonors]>n && Length[HdDonors]>n,Print["There are more hypothesized contributors than donors. Computation cancelled"]];
If[c<0 || (c>1 && DropinModel=="Linear"), Print["inadmissible choice of dropin parameter. Computation aborted"];Abort[]];
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];MixtureLL=Length/@MixtureFr ;
sortkit=Sort[kit,MixtureLL[[#1]]<MixtureLL[[#2]]&];
res=Table[{Loci[[sortkit[[L]]]],MixtureLL[[sortkit[[L]]]]-1,1,1,0,"Queued"},{L,Length[kit]}];PrependTo[res,{"Locus","distinct alleles", "LR locus","cumulative LR","Log(10,cumu-LR)", "time for locus (s)"}];
AllelesHpDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HpDonors[[i,L,1]],HpDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HpDonors]}];
Do[Do[If[Length[AllelesHpDonors[[i,L]]]==0,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHpDonors[[i,L]]]==1,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HpDonors]}];
AllelesHdDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HdDonors[[i,L,1]],HdDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HdDonors]}];
Do[Do[If[Length[AllelesHdDonors[[i,L]]]==0,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHdDonors[[i,L]]]==1,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HdDonors]}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
Do[start=SessionTime[];res[[L+1,6]]="Processing";res[[L+1,3]]=MixtureProb[
(* de replica's, alle loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes bekende donoren Hp*)AllelesHpDonors,
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropin parameter *)c]/MixtureProb[
(* de replica's, alle loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes bekende donoren Hd *)AllelesHdDonors,
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropin parameter *)c];
If[L>1,res[[L+1,4]]=res[[L,4]]*res[[L+1,3]];res[[L+1,5]]=Log[10,res[[L+1,4]]],res[[L+1,4]]=res[[L+1,3]];res[[L+1,5]]=Log[10,res[[L+1,4]]]];
res[[L+1,6]]=SessionTime[]-start,{L,Length[kit]}];
res[[Length[kit]+1,4]]=res[[Length[kit]+1,3]]res[[Length[kit],4]];If[DeleteDuplicates[NumberQ/@d]=={True}&&monitor=="Y",Print[res//MatrixForm]];Table[res[[k,3]],{k,2,Length[res]}]
]

(* same as LR but the same probability of dropout applied to all contributors *)
(* d is now a single number and equals the used probability of dropout *)
(* loci are not sorted according to number of recorded alleles but left in the standard order *)
LREqualDropoutFixedOrder[MixturesOnAllLoci_,HpDonors_,HdDonors_,n_,d_,c_,kit_,monitor_]:=Module[{MixtureAlleles,start=SessionTime[]},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
If[monitor≠"Y" && monitor ≠ "N",Print["Incorrect choice for monitor value, aborted"];Abort[]];
If[Length[HpDonors]>n && Length[HdDonors]>n,Print["There are more hypothesized contributors than donors. Computation cancelled"]];
If[c<0 || (c>1 && DropinModel=="Linear"), Print["inadmissible choice of dropin parameter. Computation aborted"];Abort[]];
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];MixtureLL=Length/@MixtureFr ;
res=Table[{Loci[[kit[[L]]]],MixtureLL[[kit[[L]]]]-1,1,1,0,"Queued"},{L,Length[kit]}];PrependTo[res,{"Locus","distinct alleles", "LR locus","cumulative LR","Log(10,cumu-LR)", "time for locus (s)"}];
AllelesHpDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HpDonors[[i,L,1]],HpDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HpDonors]}];
Do[Do[If[Length[AllelesHpDonors[[i,L]]]==0,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHpDonors[[i,L]]]==1,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HpDonors]}];
AllelesHdDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HdDonors[[i,L,1]],HdDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HdDonors]}];
Do[Do[If[Length[AllelesHdDonors[[i,L]]]==0,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHdDonors[[i,L]]]==1,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HdDonors]}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
Do[start=SessionTime[];res[[L+1,6]]="Processing";res[[L+1,3]]=MixtureProbEqualDropout[
(* de replica's, alle loci *)AllReplicates,
(* locus *)kit[[L]],
(* genotypes bekende donoren Hp*)AllelesHpDonors,
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropin parameter *)c]/MixtureProbEqualDropout[
(* de replica's, alle loci *)AllReplicates,
(* locus *)kit[[L]],
(* genotypes bekende donoren Hd *)AllelesHdDonors,
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropin parameter *)c];
If[L>1,res[[L+1,4]]=res[[L,4]]*res[[L+1,3]];res[[L+1,5]]=Log[10,res[[L+1,4]]],res[[L+1,4]]=res[[L+1,3]];res[[L+1,5]]=Log[10,res[[L+1,4]]]];
res[[L+1,6]]=SessionTime[]-start,{L,Length[kit]}];
res[[Length[kit]+1,4]]=res[[Length[kit]+1,3]]res[[Length[kit],4]];If[NumberQ[d]=={True}&&monitor=="Y",Print[res//MatrixForm]];Table[res[[k,3]],{k,2,Length[res]}]
]

(* compute a LR with different numbers of contributors *)
(* H1: the mixture has n1 donors with probabilities of dropout in the vector d1 and drop-in c1. The first donors are in HpDonors *)
(* H2: the mixture has n2 donors with probabilities of dropout in the vector d2 and drop-in c2. The first donors are in HdDonors *)
LRForNumberOfContributors[MixturesOnAllLoci_,HpDonors_,HdDonors_,n1_,d1_,c1_,n2_,d2_,c2_,kit_,monitor_]:=Module[{MixtureAlleles,start=SessionTime[],sortkit},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
If[monitor≠"Y" && monitor ≠ "N",Print["Incorrect choice for monitor value, aborted"];Abort[]];
If[Length[HpDonors]>n && Length[HdDonors]>n,Print["There are more hypothesized contributors than donors. Computation cancelled"]];
If[c<0 || (c>1 && DropinModel=="Linear"), Print["inadmissible choice of dropin parameter. Computation aborted"];Abort[]];                  
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];MixtureLL=Length/@MixtureFr ;
sortkit=Sort[kit,MixtureLL[[#1]]<MixtureLL[[#2]]&];
res=Table[{Loci[[sortkit[[L]]]],MixtureLL[[sortkit[[L]]]]-1,1,1,0,"Queued"},{L,Length[kit]}];PrependTo[res,{"Locus","distinct alleles", "LR locus","cumulative LR","Log(10,cumu-LR)", "time for locus (s)"}];
AllelesHpDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HpDonors[[i,L,1]],HpDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HpDonors]}];
Do[Do[If[Length[AllelesHpDonors[[i,L]]]==0,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHpDonors[[i,L]]]==1,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HpDonors]}];
AllelesHdDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HdDonors[[i,L,1]],HdDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HdDonors]}];
Do[Do[If[Length[AllelesHdDonors[[i,L]]]==0,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHdDonors[[i,L]]]==1,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HdDonors]}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
Do[start=SessionTime[];res[[L+1,6]]="Processing";res[[L+1,3]]=MixtureProb[
(* de replica's, alle loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes bekende donoren Hp*)AllelesHpDonors,
(* totaal aantal donoren *)n1,
(* dropoutkans *)d1,
(* dropin parameter *)c1]/MixtureProb[
(* de replica's, alle loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes bekende donoren Hd *)AllelesHdDonors,
(* totaal aantal donoren *)n2,
(* dropoutkans *)d2,
(* dropin parameter *)c2];
If[L>1,res[[L+1,4]]=res[[L,4]]*res[[L+1,3]];res[[L+1,5]]=Log[10,res[[L+1,4]]],res[[L+1,4]]=res[[L+1,3]];res[[L+1,5]]=Log[10,res[[L+1,4]]]];
res[[L+1,6]]=SessionTime[]-start,{L,Length[kit]}];
res[[Length[kit]+1,4]]=res[[Length[kit]+1,3]]res[[Length[kit],4]];If[DeleteDuplicates[NumberQ/@d1]=={True}&&monitor=="Y",Print[res//MatrixForm]];Table[res[[k,3]],{k,2,Length[res]}]
]

(* same as LRForNumberOfContributors but now equal dropout *)
LREqualDropoutForNumberOfContributors[MixturesOnAllLoci_,HpDonors_,HdDonors_,n1_,d1_,c1_,n2_,d2_,c2_,kit_,monitor_]:=Module[{MixtureAlleles,start=SessionTime[],sortkit},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
If[monitor≠"Y" && monitor ≠ "N",Print["Incorrect choice for monitor value, aborted"];Abort[]];
If[Length[HpDonors]>n && Length[HdDonors]>n,Print["There are more hypothesized contributors than donors. Computation cancelled"]];
If[c<0 || (c>1 && DropinModel=="Linear"), Print["inadmissible choice of dropin parameter. Computation aborted"];Abort[]];                  
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];MixtureLL=Length/@MixtureFr ;
sortkit=Sort[kit,MixtureLL[[#1]]<MixtureLL[[#2]]&];
res=Table[{Loci[[sortkit[[L]]]],MixtureLL[[sortkit[[L]]]]-1,1,1,0,"Queued"},{L,Length[kit]}];PrependTo[res,{"Locus","distinct alleles", "LR locus","cumulative LR","Log(10,cumu-LR)", "time for locus (s)"}];
AllelesHpDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HpDonors[[i,L,1]],HpDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HpDonors]}];
Do[Do[If[Length[AllelesHpDonors[[i,L]]]==0,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHpDonors[[i,L]]]==1,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HpDonors]}];
AllelesHdDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HdDonors[[i,L,1]],HdDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HdDonors]}];
Do[Do[If[Length[AllelesHdDonors[[i,L]]]==0,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHdDonors[[i,L]]]==1,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HdDonors]}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
Do[start=SessionTime[];res[[L+1,6]]="Processing";res[[L+1,3]]=MixtureProbEqualDropout[
(* de replica's, alle loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes bekende donoren Hp*)AllelesHpDonors,
(* totaal aantal donoren *)n1,
(* dropoutkans *)d1,
(* dropin parameter *)c1]/MixtureProbEqualDropout[
(* de replica's, alle loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes bekende donoren Hd *)AllelesHdDonors,
(* totaal aantal donoren *)n2,
(* dropoutkans *)d2,
(* dropin parameter *)c2];
If[L>1,res[[L+1,4]]=res[[L,4]]*res[[L+1,3]];res[[L+1,5]]=Log[10,res[[L+1,4]]],res[[L+1,4]]=res[[L+1,3]];res[[L+1,5]]=Log[10,res[[L+1,4]]]];
res[[L+1,6]]=SessionTime[]-start,{L,Length[kit]}];
res[[Length[kit]+1,4]]=res[[Length[kit]+1,3]]res[[Length[kit],4]];If[monitor=="Y",Print[res//MatrixForm]];Table[res[[k,3]],{k,2,Length[res]}]
]

Likelihood ratios definitions, no relatedness, with theta

ProbThetaSeveralAlleles[L_,a_,θ_]:=Module[{m=Table[Count[a[[Range[k-1]]],a[[k]]],{k,Length[a]}]},Product[(m[[k]]*θ+MixtureFr[[L,a[[k]]]]*(1-θ))/(1+(k-2)θ),{k,Length[a]}]]

(* werkt t/m 5p 8a *)
MixtureProbTheta[MixturesOnAllLoci_,L_,Kgeno_,Condgeno_,n_,d_,c_,θ_]:=Module[{u,x1,x2,x3,x4,x5,x6,x,p,e,k,mult,indices={},Ugeno={},AlleGeno={}},
k=Length[Kgeno];
u=n-Length[Kgeno];
e=Length[Condgeno];
If[u>0,indices=Range[MixtureLL[[L]]*(MixtureLL[[L]]+1)/2]];
If[u>1,Do[indices[[t]]=Table[{indices[[t]],x},{x,indices[[t]],MixtureLL[[L]]*(MixtureLL[[L]]+1)/2}],{t,Length[indices]}];indices=Flatten[indices,1]];
If[u>2,Do[Do[indices[[t]]=Table[Join[indices[[t]],{x}],{x,Last[indices[[t]]],MixtureLL[[L]]*(MixtureLL[[L]]+1)/2}],{t,Length[indices]}];indices=Flatten[indices,1],{i,u-2}]];
g=DeleteDuplicates[Sort/@Tuples[Range[MixtureLL[[L]]],2]];
x=Take[{x1,x2,x3,x4,x5,x6},n];
x=Take[x,-u];
p=Quiet[Permutations[x]/.{x1->d[[1]],x2->d[[2]],x3->d[[3]],x4->d[[4]],x5->d[[5]],x6->d[[6]]}];
p=Join[Take[d,Length[Kgeno]],#]&/@p;
If[u>0,Ugeno=Flatten/@Table[g[[x]],{x,indices}]];
mult=If[u==0,{1},If[u>1,Table[Factorial[u]/Product[Factorial[x],{x,Tally[index][[All,2]]}],{index,indices}],Table[1,{Length[indices]}]]];
Do[AppendTo[AlleGeno,Condgeno[[i,L]]],{i,Length[Condgeno]}];
Do[AppendTo[AlleGeno,Kgeno[[i,L]]],{i,Length[Kgeno]}];
Do[If[Max[Length/@Kgeno[[i]]]>2||Min[Length/@Kgeno[[i]]]<2,Print["verkeerde genotypes, match zonder statistische onderbouwing"];Abort[]],{i,Length[Kgeno]}];
AlleGeno=Flatten[AlleGeno];
If[u>0,AlleGeno=Join[AlleGeno,#]&/@Ugeno];
If[u==0,AlleGeno={Partition[AlleGeno,2]},AlleGeno=Partition[#,2]&/@AlleGeno];
If[Length[d]≠n,Print["verkeerde lengte dropout vector! daar begin ik niet aan"];Abort[]];
Tabel=Table[Product[SingleMixtureProbDonorsKnown[MixturesOnAllLoci[[t,L]],L,Take[AlleGeno[[i]],-(u+k)],Take[{x1,x2,x3,x4,x5,x6},n],c],{t,Length[MixturesOnAllLoci]}],{i,Length[AlleGeno]}];
Sum[Sum[Quiet[Tabel[[i]]/.{x1->t[[1]],x2->t[[2]],x3->t[[3]],x4->t[[4]],x5->t[[5]],x6->t[[6]]}],{t,p}]*ProbThetaSeveralAlleles[L,Flatten[AlleGeno[[i]]],θ]*mult[[i]]/Factorial[u]*10^Total[Log[10,2-KD/@Take[AlleGeno[[i]],-u]]],{i,Length[AlleGeno]}]/ProbThetaSeveralAlleles[L,Take[Flatten[AlleGeno[[1]]],2*(k+e)],θ]]

(* 1 replica
mengsel_ het geziene mengsel
L_ het locus
genotypes_ van alle donoren
d de (heterozygote) dropoutkansen
c de dropinkans *)
MixtureProbEqualDropoutTheta[MixturesOnAllLoci_,L_,Kgeno_,Condgeno_,n_,d_,c_,θ_]:=Module[{u,k,e,x1,x2,x3,x4,x5,x6,x,p,mult,g,indices,Ugeno={},AlleGeno={}},
k=Length[Kgeno];
u=n-Length[Kgeno];
e=Length[Condgeno];
If[u>0,indices=Range[MixtureLL[[L]]*(MixtureLL[[L]]+1)/2]];
If[u>1,Do[indices[[t]]=Table[{indices[[t]],x},{x,indices[[t]],MixtureLL[[L]]*(MixtureLL[[L]]+1)/2}],{t,Length[indices]}];indices=Flatten[indices,1]];
If[u>2,Do[Do[indices[[t]]=Table[Join[indices[[t]],{x}],{x,Last[indices[[t]]],MixtureLL[[L]]*(MixtureLL[[L]]+1)/2}],{t,Length[indices]}];indices=Flatten[indices,1],{i,u-2}]];
g=DeleteDuplicates[Sort/@Tuples[Range[MixtureLL[[L]]],2]];
x=Take[{x1,x2,x3,x4,x5,x6},n];
x=Take[x,-u];
If[u>0,Ugeno=Flatten/@Table[g[[x]],{x,indices}]];
mult=If[u==0,{1},If[u>1,Table[Factorial[u]/Product[Factorial[x],{x,Tally[index][[All,2]]}],{index,indices}],Table[1,{Length[indices]}]]];
Do[AppendTo[AlleGeno,Condgeno[[i,L]]],{i,Length[Condgeno]}];
Do[AppendTo[AlleGeno,Kgeno[[i,L]]],{i,Length[Kgeno]}];
AlleGeno=Flatten[AlleGeno];
Do[If[Max[Length/@Kgeno[[i]]]>2||Min[Length/@Kgeno[[i]]]<2,Print["verkeerde genotypes, match zonder statistische onderbouwing"];Abort[]],{i,Length[Kgeno]}];
If[u>0,AlleGeno=Join[AlleGeno,#]&/@Ugeno];
If[u==0,AlleGeno={Partition[AlleGeno,2]},AlleGeno=Partition[#,2]&/@AlleGeno];
Tabel=Table[Product[SingleMixtureProbDonorsKnown[MixturesOnAllLoci[[t,L]],L,Take[AlleGeno[[i]],-n],Take[{d,d,d,d,d,d},n],c],{t,Length[MixturesOnAllLoci]}],{i,Length[AlleGeno]}];
Sum[
Tabel[[i]]*ProbThetaSeveralAlleles[L,Flatten[AlleGeno[[i]]],θ]*mult[[i]]*10^Total[Log[10,2-KD/@Take[AlleGeno[[i]],-u]]],{i,Length[AlleGeno]}]/ProbThetaSeveralAlleles[L,Take[Flatten[AlleGeno[[1]]],2*(k+e)],θ]]

(* gives the LR for the *)
(*MixturesOnAllLoci: all replicates *)
(* HpDonors: the donors present under Hp *)
(* HdDonors: the donors present under Hd *)
(* n: total number of contributors *)
(* d: vector of heterozygous dropout probabilities *)
(* c: dropin parameter *)
(* monitor: if this is filled in as "Y" then the output contains a table in matrixform with the LRs per locus *)
(* intermediate results can still be monitored otherwise *)
LRArbitraryDropoutTheta[MixturesOnAllLoci_,HpDonors_,HdDonors_,n_,d_,c_,kit_,θ_,monitor_]:=Module[{MixtureAlleles,start=SessionTime[],sortkit},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
If[monitor≠"Y" && monitor ≠ "N",Print["Incorrect choice for monitor value, aborted"];Abort[]];
If[Length[HpDonors]>n && Length[HdDonors]>n,Print["There are more hypothesized contributors than donors. Computation cancelled"]];
If[c<0 || (c>1 && DropinModel=="Linear"), Print["inadmissible choice of dropin parameter. Computation aborted"];Abort[]];
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];MixtureLL=Length/@MixtureFr ;
sortkit=Sort[kit,MixtureLL[[#1]]<MixtureLL[[#2]]&];
res=Table[{Loci[[sortkit[[L]]]],MixtureLL[[sortkit[[L]]]]-1,1,1,0,"Queued"},{L,Length[kit]}];PrependTo[res,{"Locus","distinct alleles", "LR locus","cumulative LR","Log(10,cumu-LR)", "time for locus (s)"}];
AllelesHpDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HpDonors[[i,L,1]],HpDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HpDonors]}];
Do[Do[If[Length[AllelesHpDonors[[i,L]]]==0,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHpDonors[[i,L]]]==1,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HpDonors]}];
AllelesHdDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HdDonors[[i,L,1]],HdDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HdDonors]}];
Do[Do[If[Length[AllelesHdDonors[[i,L]]]==0,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHdDonors[[i,L]]]==1,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HdDonors]}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
Do[start=SessionTime[];res[[L+1,6]]="Processing";res[[L+1,3]]=MixtureProbTheta[
(* de replica's, alle loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes bekende donoren Hp*)AllelesHpDonors,
(* voor theta *)Complement[AllelesHdDonors,AllelesHpDonors],
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropinkans *)c,
(* theta *) θ]/MixtureProbTheta[
(* de replica's, alle loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes bekende donoren Hd *)AllelesHdDonors,
(* extra voor theta *)Complement[AllelesHpDonors,AllelesHdDonors],
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropinkans *)c,
(* theta *) θ];
If[L>1,res[[L+1,4]]=res[[L,4]]*res[[L+1,3]];res[[L+1,5]]=Log[10,res[[L+1,4]]],res[[L+1,4]]=res[[L+1,3]];res[[L+1,5]]=Log[10,res[[L+1,4]]]];
res[[L+1,6]]=SessionTime[]-start,{L,Length[kit]}];
res[[Length[kit]+1,4]]=res[[Length[kit]+1,3]]res[[Length[kit],4]];If[DeleteDuplicates[NumberQ/@d]=={True}&&monitor=="Y",Print[res//MatrixForm]];Table[res[[k,3]],{k,2,Length[res]}]
]

(* LR assuming equal dropout; slightly faster implementation of the same model.*)
(* d is now a single number taken as probability of dropout for all donors *)
LREqualDropoutTheta[MixturesOnAllLoci_,HpDonors_,HdDonors_,n_,d_,c_,kit_,θ_,monitor_]:=Module[{MixtureAlleles,start=SessionTime[],sortkit},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
If[monitor≠"Y" && monitor ≠ "N",Print["Incorrect choice for monitor value, aborted"];Abort[]];
If[Length[HpDonors]>n && Length[HdDonors]>n,Print["There are more hypothesized contributors than donors. Computation cancelled"]];
If[c<0 || (c>1 && DropinModel=="Linear"), Print["inadmissible choice of dropin parameter. Computation aborted"];Abort[]];
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];MixtureLL=Length/@MixtureFr ;
sortkit=Sort[kit,MixtureLL[[#1]]<MixtureLL[[#2]]&];
res=Table[{Loci[[sortkit[[L]]]],MixtureLL[[sortkit[[L]]]]-1,1,1,0,"Queued"},{L,Length[kit]}];PrependTo[res,{"Locus","distinct alleles", "LR locus","cumulative LR","Log(10,cumu-LR)", "time for locus (s)"}];
AllelesHpDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HpDonors[[i,L,1]],HpDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HpDonors]}];
Do[Do[If[Length[AllelesHpDonors[[i,L]]]==0,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHpDonors[[i,L]]]==1,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HpDonors]}];
AllelesHdDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HdDonors[[i,L,1]],HdDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HdDonors]}];
Do[Do[If[Length[AllelesHdDonors[[i,L]]]==0,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHdDonors[[i,L]]]==1,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HdDonors]}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
Do[start=SessionTime[];res[[L+1,6]]="Processing";res[[L+1,3]]=MixtureProbEqualDropoutTheta[
(* de replica's, alle loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes bekende donoren Hp*)AllelesHpDonors,
(* voor theta *)Complement[AllelesHdDonors,AllelesHpDonors],
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropinkans *)c,
(* theta *)θ]/MixtureProbEqualDropoutTheta[
(* de replica's, alle loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes bekende donoren Hd *)AllelesHdDonors,
(* extra voor theta *)Complement[AllelesHpDonors,AllelesHdDonors],
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropinkans *)c,
(* theta *) θ];
If[L>1,res[[L+1,4]]=res[[L,4]]*res[[L+1,3]];res[[L+1,5]]=Log[10,res[[L+1,4]]],res[[L+1,4]]=res[[L+1,3]];res[[L+1,5]]=Log[10,res[[L+1,4]]]];
res[[L+1,6]]=SessionTime[]-start,{L,Length[kit]}];
res[[Length[kit]+1,4]]=res[[Length[kit]+1,3]]res[[Length[kit],4]];If[(NumberQ[d]=={True}||NumberQ[d]==True)&&monitor=="Y",Print[res//MatrixForm]];Table[res[[k,3]],{k,2,Length[res]}]
]

(* the output contains three vectors *)
(* output[[1]] contains the likelihoods under Hp *)
(* output[[2]] contains the likelihoods under Hd *)
(* output[[3]] contains the likelihood ratios *)
(* the order of the loci is obtained by sorting from fewer to more recorded alleles *)
(* see also LR definition *)
LRinclLikelihoodsEqualDropoutTheta[MixturesOnAllLoci_,HpDonors_,HdDonors_,n_,d_,c_,kit_,θ_,monitor_]:=Module[{MixtureAlleles,start=SessionTime[],sortkit},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
If[monitor≠"Y" && monitor ≠ "N",Print["Incorrect choice for monitor value, aborted"];Abort[]];
If[Length[HpDonors]>n && Length[HdDonors]>n,Print["There are more hypothesized contributors than donors. Computation cancelled"]];
If[c<0 || (c>1 && DropinModel=="Linear"), Print["inadmissible choice of dropin parameter. Computation aborted"];Abort[]];
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];MixtureLL=Length/@MixtureFr ;
sortkit=Sort[kit,MixtureLL[[#1]]<MixtureLL[[#2]]&];
res=Table[{Loci[[sortkit[[L]]]],MixtureLL[[sortkit[[L]]]]-1,1,1,1,1,1,"Queued"},{L,Length[kit]}];PrependTo[res,{"Locus","distinct alleles", "Lh Hp","Lh Hd","LR locus","cumulative LR","Log(10,cumu-LR)", "time for locus (s)"}];
AllelesHpDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HpDonors[[i,L,1]],HpDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HpDonors]}];
Do[Do[If[Length[AllelesHpDonors[[i,L]]]==0,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHpDonors[[i,L]]]==1,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HpDonors]}];
AllelesHdDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HdDonors[[i,L,1]],HdDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HdDonors]}];
Do[Do[If[Length[AllelesHdDonors[[i,L]]]==0,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHdDonors[[i,L]]]==1,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HdDonors]}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
Do[start=SessionTime[];res[[L+1,3]]=MixtureProbEqualDropoutTheta[
(* de replica's, alle loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes bekende donoren Hp*)AllelesHpDonors,
(* voor theta *)Complement[AllelesHdDonors,AllelesHpDonors],
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropin parameter *)c,
(* theta *) θ];
res[[L+1,4]]=MixtureProbEqualDropoutTheta[
(* de replica's, alle loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes bekende donoren Hd *)AllelesHdDonors,
(* extra voor theta *)Complement[AllelesHpDonors,AllelesHdDonors],
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropin parameter *)c,
(* theta *) θ];
res[[L+1,5]]=res[[L+1,3]]/res[[L+1,4]];
If[L>1,res[[L+1,6]]=res[[L,6]]*res[[L+1,5]];res[[L+1,7]]=Log[10,res[[L+1,6]]],res[[L+1,6]]=res[[L+1,5]];res[[L+1,7]]=Log[10,res[[L+1,6]]]];
res[[L+1,8]]=SessionTime[]-start,{L,Length[kit]}];
res[[Length[kit]+1,6]]=res[[Length[kit]+1,5]]res[[Length[kit],6]];If[NumberQ[d]&&monitor=="Y",Print[res//MatrixForm]];{Table[res[[k,3]],{k,2,Length[res]}],Table[res[[k,4]],{k,2,Length[res]}],Table[res[[k,5]],{k,2,Length[res]}]}
]

(* the output contains three vectors *)
(* output[[1]] contains the likelihoods under Hp *)
(* output[[2]] contains 1 *)
(* output[[3]] contains 1 *)
(* the order of the loci is obtained by sorting from fewer to more recorded alleles *)
(* see also LR definition *)
HypLikelihoodEqualDropoutTheta[MixturesOnAllLoci_,HpDonors_,HdDonors_,n_,d_,c_,kit_,θ_,monitor_]:=Module[{MixtureAlleles,start=SessionTime[],sortkit},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
If[monitor≠"Y" && monitor ≠ "N",Print["Incorrect choice for monitor value, aborted"];Abort[]];
If[Length[HpDonors]>n && Length[HdDonors]>n,Print["There are more hypothesized contributors than donors. Computation cancelled"]];
If[c<0 || (c>1 && DropinModel=="Linear"), Print["inadmissible choice of dropin parameter. Computation aborted"];Abort[]];
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];MixtureLL=Length/@MixtureFr ;
sortkit=Sort[kit,MixtureLL[[#1]]<MixtureLL[[#2]]&];
res=Table[{Loci[[sortkit[[L]]]],MixtureLL[[sortkit[[L]]]]-1,1,1,1,1,1,"Queued"},{L,Length[kit]}];PrependTo[res,{"Locus","distinct alleles", "Lh Hp","Lh Hd","LR locus","cumulative LR","Log(10,cumu-LR)", "time for locus (s)"}];
AllelesHpDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HpDonors[[i,L,1]],HpDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HpDonors]}];
Do[Do[If[Length[AllelesHpDonors[[i,L]]]==0,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHpDonors[[i,L]]]==1,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HpDonors]}];
AllelesHdDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HdDonors[[i,L,1]],HdDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HdDonors]}];
Do[Do[If[Length[AllelesHdDonors[[i,L]]]==0,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHdDonors[[i,L]]]==1,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HdDonors]}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
Do[start=SessionTime[];res[[L+1,3]]=MixtureProbEqualDropoutTheta[
(* de replica's, alle loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes bekende donoren Hp*)AllelesHpDonors,
(* voor theta *)Complement[AllelesHdDonors,AllelesHpDonors],
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropin parameter *)c,
(* theta *) θ];
res[[L+1,4]]=1;
res[[L+1,5]]=1;
If[L>1,res[[L+1,6]]=res[[L,6]]*res[[L+1,5]];res[[L+1,7]]=Log[10,res[[L+1,6]]],res[[L+1,6]]=res[[L+1,5]];res[[L+1,7]]=Log[10,res[[L+1,6]]]];
res[[L+1,8]]=SessionTime[]-start,{L,Length[kit]}];
res[[Length[kit]+1,6]]=res[[Length[kit]+1,5]]res[[Length[kit],6]];If[NumberQ[d]&&monitor=="Y",Print[res//MatrixForm]];{Table[res[[k,3]],{k,2,Length[res]}],Table[res[[k,4]],{k,2,Length[res]}],Table[res[[k,5]],{k,2,Length[res]}]}
]

(* the output contains three vectors *)
(* output[[1]] contains the likelihoods under Hp *)
(* output[[2]] contains the likelihoods under Hd *)
(* output[[3]] contains the likelihood ratios *)
(* the order of the loci is obtained by sorting from fewer to more recorded alleles *)
(* see also LR definition *)
LRinclLikelihoodsTheta[MixturesOnAllLoci_,HpDonors_,HdDonors_,n_,d_,c_,kit_,θ_,monitor_]:=Module[{MixtureAlleles,start=SessionTime[],sortkit},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
If[monitor≠"Y" && monitor ≠ "N",Print["Incorrect choice for monitor value, aborted"];Abort[]];
If[Length[HpDonors]>n && Length[HdDonors]>n,Print["There are more hypothesized contributors than donors. Computation cancelled"]];
If[c<0 || (c>1 && DropinModel=="Linear"), Print["inadmissible choice of dropin parameter. Computation aborted"];Abort[]];
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];MixtureLL=Length/@MixtureFr ;
sortkit=Sort[kit,MixtureLL[[#1]]<MixtureLL[[#2]]&];
res=Table[{Loci[[sortkit[[L]]]],MixtureLL[[sortkit[[L]]]]-1,1,1,1,1,1,"Queued"},{L,Length[kit]}];PrependTo[res,{"Locus","distinct alleles", "Lh Hp","Lh Hd","LR locus","cumulative LR","Log(10,cumu-LR)", "time for locus (s)"}];
AllelesHpDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HpDonors[[i,L,1]],HpDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HpDonors]}];
Do[Do[If[Length[AllelesHpDonors[[i,L]]]==0,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHpDonors[[i,L]]]==1,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HpDonors]}];
AllelesHdDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HdDonors[[i,L,1]],HdDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HdDonors]}];
Do[Do[If[Length[AllelesHdDonors[[i,L]]]==0,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHdDonors[[i,L]]]==1,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HdDonors]}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
Do[start=SessionTime[];res[[L+1,3]]=MixtureProbTheta[
(* de replica's, alle loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes bekende donoren Hp*)AllelesHpDonors,
(* voor theta *)Complement[AllelesHdDonors,AllelesHpDonors],
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropinkans *)c,
(* theta *) θ];
res[[L+1,4]]=MixtureProbTheta[
(* de replica's, alle loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes bekende donoren Hd *)AllelesHdDonors,
(* extra voor theta *)Complement[AllelesHpDonors,AllelesHdDonors],
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropinkans *)c,
(* theta *) θ];
res[[L+1,5]]=res[[L+1,3]]/res[[L+1,4]];
If[L>1,res[[L+1,6]]=res[[L,6]]*res[[L+1,5]];res[[L+1,7]]=Log[10,res[[L+1,6]]],res[[L+1,6]]=res[[L+1,5]];res[[L+1,7]]=Log[10,res[[L+1,6]]]];
res[[L+1,8]]=SessionTime[]-start,{L,Length[kit]}];
res[[Length[kit]+1,6]]=res[[Length[kit]+1,5]]res[[Length[kit],6]];If[DeleteDuplicates[NumberQ/@d]=={True}&&monitor=="Y",Print[res//MatrixForm]];{Table[res[[k,3]],{k,2,Length[res]}],Table[res[[k,4]],{k,2,Length[res]}],Table[res[[k,5]],{k,2,Length[res]}]}
]

(* the output contains three vectors *)
(* output[[1]] contains the likelihoods under Hp *)
(* output[[2]] contains 1 *)
(* output[[3]] contains the likelihood under Hp again *)
(* the order of the loci is obtained by sorting from fewer to more recorded alleles *)
(* see also LR definition *)
HypLikelihoodTheta[MixturesOnAllLoci_,HpDonors_,HdDonors_,n_,d_,c_,kit_,θ_,monitor_]:=Module[{MixtureAlleles,start=SessionTime[],sortkit},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
If[monitor≠"Y" && monitor ≠ "N",Print["Incorrect choice for monitor value, aborted"];Abort[]];
If[Length[HpDonors]>n && Length[HdDonors]>n,Print["There are more hypothesized contributors than donors. Computation cancelled"]];
If[c<0 || (c>1 && DropinModel=="Linear"), Print["inadmissible choice of dropin parameter. Computation aborted"];Abort[]];
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];
StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureLL=Length/@MixtureFr ;
sortkit=Sort[kit,MixtureLL[[#1]]<MixtureLL[[#2]]&];
res=Table[{Loci[[sortkit[[L]]]],MixtureLL[[sortkit[[L]]]]-1,1,1,1,1,1,"Queued"},{L,Length[kit]}];PrependTo[res,{"Locus","distinct alleles", "Lh Hp","Lh Hd","LR locus","cumulative LR","Log(10,cumu-LR)", "time for locus (s)"}];
AllelesHpDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HpDonors[[i,L,1]],HpDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HpDonors]}];
Do[Do[If[Length[AllelesHpDonors[[i,L]]]==0,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHpDonors[[i,L]]]==1,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HpDonors]}];
AllelesHdDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HdDonors[[i,L,1]],HdDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HdDonors]}];
Do[Do[If[Length[AllelesHdDonors[[i,L]]]==0,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHdDonors[[i,L]]]==1,AppendTo[AllelesHdDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HdDonors]}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
Do[start=SessionTime[];res[[L+1,3]]=MixtureProbTheta[
(* de replica's, alle loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes bekende donoren Hp*)AllelesHpDonors,
(* voor theta *)Complement[AllelesHdDonors,AllelesHpDonors],
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropinkans *)c,
(* theta *) θ];
res[[L+1,4]]=1;
res[[L+1,5]]=res[[L+1,3]]/res[[L+1,4]];
If[L>1,res[[L+1,6]]=res[[L,6]]*res[[L+1,5]];res[[L+1,7]]=Log[10,res[[L+1,6]]],res[[L+1,6]]=res[[L+1,5]];res[[L+1,7]]=Log[10,res[[L+1,6]]]];
res[[L+1,8]]=SessionTime[]-start,{L,Length[kit]}];
res[[Length[kit]+1,6]]=res[[Length[kit]+1,5]]res[[Length[kit],6]];If[DeleteDuplicates[NumberQ/@d]=={True}&&monitor=="Y",Print[res//MatrixForm]];{Table[res[[k,3]],{k,2,Length[res]}],Table[res[[k,4]],{k,2,Length[res]}],Table[res[[k,5]],{k,2,Length[res]}]}
]

Final likelihood ratios definition

LR[MixturesOnAllLoci_,HpDonors_,HdDonors_,n_,d_,c_,kit_,θ_,monitor_]:=
Which[θ>0 && Length[DeleteDuplicates[d]]>1,LRArbitraryDropoutTheta[MixturesOnAllLoci,HpDonors,HdDonors,n,d,c,kit,θ,monitor],
θ>0&&Length[DeleteDuplicates[d]]==1,LREqualDropoutTheta[MixturesOnAllLoci,HpDonors,HdDonors,n,d[[1]],c,kit,θ,monitor],
θ==0&&Length[DeleteDuplicates[d]]>1,LRArbitraryDropoutNoTheta[MixturesOnAllLoci,HpDonors,HdDonors,n,d,c,kit,monitor],
θ==0&&Length[DeleteDuplicates[d]]==1,LREqualDropout[MixturesOnAllLoci,HpDonors,HdDonors,n,d[[1]],c,kit,monitor]]

LRIncludingLikelihoods[MixturesOnAllLoci_,HpDonors_,HdDonors_,n_,d_,c_,kit_,θ_,monitor_]:=
If[θ>0 ,LRinclLikelihoodsTheta[MixturesOnAllLoci,HpDonors,HdDonors,n,d,c,kit,θ,monitor],
LRinclLikelihoodsNoTheta[MixturesOnAllLoci,HpDonors,HdDonors,n,d,c,kit,monitor]]

Likelihood Ratios for Familial Searching

(* probability that a relative with IBD-coefficients listed in IBD of someone with genotype (a1,a2) on locus L has genotype (b1,b2) on this locus *)
GenotypeProbForRelative[IBD_,L_,{b1_,b2_},{a1_,a2_}]:=(2-KD[{b1,b2}])(IBD[[1]]MixtureFr[[L,b1]]MixtureFr[[L,b2]]+IBD[[2]]/4((KD[{a1,b1}]+KD[{a2,b1}])MixtureFr[[L,b2]]+(KD[{a1,b2}]+KD[{a2,b2}])MixtureFr[[L,b1]])+IBD[[3]]/2(KD[{a1,b1}]KD[{a2,b2}]+KD[{a1,b2}]KD[{a2,b1}]))

(* computes the likelihood of the mixtures given in MixturesOnAllLoci, on locus L, assuming that Kgeno are contributors, that there n contributors in total, that their dropout probabilities for heterozygous alleles are given by the vector d, that the dropin parameter is c, and that the first of the unknown donors is related to person K1 according to IBD coefficients given in the IBD vector *)
MixtureProbWithKinship[MixturesOnAllLoci_,L_,Kgeno_,K1_,n_,d_,c_,IBD_]:=Module[{u,Ugeno=Flatten/@Tuples[DeleteDuplicates[Sort/@Tuples[Range[MixtureLL[[L]]],2]],n-Length[Kgeno]],AllGeno={}},
Do[AppendTo[AllGeno,Kgeno[[i,L]]],{i,Length[Kgeno]}];AllGeno=Flatten[AllGeno];
Do[If[Max[Length/@Kgeno[[i]]]>2||Min[Length/@Kgeno[[i]]]<2,Print["impossible genotypes"];Abort[]],{i,Length[Kgeno]}];
AllGeno=Join[AllGeno,#]&/@Ugeno;
AllGeno=Partition[#,2]&/@AllGeno;
u=n-Length[Kgeno];
If[u==0,Print["there has to be at least one unknown. Aborted"]];
If[Length[d]≠n,Print["wrong length dropout vector; computation aborted"];Abort[]];
Sum[Product[SingleMixtureProbDonorsKnown[MixturesOnAllLoci[[t,L]],L,genotypesalledonoren,d,c],{t,Length[MixturesOnAllLoci]}]*10^Total[Log[10,MixtureFr[[L,Take[Flatten[genotypesalledonoren],-2*(u-1)]]]]]*10^Total[Log[10,2-KD/@Take[genotypesalledonoren,-(u-1)]]]*GenotypeProbForRelative[IBD,L,Take[genotypesalledonoren,{-u}]//Flatten,K1[[L]]],{genotypesalledonoren,AllGeno}]]

(* the first of the unknown donors is, under Hp, related via IBD to K1, and under Hd unrelated *)
LRFamilialSearch[MixturesOnAllLoci_,HpDonors_,n_,d_,c_,kit_, K1_, IBD_]:=Module[{MixtureAlleles},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
If[c<0 || (c>1 && DropinModel=="Linear"), Print["inadmissible choice of dropin parameter. Computation aborted"];Abort[]];
If[Length[HpDonors]>n ,Print["There are more hypothesized contributors than donors. Computation cancelled"]];
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];MixtureLL=Length/@MixtureFr ;
AllelesHpDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HpDonors[[i,L,1]],HpDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HpDonors]}];
Do[Do[If[Length[AllelesHpDonors[[i,L]]]==0,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHpDonors[[i,L]]]==1,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HpDonors]}];
AllK1=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{K1[[L,1]],K1[[L,2]]}],{L,NumberOfLoci}]];
Do[If[Length[AllK1[[L]]]==0,AppendTo[AllK1[[L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllK1[[L]]]==1,AppendTo[AllK1[[L]],MixtureLL[[L]]]],{L,NumberOfLoci}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
Table[MixtureProbWithKinship[
(* the replicates, all loci *)AllReplicates,
(* locus *)L,
(* genotypes donors Hp*)AllelesHpDonors,
(* Person K1, whose relative we are testing for *) AllK1,
(* total number of donors *)n,
(* dropout probabilities *)d,
(* dropin parameters *)c,
(* IBD coefficients between first unknown and K1*) IBD]/
MixtureProb[
(* the replicates, all loci *)AllReplicates,
(* locus *)L,
(* genotypes donors Hd *)AllelesHpDonors,
(* total number of donors *)n,
(* dropout probabilities *)d,
(* dropin parameter *)c],{L,kit}]
]

LRFamilialSearchLikelihoods[MixturesOnAllLoci_,HpDonors_,n_,d_,c_,kit_, K1_, IBD_]:=Module[{MixtureAlleles},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
If[c<0 || (c>1 && DropinModel=="Linear"), Print["inadmissible choice of dropin parameter. Computation aborted"];Abort[]];
If[Length[HpDonors]>n ,Print["There are more hypothesized contributors than donors. Computation cancelled"]];
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];MixtureLL=Length/@MixtureFr ;
AllelesHpDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HpDonors[[i,L,1]],HpDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HpDonors]}];
Do[Do[If[Length[AllelesHpDonors[[i,L]]]==0,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHpDonors[[i,L]]]==1,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HpDonors]}];
AllK1=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{K1[[L,1]],K1[[L,2]]}],{L,NumberOfLoci}]];
Do[If[Length[AllK1[[L]]]==0,AppendTo[AllK1[[L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllK1[[L]]]==1,AppendTo[AllK1[[L]],MixtureLL[[L]]]],{L,NumberOfLoci}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
Table[{MixtureProbWithKinship[
(* the replicates, all loci *)AllReplicates,
(* locus *)L,
(* genotypes donors Hp*)AllelesHpDonors,
(* Person K1, whose relative we are testing for *) AllK1,
(* total number of donors *)n,
(* dropout probabilities *)d,
(* dropin parameters *)c,
(* IBD coefficients between first unknown and K1*) IBD],
MixtureProb[
(* the replicates, all loci *)AllReplicates,
(* locus *)L,
(* genotypes donors Hd *)AllelesHpDonors,
(* total number of donors *)n,
(* dropout probabilities *)d,
(* dropin parameter *)c]},{L,kit}]
]

LRFamilialSearchNumerator[MixturesOnAllLoci_,HpDonors_,n_,d_,c_,kit_, K1_, IBD_]:=Module[{MixtureAlleles},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
If[c<0 || (c>1 && DropinModel=="Linear"), Print["inadmissible choice of dropin parameter. Computation aborted"];Abort[]];
If[Length[HpDonors]>n ,Print["There are more hypothesized contributors than donors. Computation cancelled"]];
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];MixtureLL=Length/@MixtureFr ;
AllelesHpDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HpDonors[[i,L,1]],HpDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HpDonors]}];
Do[Do[If[Length[AllelesHpDonors[[i,L]]]==0,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHpDonors[[i,L]]]==1,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HpDonors]}];
AllK1=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{K1[[L,1]],K1[[L,2]]}],{L,NumberOfLoci}]];
Do[If[Length[AllK1[[L]]]==0,AppendTo[AllK1[[L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllK1[[L]]]==1,AppendTo[AllK1[[L]],MixtureLL[[L]]]],{L,NumberOfLoci}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
Table[MixtureProbWithKinship[
(* the replicates, all loci *)AllReplicates,
(* locus *)L,
(* genotypes donors Hp*)AllelesHpDonors,
(* Person K1, whose relative we are testing for *) AllK1,
(* total number of donors *)n,
(* dropout probabilities *)d,
(* dropin parameters *)c,
(* IBD coefficients between first unknown and K1*) IBD],{L,kit}]
]

Likelihood Ratio for traditional (donor) searching, but Hd specifies kinship to the PoI

(* the LR is computed where Hd replaces the last contributor specified by Hpdonors with a person having relatedness defined by IBD w.r.t. this person *)
LRSuspectOpposedToKin[MixturesOnAllLoci_,HpDonors_,n_,d_,c_,kit_, IBD_]:=Module[{MixtureAlleles},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
If[c<0 || (c>1 && DropinModel=="Linear"), Print["inadmissible choice of dropin parameter. Computation aborted"];Abort[]];
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];MixtureLL=Length/@MixtureFr ;
AllelesHpDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HpDonors[[i,L,1]],HpDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HpDonors]}];
Do[Do[If[Length[AllelesHpDonors[[i,L]]]==0,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHpDonors[[i,L]]]==1,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HpDonors]}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
Table[MixtureProb[
(* the replicates, alle loci *)AllReplicates,
(* locus *)L,
(* genotypes donors Hp*)AllelesHpDonors,
(* total number of donors *)n,
(* dropout probs *)d,
(* dropin prob *)c]/MixtureProbWithKinship[
(* the replicates, all loci *)AllReplicates,
(* locus *)L,
(* genotypes donors Hp*)Delete[AllelesHpDonors,-1],
(* Person K1, whose relative we are testing for *) AllelesHpDonors[[-1]],
(* total number of donors *)n,
(* dropout probabilities *)d,
(* dropin parameters *)c,
(* IBD coefficients between first unknown and K1*) IBD],{L,kit}]
]

(* the LR is computed where Hd replaces the last contributor specified by Hpdonors with a person having relatedness defined by IBD w.r.t. this person *)
LRSuspectOpposedToKinLikelihoods[MixturesOnAllLoci_,HpDonors_,n_,d_,c_,kit_, IBD_]:=Module[{MixtureAlleles},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
If[c<0 || (c>1 && DropinModel=="Linear"), Print["inadmissible choice of dropin parameter. Computation aborted"];Abort[]];
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];MixtureLL=Length/@MixtureFr ;
AllelesHpDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HpDonors[[i,L,1]],HpDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HpDonors]}];
Do[Do[If[Length[AllelesHpDonors[[i,L]]]==0,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHpDonors[[i,L]]]==1,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HpDonors]}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
Table[{MixtureProb[
(* the replicates, alle loci *)AllReplicates,
(* locus *)L,
(* genotypes donors Hp*)AllelesHpDonors,
(* total number of donors *)n,
(* dropout probs *)d,
(* dropin prob *)c],
MixtureProbWithKinship[
(* the replicates, all loci *)AllReplicates,
(* locus *)L,
(* genotypes donors Hp*)Delete[AllelesHpDonors,-1],
(* Person K1, whose relative we are testing for *) AllelesHpDonors[[-1]],
(* total number of donors *)n,
(* dropout probabilities *)d,
(* dropin parameters *)c,
(* IBD coefficients between first unknown and K1*) IBD]},{L,kit}]//Transpose
]

(* the LR is computed where Hd replaces the last contributor specified by Hpdonors with a person having relatedness defined by IBD w.r.t. this person *)
LRSuspectOpposedToKinLikelihoodsSortedKit[MixturesOnAllLoci_,HpDonors_,n_,d_,c_,kit_, IBD_]:=Module[{MixtureAlleles,sortkit},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
If[c<0 || (c>1 && DropinModel=="Linear"), Print["inadmissible choice of dropin parameter. Computation aborted"];Abort[]];
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];sortkit=Sort[kit,MixtureLL[[#1]]<MixtureLL[[#2]]&];
Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];MixtureLL=Length/@MixtureFr ;
AllelesHpDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HpDonors[[i,L,1]],HpDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HpDonors]}];
Do[Do[If[Length[AllelesHpDonors[[i,L]]]==0,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHpDonors[[i,L]]]==1,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HpDonors]}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
Table[{MixtureProb[
(* the replicates, alle loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes donors Hp*)AllelesHpDonors,
(* total number of donors *)n,
(* dropout probs *)d,
(* dropin prob *)c],
MixtureProbWithKinship[
(* the replicates, all loci *)AllReplicates,
(* locus *)sortkit[[L]],
(* genotypes donors Hp*)Delete[AllelesHpDonors,-1],
(* Person K1, whose relative we are testing for *) AllelesHpDonors[[-1]],
(* total number of donors *)n,
(* dropout probabilities *)d,
(* dropin parameters *)c,
(* IBD coefficients between first unknown and K1*) IBD]},{L,Length[kit]}]//Transpose
]

Likelihood Ratio distributions for familial searching

(* the person whose relative may be the first unknown in the mixture is replaced by all possibe genotypes, on every locus  *)
(* output is a table with for every locus in kit, the probability distribution of the LRs for randomly drawn PoI under Hp*)
(* this functionality should only be used for familial searches, where the known donors (if any) are the same under both hypotheses *)
MixtureProbsForRandomAllegedRelative[MixturesOnAllLoci_,Donors_,n_,d_,c_,kit_,IBD_]:=Module[{MixtureAlleles},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];MixtureLL=Length/@MixtureFr ;
AllelesDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{Donors[[i,L,1]],Donors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[Donors]}];
Do[Do[If[Length[AllelesDonors[[i,L]]]==0,AppendTo[AllelesDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesDonors[[i,L]]]==1,AppendTo[AllelesDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[Donors]}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
PossibleGenotypes=Table[DeleteDuplicates[Sort/@Tuples[Range[MixtureLL[[L]]],2]],{L,NumberOfLoci}];
Table[Transpose[Table[{MixtureFr[[L,x[[1]]]]MixtureFr[[L,x[[2]]]](2-KroneckerDelta[x[[1]],x[[2]]]),MixtureProbWithKinship[
(* de replica's, alle loci *)AllReplicates,
(* locus *)L,
(* genotypes bekende donoren Hd *)AllelesDonors,
(* Persoon K1, waarvan we de verwant willen weten *) ReplacePart[Partition[RandomProfile[NumberOfLoci],2],L->x],
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropin parameter *)c,
(*verwantschapsvector*) IBD]},{x,PossibleGenotypes[[L]]}]],{L,kit}]]

(* WE GET A SAMPLE OF FAMILIAL SEARCHING LR'S WHERE DONORS_ ARE SUPPOSED TO BE KNOWN AND THE FIRST OF THE UNKNOWN TRACE DONORS IS EVALUATED FOR RELATEDNESS BY IBD WITH A PERSON WHO IS ACTUALLY UNRELATED (FOR EXAMPLE TO BE USED TO SIMULATE A DATABASE SEARCH WHERE EVERYONE IS UNRELATED TO THE TRACE DONORS) *)
LRSampleFromFamilialHd[MixturesOnAllLoci_,Donors_,n_,d_,c_,kit_,IBD_,samplesize_]:=Module[{HdProbs,Hd,simu},
HdProbs=MixtureProbsForRandomAllegedRelative[MixturesOnAllLoci,Donors,n,d,c,kit,IBD];
Hd=Total[Log[10,Table[MixtureProb[
(* de replica's, alle loci *)AllReplicates,
(* locus *)L,
(* genotypes bekende donoren Hd *)AllelesDonors,
(* totaal aantal donoren *)n,
(* dropoutkansen *)d,
(* dropin parameter *)c],{L,kit}]]];
simu=Table[Total[Log[10,Table[RandomChoice[HdProbs[[L,1]]->HdProbs[[L,2]]],{L,Length[kit]}]]],{i,samplesize}]-Hd;simu]

(* WE GET A SAMPLE OF FAMILIAL SEARCHING LR'S WHERE DONORS_ ARE SUPPOSED TO BE KNOWN AND THE FIRST OF THE UNKNOWN TRACE DONORS IS EVALUATED FOR RELATEDNESS BY IBD WITH A PERSON WHO IS ACTUALLY RELATED IN THIS WAY  *)
(* WE DO NOT NEED TO KNOW THE ACTUAL DONORS FOR THIS.
FROM THE MIXTURE WE OBTAIN THE PROBABILITY DISTRIBUTION ON THE DONORS AND HENCE ON THEIR RELATIVES *)
LRSampleFromFamilialHp[MixturesOnAllLoci_,Donors_,n_,d_,c_,kit_,IBD_,samplesize_]:=Module[{HdProbs,Hd,HpProbs,simu},
HdProbs=MixtureProbsForRandomAllegedRelative[MixturesOnAllLoci,Donors,n,d,c,kit,IBD];
Hd=Total[Log[10,Table[MixtureProb[
(* de replica's, alle loci *)AllReplicates,
(* locus *)L,
(* genotypes bekende donoren Hd *)AllelesDonors,
(* totaal aantal donoren *)n,
(* dropoutkansen *)d,
(* dropin parameter *)c],{L,kit}]]];
HpProbs=Table[HdProbs[[i,1]]HdProbs[[i,2]],{i,Length[kit]}];
simu=Table[Total[Log[10,Table[RandomChoice[HpProbs[[L]]->HdProbs[[L,2]]],{L,Length[kit]}]]],{i,samplesize}]-Hd;simu]

(* produces an ROC curve *)
(* assume that Donors (second argument) are  known; they are used as the first donors *)
(* we look for a relative of the first unknwown donor with IBD coefficients specified *)
(* based on a sample from Hp of size samplesize *)
(* plot is based on LR thresholds going up from -5 to the largest integer LR in the sample (rounded below) with stepsize 0.2 *)
FamilialSearchingROC[MixturesOnAllLoci_,Donors_,n_,d_,c_,kit_,IBD_,samplesize_,color_]:=
Module[{sample,max,data},
sample=LRSampleFromFamilialHp[MixturesOnAllLoci,Donors,n,d,c,kit,IBD,samplesize]; 
max=Floor[Max[sample]];
data={};
Do[sample=Select[sample,#>x&];AppendTo[data,{Log[10,Sum[10^-sample[[i]],{i,Length[sample]}]/samplesize],Length[sample]/samplesize}],{x,-5,max,0.2}];
ListPlot[data,PlotRange->All,GridLines->Automatic,Joined->True,PlotMarkers->Automatic,AxesLabel->{"Log(10,FPR)","TPR"},LabelStyle->Larger,PlotStyle->color]
]

Likelihood Ratio distributions for traditional (donor) searching, where PoI may be related to the actual donor

(* probability that a relative with IBD-coefficients listed in IBD of someone with genotype (a1,a2) on locus L has genotype (b1,b2) on this locus *)
GenotypeProbForRelative[IBD_,L_,{b1_,b2_},{a1_,a2_}]:=(2-KD[{b1,b2}])(IBD[[1]]MixtureFr[[L,b1]]MixtureFr[[L,b2]]+IBD[[2]]/4((KD[{a1,b1}]+KD[{a2,b1}])MixtureFr[[L,b2]]+(KD[{a1,b2}]+KD[{a2,b2}])MixtureFr[[L,b1]])+IBD[[3]]/2(KD[{a1,b1}]KD[{a2,b2}]+KD[{a1,b2}]KD[{a2,b1}]))

(* only works for the situation where one PoI is tested, i.e. the p^th of the donors *)
(* this PoI is replaced by all possibe genotypes according to the probabilty that his relative according to IBD has this genotype, on every locus  *)
(* output is a table with for every locus in kit, the probability distribution of the mixture for randomly drawn PoI under Hp*)
MixtureProbsForRandomRelative[MixturesOnAllLoci_,Donors_,p_,n_,d_,c_,kit_,IBD_]:=Module[{MixtureAlleles},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
If[Not[MemberQ[Range[Length[Donors]],p]],Print["wrong index of donor to be replaced"];Abort[]];
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];MixtureLL=Length/@MixtureFr ;
AllelesDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{Donors[[i,L,1]],Donors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[Donors]}];
Do[Do[If[Length[AllelesDonors[[i,L]]]==0,AppendTo[AllelesDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesDonors[[i,L]]]==1,AppendTo[AllelesDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[Donors]}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
PossibleGenotypes=Table[DeleteDuplicates[Sort/@Tuples[Range[MixtureLL[[L]]],2]],{L,NumberOfLoci}];
Table[Transpose[Table[{GenotypeProbForRelative[IBD,L,{x[[1]],x[[2]]},{AllelesDonors[[p,L,1]],AllelesDonors[[p,L,2]]}],MixtureProb[
(* de replica's, alle loci *)AllReplicates,
(* locus *)L,
(* genotypes bekende donoren *)ReplacePart[AllelesDonors,{p,L}->x],
(* totaal aantal donoren *)n,
(* dropoutkans *)d,
(* dropin parameter *)c]},{x,PossibleGenotypes[[L]]}]],{L,kit}]]

(* under Hd the PoI is kept *)
(* so we get the LR in favour of the PoI being a donor versus an IBD-relative being a donor *)
LRSampleFromHdRelatives[MixturesOnAllLoci_,Donors_,p_,n_,d_,c_,kit_,IBD_,samplesize_]:=Module[{HdProbs,Hd,simu},
HdProbs=MixtureProbsForRandomRelative[MixturesOnAllLoci,Donors,p,n,d,c,kit,IBD];
Hd=Total[Log[10,Table[MixtureProb[
(* de replica's, alle loci *)AllReplicates,
(* locus *)L,
(* genotypes bekende donoren Hd *)AllelesDonors,
(* totaal aantal donoren *)n,
(* dropoutkansen *)d,
(* dropin parameter *)c],{L,kit}]]];
simu=Table[Total[Log[10,Table[RandomChoice[HdProbs[[L,1]]->HdProbs[[L,2]]],{L,Length[kit]}]]],{i,samplesize}]-Hd;simu]

(* under Hd the PoI is not kept *)
(* so we get a sample of LR's in favour of being a donor that IBD-relatives of the PoI get, where the PoI is the last one of the specified donors *)
LRSampleFromHdForRelatives[MixturesOnAllLoci_,Donors_,n_,d_,c_,kit_,IBD_,samplesize_]:=Module[{HdProbs,Hd,simu},
HdProbs=MixtureProbsForRandomRelative[MixturesOnAllLoci,Donors,Length[Donors],n,d,c,kit,IBD];
Hd=Total[Log[10,Table[MixtureProb[
(* de replica's, alle loci *)AllReplicates,
(* locus *)L,
(* genotypes bekende donoren Hd *)Delete[AllelesDonors,-1],
(* totaal aantal donoren *)n,
(* dropoutkansen *)d,
(* dropin parameter *)c],{L,kit}]]];
simu=Table[Total[Log[10,Table[RandomChoice[HdProbs[[L,1]]->HdProbs[[L,2]]],{L,Length[kit]}]]],{i,samplesize}]-Hd;simu]

(* non-contributor performance tests *)
(* we replace the last r of the Donors by unknowns under Hd *)
(* the PoI being replaced is the last one of the donors *)
(* this PoI is replaced by IBD-relatives *)
LRNCPTRelatives[MixturesOnAllLoci_,Donors_,r_,n_,d_,c_,kit_,IBD_,samplesize_]:=Module[{HdProbs,Hd,simu},
HdProbs=MixtureProbsForRandomRelative[MixturesOnAllLoci,Donors,Length[Donors],n,d,c,kit,IBD];
Hd=Total[Log[10,Table[MixtureProb[
(* de replica's, alle loci *)AllReplicates,
(* locus *)L,
(* genotypes bekende donoren Hd *)Delete[AllelesDonors,-Partition[Range[r],1]],
(* totaal aantal donoren *)n,
(* dropoutkansen *)d,
(* dropin parameter *)c],{L,kit}]]];
simu=Table[Total[Log[10,Table[RandomChoice[HdProbs[[L,1]]->HdProbs[[L,2]]],{L,Length[kit]}]]],{i,samplesize}]-Hd;simu]

(* LRs we obtain for actual donors, where we look for the first of the unknown donors. Donors are the ones we condition on *)
LRSampleFromHp[MixturesOnAllLoci_,Donors_,n_,d_,c_,kit_,samplesize_]:=Module[{HdProbs,Hd,HpProbs,simu,CalcDonors=Donors},
AppendTo[CalcDonors,Partition[RandomProfile[NumberOfLoci],2]];HdProbs=MixtureProbsForRandomRelative[MixturesOnAllLoci,CalcDonors,Length[CalcDonors],n,d,c,kit,{1,0,0}];
Hd=Total[Log[10,Table[MixtureProb[
(* de replica's, alle loci *)AllReplicates,
(* locus *)L,
(* genotypes bekende donoren Hd *)Delete[AllelesDonors,-1],
(* totaal aantal donoren *)n,
(* dropoutkansen *)d,
(* dropin parameter *)c],{L,kit}]]];
HpProbs=Table[HdProbs[[i,1]]HdProbs[[i,2]],{i,Length[kit]}];
simu=Table[Total[Log[10,Table[RandomChoice[HpProbs[[L]]->HdProbs[[L,2]]],{L,Length[kit]}]]],{i,samplesize}]-Hd;simu]

(* ROC curve if we look for the first unknown person. Donors are the persons we condition on *)
MixtureROC[MixturesOnAllLoci_,Donors_,n_,d_,c_,kit_,samplesize_,color_]:=
Module[{sample,max,data},
sample=LRSampleFromHp[MixturesOnAllLoci,Donors,n,d,c,kit,samplesize]; 
max=Floor[Max[sample]];
data={};
Do[sample=Select[sample,#>x&];AppendTo[data,{Log[10,Sum[10^-sample[[i]],{i,Length[sample]}]/samplesize],Length[sample]/samplesize}],{x,-5,max,0.2}];
ListPlot[data,PlotRange->All,GridLines->Automatic,Joined->True,PlotMarkers->Automatic,AxesLabel->{"Log(10,FPR)","TPR"},LabelStyle->Larger,PlotStyle->color]
]

(* again we look for the first of the unknowns; Donors are the known contributors that we condition on. The first parameters are typically the true ones that we have generated the mixture with and the second set gives the results for donors drawn from this distribution calculated according to the second set of parameters *)
LRSampleFromHpForDifferentParameters[MixturesOnAllLoci_,Donors_,ntrue_,dtrue_,ctrue_,ncalc_,dcalc_,ccalc_,kit_,samplesize_]:=
Module[{CalcDonors=Donors,HdProbsTrue,HdProbsCalc,HdTrue,HdCalc,HpProbsTrue,TrueDonorSample,HpProbsCalc,CalcDonorLRSample,TrueDonorLRSample,maxcalc,maxtrue,datatrue,datacalc},
AppendTo[CalcDonors,Partition[RandomProfile[NumberOfLoci],2]];HdProbsTrue=MixtureProbsForRandomRelative[MixturesOnAllLoci,CalcDonors,Length[CalcDonors],ntrue,dtrue,ctrue,kit,{1,0,0}];
HdProbsCalc=MixtureProbsForRandomRelative[MixturesOnAllLoci,CalcDonors,Length[CalcDonors],ncalc,dcalc,ccalc,kit,{1,0,0}];
HpProbsTrue=Table[HdProbsTrue[[i,1]]HdProbsTrue[[i,2]],{i,Length[kit]}];
HpProbsCalc=Table[HdProbsCalc[[i,1]]HdProbsCalc[[i,2]],{i,Length[kit]}];Do[HpProbsCalc[[i]]=HpProbsCalc[[i]]/Total[HpProbsCalc[[i]]],{i,Length[kit]}];
Do[HpProbsTrue[[i]]=HpProbsTrue[[i]]/Total[HpProbsTrue[[i]]],{i,Length[kit]}];
TrueDonorSample=Table[Table[RandomChoice[HpProbsTrue[[L]]->Range[Length[HpProbsTrue[[L]]]]],{L,Length[kit]}],{samplesize}];
CalcDonorLRSample=Log[10,Table[Product[HpProbsCalc[[L,TrueDonorSample[[i,L]]]],{L,Length[kit]}]/Product[HdProbsCalc[[L,1,TrueDonorSample[[i,L]]]],{L,Length[kit]}],{i,samplesize}]];
TrueDonorLRSample=Log[10,Table[Product[HpProbsTrue[[L,TrueDonorSample[[i,L]]]],{L,Length[kit]}]/Product[HdProbsTrue[[L,1,TrueDonorSample[[i,L]]]],{L,Length[kit]}],{i,samplesize}]];
Transpose[{TrueDonorLRSample,CalcDonorLRSample}]]

(* Donors are the known ones we condition on. We look for the first unknown *)
MixtureROCWithDifferentParameters[MixturesOnAllLoci_,Donors_,ntrue_,dtrue_,ctrue_,ncalc_,dcalc_,ccalc_,kit_,samplesize_,options_]:=
Module[{HdProbsTrue,HdProbsCalc,HdTrue,HdCalc,HpProbsTrue,TrueDonorSample,HpProbsCalc,CalcDonorLRSample,TrueDonorLRSample,maxcalc,maxtrue,datatrue,datacalc},
{TrueDonorLRSample,CalcDonorLRSample}=Transpose[LRSampleFromHpForDifferentParameters[MixturesOnAllLoci,Donors,ntrue,dtrue,ctrue,ncalc,dcalc,ccalc,kit,samplesize]];
Duo=Transpose[{CalcDonorLRSample,TrueDonorLRSample}];
(*Print[Duo[[Range[10]]]];*)
maxcalc=Floor[Max[CalcDonorLRSample]];
maxtrue=Floor[Max[TrueDonorLRSample]];
datatrue={};
Do[TrueDonorLRSample=Select[TrueDonorLRSample,#>x&];AppendTo[datatrue,{Log[10,Sum[10^-TrueDonorLRSample[[i]],{i,Length[TrueDonorLRSample]}]/samplesize],Length[TrueDonorLRSample]/samplesize}],{x,-5,maxtrue,0.2}];
datacalc={};
Do[Duo=Select[Duo,#[[1]]>x&];AppendTo[datacalc,{Log[10,Sum[10^-Duo[[i,2]],{i,Length[Duo]}]/samplesize],Length[Duo]/samplesize}],{x,-5,maxcalc,0.2}];
(*Print[datacalc];*)
ListPlot[{datatrue,datacalc},PlotRange->All,GridLines->Automatic,Joined->True,PlotMarkers->Automatic,AxesLabel->{"Log(10,FPR)","TPR"},LabelStyle->Larger,PlotStyle->options]
]

Genotype probability distributions of the donors

(* We calculate the genotype probabilities for the first of the unknown donors. The specified donors are the ones that we condition on *)
GenotypeDistribution[MixturesOnAllLoci_,Donors_,n_,d_,c_,kit_]:=Module[{HdProbs,Hd,HpProbs,simu,CalcDonors=Donors,PossibleGenotypes,MixtureAllelesRepeatNumber},
MixtureAllelesRepeatNumber=Table[Alleles[[L,#]]&/@MixturesOnAllLoci[[i,L]],{L,NumberOfLoci},{i,Length[MixturesOnAllLoci]}];
MixtureAllelesRepeatNumber=DeleteDuplicates/@Sort/@Flatten/@MixtureAllelesRepeatNumber;
Do[AppendTo[MixtureAllelesRepeatNumber[[L]],"Unseen"],{L,NumberOfLoci}];
MixtureAllelesRepeatNumber=Flatten/@MixtureAllelesRepeatNumber;
AppendTo[CalcDonors,Partition[RandomProfile[NumberOfLoci],2]];
HdProbs=MixtureProbsForRandomRelative[MixturesOnAllLoci,CalcDonors,Length[CalcDonors],n,d,c,kit,{1,0,0}];
PossibleGenotypes=Table[DeleteDuplicates[Sort/@Tuples[MixtureAllelesRepeatNumber[[kit[[L]]]],2]],{L,Length[kit]}];
HpProbs=Table[HdProbs[[i,1]]HdProbs[[i,2]],{i,Length[kit]}];
Do[HpProbs[[L]]=HpProbs[[L]]/Total[HpProbs[[L]]],{L,Length[kit]}];
Do[HpProbs[[L]]=Sort[Transpose[{PossibleGenotypes[[L]],HpProbs[[L]]}],#1[[2]]>#2[[2]]&],{L,Length[kit]}];
Do[HpProbs[[L]]=Select[HpProbs[[L]],#[[2]]≠0&],{L,Length[kit]}];
Do[PrependTo[HpProbs[[L]],{Loci[[kit[[L]]]],{"Prob for donor with dropout" ,d[[Length[Donors]+1]]}}],{L,Length[kit]}];HpProbs]

(* We calculate the genotype probabilities for the first of the unknown donors. The specified donors are the ones that we condition on *)
GenotypeDistributionWithLR[MixturesOnAllLoci_,Donors_,n_,d_,c_,kit_]:=Module[{HdProbs,Hd,HpProbs,simu,CalcDonors=Donors,PossibleGenotypes,MixtureAllelesRepeatNumber},
MixtureAllelesRepeatNumber=Table[Alleles[[L,#]]&/@MixturesOnAllLoci[[i,L]],{L,NumberOfLoci},{i,Length[MixturesOnAllLoci]}];
MixtureAllelesRepeatNumber=DeleteDuplicates/@Sort/@Flatten/@MixtureAllelesRepeatNumber;
Do[AppendTo[MixtureAllelesRepeatNumber[[L]],"Unseen"],{L,NumberOfLoci}];
MixtureAllelesRepeatNumber=Flatten/@MixtureAllelesRepeatNumber;
AppendTo[CalcDonors,Partition[RandomProfile[NumberOfLoci],2]];
HdProbs=MixtureProbsForRandomRelative[MixturesOnAllLoci,CalcDonors,Length[CalcDonors],n,d,c,kit,{1,0,0}];
PossibleGenotypes=Table[DeleteDuplicates[Sort/@Tuples[MixtureAllelesRepeatNumber[[kit[[L]]]],2]],{L,Length[kit]}];
HpProbs=Table[HdProbs[[i,1]]HdProbs[[i,2]],{i,Length[kit]}];
Do[HpProbs[[L]]=HpProbs[[L]]/Total[HpProbs[[L]]],{L,Length[kit]}];
Do[HpProbs[[L]]=Sort[Transpose[{PossibleGenotypes[[L]],HpProbs[[L]],HdProbs[[L,1]],HpProbs[[L]]/HdProbs[[L,1]]}],#1[[2]]>#2[[2]]&],{L,Length[kit]}];
Do[HpProbs[[L]]=Select[HpProbs[[L]],#[[2]]≠0&],{L,Length[kit]}];
Do[PrependTo[HpProbs[[L]],{Loci[[kit[[L]]]],{"Prob for donor with dropout" ,d[[Length[Donors]+1]]}," a priori prob", "LR"}],{L,Length[kit]}];HpProbs]

(* We calculate the genotype probabilities for the first of the unknown donors. The specified donors are the ones that we condition on *)
(* output is in terms of the allelen numbered 1,2,3 etc. An unseen allele is coded 0 *).
GenotypeDistributionNumbered[MixturesOnAllLoci_,Donors_,n_,d_,c_,kit_]:=Module[{HdProbs,Hd,HpProbs,simu,CalcDonors=Donors,PossibleGenotypes,MixtureAlleles},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
If[c<0 || (c>1 && DropinModel=="Linear"), Print["inadmissible choice of dropin parameter. Computation aborted"];Abort[]];
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;
Do[AppendTo[MixtureAlleles[[L]],0],{L,NumberOfLoci}];AppendTo[CalcDonors,Partition[RandomProfile[NumberOfLoci],2]];
HdProbs=MixtureProbsForRandomRelative[MixturesOnAllLoci,CalcDonors,Length[CalcDonors],n,d,c,kit,{1,0,0}];
PossibleGenotypes=Table[DeleteDuplicates[Sort/@Tuples[MixtureAlleles[[kit[[L]]]],2]],{L,Length[kit]}];
HpProbs=Table[HdProbs[[i,1]]HdProbs[[i,2]],{i,Length[kit]}];
Do[HpProbs[[L]]=HpProbs[[L]]/Total[HpProbs[[L]]],{L,Length[kit]}];
Do[HpProbs[[L]]=Sort[Transpose[{PossibleGenotypes[[L]],HpProbs[[L]]}],#1[[2]]>#2[[2]]&],{L,Length[kit]}];
Do[HpProbs[[L]]=Select[HpProbs[[L]],#[[2]]≠0&],{L,Length[kit]}];
HpProbs]

Likelihood ratios for two mixtures to have a common donor

(* We calculate the genotype probabilities for the first of the unknown donors. The specified donors are the ones that we condition on *)
(* the output is a table involving all the genotypes on the locus *)
GenotypeDistributionAllGenotypes[MixturesOnAllLoci_,Donors_,n_,d_,c_,kit_]:=Module[{,HpProbs,AllHpProbs,HdProbs,Hd,simu,CalcDonors=Donors,PossibleGenotypes,MixtureAllelesRepeatNumber,AllGenotypesOnLoci},
MixtureAllelesRepeatNumber=Table[Alleles[[L,#]]&/@MixturesOnAllLoci[[i,L]],{L,NumberOfLoci},{i,Length[MixturesOnAllLoci]}];
MixtureAllelesRepeatNumber=DeleteDuplicates/@Sort/@Flatten/@MixtureAllelesRepeatNumber;
Do[AppendTo[MixtureAllelesRepeatNumber[[L]],"Unseen"],{L,NumberOfLoci}];
MixtureAllelesRepeatNumber=Flatten/@MixtureAllelesRepeatNumber;
AppendTo[CalcDonors,Partition[RandomProfile[NumberOfLoci],2]];
HdProbs=MixtureProbsForRandomRelative[MixturesOnAllLoci,CalcDonors,Length[CalcDonors],n,d,c,kit,{1,0,0}];
PossibleGenotypes=Table[DeleteDuplicates[Sort/@Tuples[MixtureAllelesRepeatNumber[[kit[[L]]]],2]],{L,Length[kit]}];
HpProbs=Table[HdProbs[[i,1]]HdProbs[[i,2]],{i,Length[kit]}];
Do[HpProbs[[L]]=HpProbs[[L]]/Total[HpProbs[[L]]],{L,Length[kit]}];
Do[HpProbs[[L]]=Sort[Transpose[{PossibleGenotypes[[L]],HpProbs[[L]]}]],{L,Length[kit]}];
AllGenotypesOnLoci=Table[DeleteDuplicates[Sort/@Tuples[Alleles[[kit[[L]]]],2]],{L,Length[kit]}];
AllHpProbs=Table[Table[{AllGenotypesOnLoci[[L,t]],},{t,Length[AllGenotypesOnLoci[[L]]]}],{L,Length[kit]}];
Do[Do[AllHpProbs[[L,x,2]]=Which[
MemberQ [MixtureAllelesRepeatNumber[[kit[[L]]]],AllHpProbs[[L,x,1,1]]]&&MemberQ [MixtureAllelesRepeatNumber[[kit[[L]]]],AllHpProbs[[L,x,1,2]]],
Sum[KroneckerDelta[AllHpProbs[[L,x,1]],HpProbs[[L,t,1]]]*HpProbs[[L,t,2]],{t,Length[HpProbs[[L]]]}],
MemberQ [MixtureAllelesRepeatNumber[[kit[[L]]]],AllHpProbs[[L,x,1,1]]]&&Not[MemberQ [MixtureAllelesRepeatNumber[[kit[[L]]]],AllHpProbs[[L,x,1,2]]]],Sum[KroneckerDelta[{AllHpProbs[[L,x,1,1]],"Unseen"},HpProbs[[L,t,1]]]*HpProbs[[L,t,2]]*Fr[[kit[[L]],Flatten[Position[Alleles[[kit[[L]]]],AllHpProbs[[L,x,1,2]]]]]][[1]]/Last[MixtureFr[[kit[[L]]]]],{t,Length[HpProbs[[L]]]}],MemberQ [MixtureAllelesRepeatNumber[[kit[[L]]]],AllHpProbs[[L,x,1,2]]]&&Not[MemberQ [MixtureAllelesRepeatNumber[[kit[[L]]]],AllHpProbs[[L,x,1,1]]]],Sum[KroneckerDelta[{AllHpProbs[[L,x,1,2]],"Unseen"},HpProbs[[L,t,1]]]*HpProbs[[L,t,2]]*Fr[[kit[[L]],Flatten[Position[Alleles[[kit[[L]]]],AllHpProbs[[L,x,1,1]]]]]][[1]]/Last[MixtureFr[[kit[[L]]]]],{t,Length[HpProbs[[L]]]}],Not[MemberQ [MixtureAllelesRepeatNumber[[kit[[L]]]],AllHpProbs[[L,x,1,1]]]]&&Not[MemberQ [MixtureAllelesRepeatNumber[[kit[[L]]]],AllHpProbs[[L,x,1,2]]]],(2-KroneckerDelta[AllHpProbs[[L,x,1,1]],AllHpProbs[[L,x,1,2]]])Fr[[kit[[L]],Flatten[Position[Alleles[[kit[[L]]]],AllHpProbs[[L,x,1,1]]]]]][[1]]Fr[[kit[[L]],Flatten[Position[Alleles[[kit[[L]]]],AllHpProbs[[L,x,1,2]]]]]][[1]]*Last[HpProbs[[L]]][[2]]/(Last[MixtureFr[[kit[[L]]]]])^2],{x,Length[AllHpProbs[[L]]]}],{L,Length[kit]}];AllHpProbs]

(* LR for the first of the unknown persons in each mixture to be the same versus unrelated *)
(* each mixture with its own number of persons, dropout probabilities, number of replicates *) 
LRSameDonor[Mixture1OnAllLoci_,Mixture2OnAllLoci_,Donors1_,Donors2_,n1_,d1_,c1_,n2_,d2_,c2_,kit_]:=Module[{genomix1=GenotypeDistributionAllGenotypes[Mixture1OnAllLoci,Donors1,n1,d1,c1,kit],genomix2=GenotypeDistributionAllGenotypes[Mixture2OnAllLoci,Donors2,n2,d2,c2,kit],result},result=Table[Sum[If[genomix1[[L,x,2]]>0&& genomix2[[L,x,2]]>0,genomix1[[L,x,2]]*genomix2[[L,x,2]]/((2-KroneckerDelta[genomix1[[L,x,1,1]],genomix1[[L,x,1,2]]])*Fr[[kit[[L]],Flatten[Position[Alleles[[kit[[L]]]],genomix1[[L,x,1,1]]]]]][[1]]*Fr[[kit[[L]],Flatten[Position[Alleles[[kit[[L]]]],genomix1[[L,x,1,2]]]]]][[1]]),0],{x,Length[genomix1[[L]]]}],{L,Length[kit]}];result]

Likelihood ratios for an unknown donor to be sibling of two typed siblings

T3[a_,b_,c_]:=KroneckerDelta[a,b]KroneckerDelta[a,c]

SibProb[L_,a_,b_,c_,d_]:=1/8*(KD[{a,c}]MixtureFr[[L,a]]KD[{b,d}]MixtureFr[[L,b]]+KD[{a,c}]MixtureFr[[L,a]]MixtureFr[[L,b]]MixtureFr[[L,d]]+KD[{a,d}]MixtureFr[[L,a]]KD[{b,c}]MixtureFr[[L,b]]+KD[{a,d}]MixtureFr[[L,a]]MixtureFr[[L,b]]MixtureFr[[L,c]]+KD[{b,c}]MixtureFr[[L,b]]MixtureFr[[L,a]]MixtureFr[[L,d]]+KD[{b,d}]MixtureFr[[L,b]]MixtureFr[[L,a]]MixtureFr[[L,c]]+2MixtureFr[[L,a]]MixtureFr[[L,b]]MixtureFr[[L,c]]MixtureFr[[L,d]])

P1[L_,a_,b_,c_,d_,e_,f_]:=1/(64SibProb[L,a,b,c,d])*(T3[a,c,e]MixtureFr[[L,a]]T3[b,d,f]MixtureFr[[L,b]]+T3[a,c,e]MixtureFr[[L,a]]KD[{b,d}]MixtureFr[[L,b]]MixtureFr[[L,f]]+T3[a,c,f]MixtureFr[[L,a]]T3[b,d,e]MixtureFr[[L,b]]+T3[a,c,f]MixtureFr[[L,a]]KD[{b,d}]MixtureFr[[L,b]]MixtureFr[[L,e]]+T3[b,d,f]MixtureFr[[L,b]]KD[{a,c}]MixtureFr[[L,a]]MixtureFr[[L,e]]+T3[b,d,e]MixtureFr[[L,b]]KD[{a,c}]MixtureFr[[L,a]]MixtureFr[[L,f]]+2MixtureFr[[L,e]]MixtureFr[[L,f]]KD[{a,c}]MixtureFr[[L,a]]KD[{b,d}]MixtureFr[[L,b]])

P2[L_,a_,b_,c_,d_,e_,f_]:=1/(64SibProb[L,a,b,c,d])*(   T3[a,c,e]MixtureFr[[L,a]]KD[{b,f}]MixtureFr[[L,b]]MixtureFr[[L,d]]+T3[a,c,e]MixtureFr[[L,a]]KD[{d,f}]MixtureFr[[L,d]]MixtureFr[[L,b]]+T3[a,c,f]MixtureFr[[L,a]]KD[{b,e}]MixtureFr[[L,b]]MixtureFr[[L,d]]+T3[a,c,f]MixtureFr[[L,a]]KD[{d,e}]MixtureFr[[L,d]]MixtureFr[[L,b]]+KD[{a,c}]MixtureFr[[L,a]]KD[{b,f}]MixtureFr[[L,b]]MixtureFr[[L,d]]MixtureFr[[L,e]]+KD[{a,c}]MixtureFr[[L,a]]KD[{d,f}]MixtureFr[[L,d]]MixtureFr[[L,b]]MixtureFr[[L,e]]+KD[{a,c}]MixtureFr[[L,a]]KD[{b,e}]MixtureFr[[L,b]]MixtureFr[[L,d]]MixtureFr[[L,f]]+KD[{a,c}]MixtureFr[[L,a]]KD[{d,e}]MixtureFr[[L,d]]MixtureFr[[L,b]]MixtureFr[[L,f]])

P7[L_,a_,b_,c_,d_,e_,f_]:=1/(64SibProb[L,a,b,c,d])*(KD[{a,e}]MixtureFr[[L,a]]KD[{b,f}]MixtureFr[[L,b]]MixtureFr[[L,c]]MixtureFr[[L,d]]+KD[{a,e}]MixtureFr[[L,a]]KD[{d,f}]MixtureFr[[L,d]]MixtureFr[[L,c]]MixtureFr[[L,b]]+KD[{a,f}]MixtureFr[[L,a]]KD[{b,e}]MixtureFr[[L,b]]MixtureFr[[L,c]]MixtureFr[[L,d]]+KD[{a,f}]MixtureFr[[L,a]]KD[{d,e}]MixtureFr[[L,d]]MixtureFr[[L,b]]MixtureFr[[L,c]]+KD[{c,e}]MixtureFr[[L,c]]KD[{b,f}]MixtureFr[[L,b]]MixtureFr[[L,a]]MixtureFr[[L,d]]+KD[{c,e}]MixtureFr[[L,c]]KD[{d,f}]MixtureFr[[L,d]]MixtureFr[[L,a]]MixtureFr[[L,b]]+KD[{c,f}]MixtureFr[[L,c]]KD[{b,e}]MixtureFr[[L,b]]MixtureFr[[L,a]]MixtureFr[[L,d]]+KD[{c,f}]MixtureFr[[L,c]]KD[{d,e}]MixtureFr[[L,d]]MixtureFr[[L,a]]MixtureFr[[L,b]])

(* probability that a sibling of persons with genotypes (a1,a2) and (b1,b2) on locus L has genotype (c1,c2) on this locus *)
GenotypeProbForThirdSib[L_,{e_,f_},{a_,b_},{c_,d_}]:=(2-KD[{e,f}])(P1[L,a,b,c,d,e,f]+P1[L,a,b,d,c,e,f]+P2[L,a,b,c,d,e,f]+P2[L,a,b,d,c,e,f]+P2[L,b,a,d,c,e,f]+P2[L,b,a,c,d,e,f]+P7[L,a,b,c,d,e,f]+P7[L,a,b,d,c,e,f])

(* computes the likelihood of the mixtures given in MixturesOnAllLoci, on locus L, assuming that Kgeno are contributors, that there n contributors in total, that their dropout probabilities for heterozygous alleles are given by the vector d, that the dropin parameter is c, and that the first of the unknown donors is a sibling of known sibs K1 and K2  *)
MixtureProbThirdSib[MixturesOnAllLoci_,L_,Kgeno_,K1_,K2_,n_,d_,c_]:=Module[{u,Ugeno=Flatten/@Tuples[DeleteDuplicates[Sort/@Tuples[Range[MixtureLL[[L]]],2]],n-Length[Kgeno]],AllGeno={}},
Do[AppendTo[AllGeno,Kgeno[[i,L]]],{i,Length[Kgeno]}];AllGeno=Flatten[AllGeno];
Do[If[Max[Length/@Kgeno[[i]]]>2||Min[Length/@Kgeno[[i]]]<2,Print["impossible genotypes"];Abort[]],{i,Length[Kgeno]}];
AllGeno=Join[AllGeno,#]&/@Ugeno;
AllGeno=Partition[#,2]&/@AllGeno;
u=n-Length[Kgeno];
If[u==0,Print["there has to be at least one unknown. Aborted"]];
If[Length[d]≠n,Print["wrong length dropout vector; computation aborted"];Abort[]];
Sum[Product[SingleMixtureProbDonorsKnown[MixturesOnAllLoci[[t,L]],L,genotypesalledonoren,d,c],{t,Length[MixturesOnAllLoci]}]*10^Total[Log[10,MixtureFr[[L,Take[Flatten[genotypesalledonoren],-2*(u-1)]]]]]*10^Total[Log[10,2-KD/@Take[genotypesalledonoren,-(u-1)]]]*GenotypeProbForThirdSib[L,Take[genotypesalledonoren,{-u}]//Flatten,K1[[L]],K2[[L]]],{genotypesalledonoren,AllGeno}]]

(* H1: the first of the unknown donors is a sibling of both K1 and K2 *)
(* H2: unrelated *) 
LRFirstUnknownThirdSib[MixturesOnAllLoci_,HpDonors_,n_,d_,c_,kit_, K1_, K2_]:=Module[{MixtureAlleles},
MixtureAlleles=Table[{},{L,NumberOfLoci}];
(*If[monitor≠"Y" && monitor ≠ "N",Print["Incorrect choice for monitor value, aborted"];Abort[]];*)
If[c<0 || (c>1 && DropinModel=="Linear"), Print["inadmissible choice of dropin parameter. Computation aborted"];Abort[]];
If[Length[HpDonors]>n >n,Print["There are more hypothesized contributors than donors. Computation cancelled"]];
Do[Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],MixturesOnAllLoci[[r,L]]]],{L,NumberOfLoci}],{r,Length[MixturesOnAllLoci]}];
Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],K1[[L]]]],{L,NumberOfLoci}];
Do[MixtureAlleles[[L]]=Flatten[AppendTo[MixtureAlleles[[L]],K2[[L]]]],{L,NumberOfLoci}];
MixtureAlleles=DeleteDuplicates/@Sort/@MixtureAlleles;StutterPositions=Table[ParentOfMixAllele[L,#,MixtureAlleles]&/@Range[Length[MixtureAlleles[[L]]]],{L,NumberOfLoci}];
Do[AppendTo[StutterPositions[[L]],0],{L,NumberOfLoci}];
MixtureFr=Table[Table[Fr[[L,MixtureAlleles[[L,i]]]],{i,Length[MixtureAlleles[[L]]]}],{L,NumberOfLoci}];Do[AppendTo[MixtureFr[[L]],1-Total[MixtureFr[[L]]]],{L,NumberOfLoci}];MixtureLL=Length/@MixtureFr ;
AllelesHpDonors=Table[Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{HpDonors[[i,L,1]],HpDonors[[i,L,2]]}],{L,NumberOfLoci}]],{i,Length[HpDonors]}];
Do[Do[If[Length[AllelesHpDonors[[i,L]]]==0,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllelesHpDonors[[i,L]]]==1,AppendTo[AllelesHpDonors[[i,L]],MixtureLL[[L]]]],{L,NumberOfLoci}],{i,Length[HpDonors]}];
AllK1=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{K1[[L,1]],K1[[L,2]]}],{L,NumberOfLoci}]];
Do[If[Length[AllK1[[L]]]==0,AppendTo[AllK1[[L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllK1[[L]]]==1,AppendTo[AllK1[[L]],MixtureLL[[L]]]],{L,NumberOfLoci}];
AllK2=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@{K2[[L,1]],K2[[L,2]]}],{L,NumberOfLoci}]];
Do[If[Length[AllK2[[L]]]==0,AppendTo[AllK2[[L]],MixtureLL[[L]]]],{L,NumberOfLoci}];Do[If[Length[AllK2[[L]]]==1,AppendTo[AllK2[[L]],MixtureLL[[L]]]],{L,NumberOfLoci}];
AllReplicates=Table[{},{r,Length[MixturesOnAllLoci]}];
Do[AllReplicates[[r]]=Quiet[Table[Flatten[Position[MixtureAlleles[[L]],#]&/@MixturesOnAllLoci[[r,L]]],{L,NumberOfLoci}]],{r,Length[MixturesOnAllLoci]}];
Table[MixtureProbThirdSib[
(* the replicates, all loci *)AllReplicates,
(* locus *)L,
(* genotypes donors Hp*)AllelesHpDonors,
(* Person K1, first of the known siblings *) AllK1,
(* Person K2, second known sibling *) AllK2,
(* total number of donors *)n,
(* dropout probabilities *)d,
(* dropin parameters *)c]/
MixtureProb[
(* the replicates, all loci *)AllReplicates,
(* locus *)L,
(* genotypes donors Hd *)AllelesHpDonors,
(* total number of donors *)n,
(* dropout probabilities *)d,
(* dropin parameter *)c],{L,kit}]
]