Simulering av sammanträffanden - I

Författare: Håkan Kjellerstrand (hakank@bonetmail.com)
Dokumentet påbörjades: 2003-10-14

Eventuellt finns kommentarer på detta dokument i blogg-anteckningen Simulering av sammanträffanden.

Inledning

När jag läste Persi Diaconis och Frederick Mostellers underbara paper "Methods for Studying Coincidences" skrev jag lite R-kod för att exemplifiera vissa avsnitt. Några exempel publicerades i Sammanträffanden - anteckningar vid läsning av Diaconis och Mosteller 'Methods for Studying Coincidences'.

Dock skrev jag då ingen simulering av hur sammanträffanden slumpmässigt kan uppstå, vilket nu åtgärdas. Så här är lite R-kod för att simulera sammanträffanden. Mer om R finns på www.r-project.org. R-koden till programmet finns nedan.

Maila mig gärna om du har frågor eller kommentarer.

Simulering - vad gör programmet

Programmet simulerar en grupp av n agenter (personer etc) som har ett antal attribut med ett visst antal olika värden, t.ex. födelsedagar (bland 365 dagar), gått i en viss skola (bland, säg, 500 möjliga skolor) etc. Detta görs genom att för varje agent och respektive attribut slumpas ett heltalsvärde som utgör agentens värde för detta attributet.

När två agenter har samma attribut + värde (samma attribut-värde-par, AVP) räknas det som ett sammanträffande. Exempel: Om agenterna 12 och 17 båda råkar ha värde 24 för attributet 2 innebär det ett sammanträffande mellan dessa agenter.

Det kan även förekomma att dessa agenter har flera sådana AVP gemensamma. Då räknas det fortfarande som ett sammanträffande. Antalet matchningar av AVP för agenterna visas speciellt i en lista (ss.list).

(Vi räknar ett sammanträffande mellan två agenter endast en gång. Dvs "agent X har samma födelsedag som agent Y" är samma som "agent Y har samma födelsedag som agent X".)

Ett litet exempel
Låt oss anta att vi vill studera sannolikheten för sammanträffande i en grupp med 100 personer, där attributen är följande: Vad är då sannolikheten att det finns några personer som har något av följande sammanträffanden: Naturligtvis kan man kombinera dessa sammanträffanden, t.ex. har samma födelsedag och gått på samma skola.

Det berömda födelsedagsproblemet utgörs av det enkla sammanträffande-attributet samma födelsedag. Se nedan under Födelsedagsproblemet.

Enkla attribut
Notera att vi här endast hanterar enkla attribut hos en agent. Mer komplicerade relationer mellan agenterna X och Y och andra agenter såsom X känner Z som känner Y hanteras alltså inte. Vi har endast "samma AVP" som kriterium för ett sammanträffande.

Vill man kan man översätta de mer komplexa relationer till enkla attribut, men min känsla är att komplexa (transitiva) relationer skapar ännu fler sammanträffanden än rena AVP-sammanträffanden. Detta är något man kan fundera vidare på. Se även slutkommentaren om nätverksanalys nedan.

Här är en minimal skiss på hur man kan studera transitiva relationer (kursivt avsnitt)

Programmets parametrar

Programmet tar följande parametrar. (Om värdet inte finns med i anropet till programmet används de default-parametrar som anges i funktionsdefinitionen. Man behöver alltså inte ange t.ex. min.coincidences=1 i anropet.)

En exempelkörning

Här simulerar vi med 100 agenter, tre attribut med respektive 365, 500, och 1000 möjliga värden. Låt anta att det är respektive födelsedagar på ett år, antal teaterföreställningar samt antal skolor i landet. Frågan är hur sannolikt det är att det uppstår (slumpmässigt) åtminstone ett sammanträffande bland denna grupp av agenter. range är 0, och vi studerar de sammanträffanden som är lika med eller större än 1 (dvs min.coincidences=1). 
> coincidence.test(num.agents=100, num.att=3, num.vals=c(365,500,1000), min.coincidences=1,range=0)
$ss.list
 [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

$counts
[1] 28

$coincidence.limit
[1] 1

$p.50
[1] 15.83929

$p.95
[1] 32.99852

$range
[1] 0

$num.connections
[1] 10000
Resultat ska tolkas enligt följande: De teoretiska värdena p.50 och p.95 är uträknade enligt den approximation som finns i Diaconis och Mostellers "Methods for Studying Coincidences" och som även finns i min anteckning Sammanträffanden - anteckningar vid läsning av Diaconis och Mosteller 'Methods for Studying Coincidences'. Notera att de gäller för range=0.

Så, ovanstående exempel visar att av de 100 agenterna som fannns i gruppen var det (i denna körning) 28 par som hade ett sammanträffande. Via listan ss.list ser vi att samtliga dessa par bara hade ett attribut-värde-par som var samma.

Naturligvis utgör ingen enstaka simulering någon riktigt tillförlitlig indikation på sannolikheterna. För att få lite mer statistisk stuns på det hela måste vi simulera lite fler fall. Här görs en simulering av 100 körningar, där vi endast bryr oss om hur många sammanträffanden som fanns, dvs resultatet för count (vi bryr oss alltså inte om antalet AVP per sammanträffande). 
> xx <- replicate(100, coincidence.test(num.agents=100, num.att=3, num.vals=c(365,500,1000), min.coincidences=1)$counts)
> xx
  [1] 27 31 27 31 28 19 35 29 29 44 24 35 24 33 45 23 33 21 31 39 30 30 32 30 32
 [26] 24 28 32 28 35 26 24 36 26 30 21 29 26 15 29 24 34 31 27 20 26 20 26 19 30
 [51] 21 24 39 28 22 30 24 33 25 27 24 25 27 29 34 23 27 30 26 33 31 24 32 37 31
 [76] 25 42 33 27 32 29 20 27 39 28 25 25 33 34 29 37 28 30 23 30 25 31 32 43 26
Vi gör nu en tabulering som visar att det är väldigt många olika värden som kan uppkomma. Dvs det finns 1 körning som hade 15 sammanträffanden, 6 körningar som hade 25 sammanträffande etc. 
> table(xx)
xx
15 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 39 42 43 44 45 
 1  2  3  3  1  3  9  6  7  8  6  7  9  7  6  6  3  3  1  2  3  1  1  1  1
En lite mer överskådlig uppställning får vi via: 
> describe(xx)
xx 
      n missing  unique    Mean     .05     .10     .25     .50     .75     .90 
    100       0      25   28.87    20.0    22.9    25.0    29.0    32.0    35.1 
    .95 
   39.0 

lowest : 15 19 20 21 22, highest: 39 42 43 44 45
Dvs, medelvärdet är cirka 29 sammanträffanden. I 5% av fallen får vi under 20 sammanträffande respektive 5% över 39 sammanträffanden. Detta innebär att om det i en riktig undersökning av en (slumpvis utvald) grupp skulle finnas fler än cirka 39 sammanträffanden är det troligen ingen slump utan det finns någon annan förklaring bakom (tvillingar, vänner etc). Om det är färre än 20 sammanträffanden är det också något som bör kontrolleras vidare. Men värden mellan 20 och 39 kan man förvänta sig utan att behöva höja på ögronbrynen.

Låt oss nu titta lite på de teoretiska värdena (p.50 respektive p.95). Vi ser i resultatet ovan att p.50 är cirka 16 agenter och p.95 är cirka 33, dvs för att det ska vara 50% (respektive 95%) chans att det ska bli åtminstone ett sammanträffande krävs det cirka 16 (respektive 33) agenter. Vi börjar med 50%-gränsen. Låt oss se vad som händer om vi endast har 16 agenter i gruppen. Så, av de 100 körningarna fördelar sig antalet sammanträffanden enligt följande. Här läggs även in den procentuella andelen (table(xx)/100)) för respektive antal sammanträffanden. 
> xx <- replicate(100,coincidence.test(num.agents=16, num.att=3, num.vals=c(365,500,1000), min.coincidences=1)$counts)
> xx
  [1] 0 0 1 0 1 0 1 1 0 2 0 0 0 0 1 0 1 1 1 0 1 1 3 0 0 1 0 1 0 0 1 0 1 1 1 1 1
 [38] 0 0 1 0 1 1 1 2 1 0 0 0 1 1 1 0 0 1 0 0 0 0 1 1 0 0 0 1 2 1 2 0 3 0 1 1 0
 [75] 1 0 1 2 0 0 1 0 0 2 2 1 1 0 1 1 0 0 0 2 2 2 1 0 1 1

> table(xx)
xx
 0  1  2  3 
45 43 10  2 

> describe(xx)
xx 
      n missing  unique    Mean 
    100       0       4    0.69 

0 (45, 45%), 1 (43, 43%), 2 (10, 10%), 3 (2, 2%)
Här ser vi att det finns 45 fall där det inte fanns något sammanträffande alls, vilket ju stämmer rätt bra med teorin. När jag gjorde en simulering med 15 agenter var det 51% körningar som inte hade något sammanträffande.

Om vi nu kollar 95%-gränsen (dvs 33 agenter) får vi följande: 
> xxx <- replicate(100,coincidence.test(num.agents=33, num.att=3, num.vals=c(365,500,1000), min.coincidences=1)$counts)

> table(xxx)
xxx
 0  1  2  3  4  5  6  7  8  9 
 4 12 20 23 18 13  6  2  1  1 

> table(xxx)/100
xxx
   0    1    2    3    4    5    6    7    8    9 
0.04 0.12 0.20 0.23 0.18 0.13 0.06 0.02 0.01 0.01 


> describe(xxx)
xxx 
      n missing  unique    Mean     .05     .10     .25     .50     .75     .90 
    100       0      10    3.25     1.0     1.0     2.0     3.0     4.0     5.1 
    .95 
    6.0 

          0  1  2  3  4  5 6 7 8 9
Frequency 4 12 20 23 18 13 6 2 1 1
%         4 12 20 23 18 13 6 2 1 1
Vi kan notera att det är 4 av 100 körningar (dvs 4%) som inte har något sammanträffande alls, och det stämmer bra med det teoretiska värdet 5% (dvs 100% - 95%). Medelvärdet är drygt 3 sammanträffanden. Vi kan här även notera att det rent slumpmässigt kan uppstå ända upp till 9 sammanträffanden.

Simulerar man med fler körningar (t.ex. 1000) får man ett mer stabilt resultat. Det tar dock en hel del mer tid. Vilket härmed görs: 
# 15 agenter
> xx<-replicate(1000, coincidence.test(num.agents=15, num.att=3, num.vals=c(365,500,1000), min.coincidences=1)$counts)

> table(xx)/1000
xx
    0     1     2     3     4     6 
0.549 0.326 0.102 0.020 0.002 0.001 

> describe(xx)
xx 
      n missing  unique    Mean 
   1000       0       6   0.604 

            0   1   2  3 4 6
Frequency 549 326 102 20 2 1
%          55  33  10  2 0 0

# 33 agenter
> xxx<-replicate(1000, coincidence.test(num.agents=33, num.att=3, num.vals=c(365,500,1000), min.coincidences=1)$counts)

# fördelningen i procent
> table(xxx)/1000
xxx
    0     1     2     3     4     5     6     7     8     9    10 
0.045 0.149 0.235 0.236 0.160 0.089 0.053 0.017 0.010 0.003 0.003 

> describe(xxx)
xxx 
      n missing  unique    Mean     .05     .10     .25     .50     .75     .90 
   1000       0      11   2.986       1       1       2       3       4       5 
    .95 
      6 

           0   1   2   3   4  5  6  7  8  9 10
Frequency 45 149 235 236 160 89 53 17 10  3  3
%          4  15  24  24  16  9  5  2  1  0  0
Om man järmför 100-simuleringen med 1000-simuleringen ser vi inte så stora skillnader. Exakt antal simuleringar man måste köra beror på problemets art, precisionen man vill ernå samt den tid man tycker man har att vänta på resultatet.

Födelsedagsproblemet

Här är en simulering av födelsedagsproblemet, som också nämns i blogganteckningen Sammanträffanden - anteckningar vid läsning av Diaconis och Mosteller 'Methods for Studying Coincidences'.

Sannolikheten för att det åtminstone finns två personer i en grupp av k personer som delar samma attributvärde av N möjliga är (pseudo-kod): 
   1 - Prod( 1 - i/ N, i = 1..k-1)
I R skrivs födelsedagsproblemet (N=365, k=23 resp k=22) på detta sätt 
> N <- 365
> k <- 23
> x <- 1:(k-1)
> 1 - prod(1 - x / N)
[1] 0.5072972

> k <- 22
> x <- 1:(k-1)
> 1 - prod(1 - x / N)
[1] 0.4756953
OK, låt oss nu simulera födelsedagsproblemet med coincidence.test() för respektive 22 och 23 personer. 
> xx <- replicate(100, coincidence.test(num.agents=22, num.att=1, num.vals=365)$count)
table(xx)
> table(xx)
xx
 0  1  2  3 
53 33 13  1

> xx <- replicate(100, coincidence.test(num.agents=23, num.att=1, num.vals=365)$count);
> table(xx)
xx
 0  1  2  3  4 
48 32 17  2  1
Vilket stöder de teoretiska värdena rätt bra.

Multiple endpoints

Parametern range styr nästan-matchningar, dvs att attributvärdena inte behöver vara exakta utan kan utgöra ett intervall. Se även här Sammanträffanden ... .

Som exempel på effekten av detta kör vi födelsedagsproblemet med lite olika värden på range, och för konstant antal agenter (22). Det räcker att se hur värdet för 0 sammanträffanden ändras. 
# range = 1
>  xx <- replicate(100, coincidence.test(num.agents=22, num.att=1, num.vals=365,range=1)$count)
table(xx)
> table(xx)
xx
 0  1  2  3  4  5  6  8 
15 30 24 14  8  7  1  1 

# range = 2
>  xx <- replicate(100, coincidence.test(num.agents=22, num.att=1, num.vals=365,range=2)$count)
table(xx)
> table(xx)
xx
 0  1  2  3  4  5  6  7  8 10 
 2 12 18 28 18 11  4  3  3  1 


# range = 10 (dvs, 10 dagar före eller 10 dagar efter)
> table(xx)
xx
 7  8  9 10 11 12 13 14 15 16 17 18 19 22 
 1  4  6 10 13 11 13 13 13  3  7  2  3  1
Låt oss nu avsluta med endast 6 agenter och range=10 : 
>  xx <- replicate(100, coincidence.test(num.agents=6, num.att=1, num.vals=365,range=10)$count)
> table(xx)
xx
 0  1  2  3 
38 43 13  6
Det approximativa teoretiska antalet agenter som krävs för en fifty-fifty-chans med range 10 är cirka 5 (se formel i den refererade blogganteckningen): 
> range <- 10
> c     <- 365
> 1.2 * sqrt(c/(2*range + 1))
[1] 5.002856

Matematisk analys och resampling

Denna typ av problem går säkerligen att lösa med hjälp av rena matematiska metoder (vissa av dessa har vi sett). Genom att skriva simuleringsprogram (av vissa kallat resampling) är det enkelt att lägga till funktioner eller införa små förändringar av problemet, även där matematiska metoder blir mycket svåra (måhända undantaget professorer i matematisk statistik :-).

Även om man löser med matematiska metoder är simuleringar ofta ett bra sätt att få en känsla för problemet, kontrollera att ens uträkningar är korrekta och ibland helt enkelt lättare att göra. Det kan också vara ett bra sätt att generera data för vidare analyser.

Lite mer om resampling och simulering finns på Lite om resampling, simulering, sannolikhetsproblem etc..

Nätverksanalys

Det är naturligtvis ingen tillfällighet att problemet med sammanträffanden påminner om Small World-problemet och social nätverksanalys, dvs forskningen kring olika typer av nätverk, såväl sociala som andra, t.ex. webben, författare som skrivit artiklar tillsammans etc. Se vidare t.ex. Social Network Analysis och Complex networks.

Troligen kommer jag kika lite mer på detta vid tillfälle. Förhoppnings kan jag då använda min minimala skiss.

R-funktionen coincidence.test()

Här är R-funktionen coincidence.test(). Parametrarna för programet förklaras ovan. 
coincidence.test <- function(num.agents=10, num.att=4, num.vals=4, range=NULL, min.coincidences=1) {
  # some flexibility for num.vals
  if (length(num.vals) == 1) {
     num.vals <- rep(num.vals, num.att)
  }

  # theoretical values (number of persons required for 50% and 95% chance)
  # according to Diaconics & Mosteller "Methods for Studying Coincidences"
  p.50 <- 1.2*sqrt(1/sum(1/num.vals))
  p.95 <- 2.5*sqrt(1/sum(1/num.vals))
  num.connections <- num.agents*num.agents

  # build matrix
  mat <- matrix(0, nrow=num.agents, ncol=num.att)
  for (i in 1:num.att) {
    cc <- sample(1:num.vals[i], num.agents, replace=T)
    mat[,i] <- cc
  }

  ss <- NULL
  # just check the lower triangle
  if (identical(range, NULL)) {
    ss <- sapply(1:num.agents, function(x) sapply(1:num.agents, function(y)  if (y > x) {sum(mat[x,]==mat[y,])} else {0}))
  } else {
    ss <- sapply(1:num.agents, function(x) 
              sapply(1:num.agents, function(y)  
                  if (y > x) {
                     sum( (abs(mat[x,] - mat[y]) <= range) )
                  } else {
                     0
                  }
              ) # end sapply y
          ) # end sapply x
  }

  # just those with proper number of coincidences
  ss.list <- ss[ss>=min.coincidences] 

  counts <- length(ss.list)

  list(ss.list=ss.list,counts=counts,coincidence.limit=min.coincidences,p.50=p.50,p.95=p.95, range=range, num.connections=num.connections)

} # end coincidence.test
Se också: Created by Hakan Kjellerstrand (hakank@bonetmail.com)
Last modified: Sat May 2 08:24:03 CEST 2009