Lite om resampling, simulering, sannolikhetsproblem etc.

Författad av Håkan Kjellerstrand (hakank@bonetmail.com).

Här är lite spridda saker om simulering av sannolikhetsproblem, statistik etc.

Resampling

Resampling är en metod att göra simuleringar på olika typer av sannolikhets-/ statistiska problem. Andra begrepp som betyder ungefär samma sak som resampling eller använder resampling

Dessa metoder tillhör dock mer den formella statistiken än vad som avses med resampling här. T.ex. används begreppet Monte Carlo-simulering av vissa författare endast för att lösa rena matematiska problem, men resampling används för "alla" typer av problem. Förvirrande? Visst, men så är det ibland:-)

www.resample.com

http://www.resample.com/

Här finns exempel, böcker etc hur man "resamplar" olika typer av statistiska/ sannolikhetsproblem. Mycket inspirerande läsning.

http://www.resample.com/content/text/index.shtml: En onlinebok som är en väldigt grundlig genomgång om principerna bakom resampling och statistik/sannolikhetsanalys överhuvudtaget.

Programmet ("demo standalone") Resampling Stat. Finns även versioner för Excel och Matlab:

En Javaversion med ungefär samma funktionalitet och syntax som standaloneversionen av Resampling Stats finns på http://www.ipd.uka.de/~prechelt/sw/ (rubriken Resampling Statistics).

Ingen av ovanstående har jag kollat in så väldigt mycket eftersom jag skriver mina saker i Java och framförallt (det GNU:ade) statistiksystemet R, http://cran.r-project.org/. Jag har ännu inte gjort någon systematisk library för resampling, men det är inte så många speciella funktioner man behöver för att få samma funktionalitet som i boken. (Se nedan för lite exempel i R.)

Andra intressanta länkar på resample.com, mestadels introduktionsartiklar av olika slag och grad:

Julian Simon

Julian Simon, som är en av grundarna av Resampling Stats Inc. och skrev ovanstående online-bok ("Resampling: The New Statistics"), har också skrivit en massa andra saker också, främst inom ekonomi (han var ekonom). Han är mycket kritisk till traditionell statistisk analys som han tycker är alldeles för krånglig och formelbaserad.

http://www.juliansimon.com/writings/ finns flera av hans andra texter online. Tyvärr är de alla i ASCII, men det går att läsa efter viss tids tillvänjning.

De mest intressanta vad gäller simulering/resampling är:

Andra intressanta sajter/personer

En av de saker jag speciellt intresserar mig för - förutom sannolikhetsproblem - är slump/tillfällighet/skrock (coincidence, superstition). Det finns en del roliga sajter om detta.

Några exempel på simuleringar, främst i R

I det följande gör jag lite simuleringar av olika typer av problem/gåtor. Jag kommer att blanda engelsk och svensk text eftersom jag har snott en del problemformuleringar från utlandet (och det är i princip kopierat rakt av från mina anteckningar).

Notera att för mer enkla problem finns det rätt enkla formler att tillämpa, men - som Julian Simon skriver - för mer avancerade problem kan det dock vara det bästa om man inte är expert i matematiskt statistik; om inte annat för att kontrollera att ens analys är korrekt. Att programmera (skriva simuleringskoden) ett sådant problem är för övrigt ett bra sätt att hitta konstigheter eller undantag i formuleringen av problemet.

Jag förklarar inte R här. Mer info om detta utsökta program finns på http://cran.r-project.org/

Eller maila mig (hakank@bonetmail.com) om det är några konstigheter.

Födelsedagsproblemet

Se till exempel denna text (sid 23) för formulering av problemet.

Här är exempel på simulering dels i Resampling Stats dels i R.

Simulering i Resampling Stats

Från texten ovan:
What is the probability that two or more people among a roomful of 25 have the
same birthday?

The birthday problem is amazingly simple with RESAMPLING STATS. 
First, GENERATE 25 numbers between 1 and 365 (the number of days in a typical 
year) and store the results in a. 

   GENERATE 25 1,365 a 

Next, using the MULTIPLES command, check to see whether any two or more people
have the same birthday by searching for duplicate (or triplicate, 
quadruplicate, etc.) numbers and put the result in b. 

   MULTIPLES a >=2 b 

We want to keep track of the individual results b and keep them in our special
 score  vector scrboard. (Vector names can have up to 8 characters.) 

  SCORE b scrboard 

Let s review the program to see where the REPEAT, END  loop  goes. 

REPEAT 1000
  GENERATE 25 1,365 a 
  MULTIPLES a >=2 b 
  SCORE b scrboard
END
COUNT scrboard >=1 k 
DIVIDE k 1000 prob 
PRINT prob

Simulering i R

Några korta förklaringar av viktiga funktioner i R:

Här är 25 tal från intervallet 1 till 365 för att simulera 25 personers födelsedagar. Notera konstrukten sum(...)>=1 .
> sum(table(sample(1:365, 25,replace=T))>=2)
[1] 1

> sum(sapply(1:1000, function(x) sum(table(sample(1:365, 25,replace=T))>=2))>=1)/1000
[1] 0.563

Detta visar alltså att sannolikheten för att det ska finnas två personer med samma födelsedag bland 25 personer är något över 50%.

Vad är sannolikheten för två näraliggande födelsedagar, dvs samma datum eller de närmast intilliggande? (Se t.ex. http://ourworld.compuserve.com/homepages/rajm/tscoin.htm).

Sannolikheten att någon ska fylla år samma dag eller en dag före eller en dag efter är i häradet 13-14:

> sum(sapply(1:1000, function(x) sum(diff(sort(sample(1:365,14,replace=T)),1)<=1))>0)/1000
[1] 0.523

> sum(sapply(1:1000, function(x) sum(diff(sort(sample(1:365,13,replace=T)),1)<=1))>0)/1000
[1] 0.487

Hur många krävs för att de ska fylla år samma månad med 50% chans? Någonstans i häradet 4-5:

> sum(sapply(1:10000, function(x) sum(table(sample(1:12,4,replace=T))>1))>0)/10000
[1] 0.4258
> sum(sapply(1:10000, function(x) sum(table(sample(1:12,5,replace=T))>1))>0)/10000
[1] 0.6212

Vacker granne

(Ett eget litet problem jag haft:) Jag och min vackra granne tittar ut x resp y gånger dagligen. Problemet jag simulerar är hur många gånger kommer vi att titta på varandra. Jag förutsätter just nu att vi tittar ut en enda sekund. En senare fråga blir att undersöka om vi tittar ut längre och möjligen olika långa stunder varje gång.

Simulering

Varje dag har 24*60*60 sekunder = 86400 sekunder. Jag förenklar det och tar bara minuterna, dvs 24*60 = 1440 minuter. Dessa representeras av en vektor av 1:1400 .

För tillfället hårdkodar jag värdena hur ofta vi tittar ut och förutsätter att vi tittar ut lika många gånger (det blir lite enklare då), t.ex. n=200 gånger per dag:


# R funktion
 granne.sim <- function(n=200) {
    hakank <- sample(1:1400, n, replace=F)
    vacker.granne <- sample(1:1400, n, replace=F)

    # summerar hur många gånger det är samma värde i de två vektorerna
    sum(hakank==vacker.granne)
 }

Här simuleras 10 olika värden för n, där vi ser att för n = 1400 kommer vi att titta på varandra cirka 1 gång per dag.

> sapply(seq(1,1400,length=10), function(n) c(n, sum(sapply(1:1000, function(z) sum(granne.sim(n)))/1000)))

     [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]
[1,]    1 156.4444 311.8889 467.3333 622.7778 778.2222 933.6667 1089.111
[2,]    0   0.1070   0.2300   0.3310   0.4300   0.5890   0.6770    0.779
         [,9]   [,10]
[1,] 1244.556 1400.00
[2,]    0.901    1.04

Nästa version är att vi kikar olika många gånger ut.

granne.mult.sim <- function(hakank.times=200,granne.times=250) {
     secs <- 24*60*60
     hakank <- sample(secs, hakank.times, replace=F)
     granne <- sample(secs, granne.times, replace=F)
     length(intersect(hakank,granne))
}

> sum(sapply(1:1000,function(x) granne.mult.sim(200,250))>0)/1000
[1] 0.422

Dvs om vi ser ut 200 resp 250 gånger per dag, kommer vi att - under en dag - se varandra med cirka 42% chans. Vilket är

> 1/0.422
[1] 2.369668
cirka lite mer än en gång varannan dag.

Testar med 25 gånger per dag:

> 1/(sum(sapply(1:10000,function(x) granne.mult.sim(25,25)))/10000)
[1] 125

Tolkning: om vi tittar ut 25 gånger om dagen i en sekund så kommer vi att titta på varandra en gång per 125 dagar.

Här försöker jag även få med att vi inte tittar ut lika länge varje gång. Denna version kan råka ut för att vi tittar ut "flera gånger samma gång" eftersom jag - för varje slumpmässigt startsekund - lägger till slumpmässigt antal sekunder (mellan 1 och *.dur sekunder)

granne.mult2.sim <- function(hakank.times=200,granne.times=250,hakank.dur=5,granne.dur=5) {
     secs <- 24*60*60
     hakank <- unlist(sapply(1:hakank.times, function(x) {
                xx<-sample(secs, 1, replace=F);
                c(xx:(xx+sample(hakank.dur,1)))
               }))
     granne <- unlist(sapply(1:granne.times, function(x) {
                 xx<-sample(secs, 1, replace=F);
                 c(xx:(xx+sample(granne.dur,1)))
               }))

     length(intersect(hakank,granne))
     
}

> sum(sapply(1:100, function(x) granne.mult2.sim(25,25,10,10))>0)
[1] 13

Dvs: i snitt .13 gemensamma tittar per dag.

> describe(sapply(1:1000, function(x) granne.mult2.sim(25,25,10,10))/1000)
sapply(1:1000, function(x) granne.mult2.sim(25, 25, 10, 10))/1000 
       n  missing   unique     Mean      .05      .10      .25      .50 
    1000        0       11 0.000294    0.000    0.000    0.000    0.000 
     .75      .90      .95 
   0.000    0.000    0.003 

          0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.010 0.015
Frequency   913    13    18    27    12     8     1     3     3     1     1
%            91     1     2     3     1     1     0     0     0     0     0


> describe(sapply(1:1000, function(x) granne.mult2.sim(25,25,10,10)))
sapply(1:1000, function(x) granne.mult2.sim(25, 25, 10, 10)) 
      n missing  unique    Mean     .05     .10     .25     .50     .75     .90 
   1000       0      12   0.409    0.00    0.00    0.00    0.00    0.00    1.00 
    .95 
   3.05 

            0  1 10 11  2  3  4  5  6  7  8  9
Frequency 894 17  1  3 27 12  9 15 12  6  2  2
%          89  2  0  0  3  1  1  2  1  1  0  0

Samlarkortproblemet (Coupon collector's problem)

Hur länge måste man slå en tärning innan alla sidor har visats? Se t.ex. http://www.maa.org/features/mathchat/mathchat_1_7_99.html för en formulering och diskussioner. Från denna sida:

OLD CHALLENGE (Steve Jabloner). How many times do you think you need to roll a
normal die to be 90% sure that each of the six faces has appeared at least once? Why?

ANSWER. It takes 23 rolls, as Al Zimmerman discovered in a computer experiment 
rolling 5000 dice until 90% of them had shown each face. In a less accurate 
experiment with just 300 dice, Eric Brahinsky wrongly concluded that 22 rolls 
provide a 90.7% probability. 

Jag testar detta:

dice.sim <- function() {   
    n <- 6
    total <- c(rep(0,n))
    num.sides <- FALSE
    num.throws <- 0
    while (!num.sides) {
       dice<- sample(1:n,1)
       total[dice] <- TRUE
       num.sides <- sum(total)==n
       num.throws <- num.throws + 1
    }
    # returnera antal kast som krävdes
    num.throws
}

Gör nu en simulering av 5000 sådana körningar eftersom man gjorde det i referensen som nämns i artikeln:
> describe(sapply(1:5000,function (x) dice.sim()))
      n missing  unique    Mean     .05     .10     .25     .50     .75     .90 
   5000       0      44   14.75       8       8      10      13      18      23 
    .95 
     26 

lowest :  6  7  8  9 10, highest: 46 47 48 49 50 

Dvs det krävs i medelvärde cirka 15 kast. Om man vill vara säker till 90% krävs det cirka 23 kast, vilket alltså överensstämmer med texten ovan. För att vara säker till 95% krävs cirka 26 kast.

Notera här att man bör ta värdena man får via en simulering med en viss nypa skepsis (och salt). Ju större simuleringar man gör, desto säkrare kan man vara på att värdena är korrekta (givet att programmet är korrekt:-).

Detta problem är samma som "samlarkortproblemet" (coupon collector problem, dvs hur många paket frukostflingor måste man köpa för att få samtliga bilder i en samlarserie. Jag antar att bilderna har lika stor chans att finnas i ett paket.

collector.sim <- function(n=6) {   
    total <- c(rep(0,n))  # det vore trevligt med en hashfunktion i R!
    all.collected <- FALSE
    num <- 0
    while (!all.collected) {
       samp<- sample(n,1)
       total[samp] <- TRUE
       all.collected <- sum(total)==n
       num <- num + 1
    }
    num
}

Hur många paket måste man köpa för att vara säker till 97.5% att få hela serien. Vi kollar olika samlarserier där antalet bilder i serien är 2 till 20.

# detta tar en stund
> sapply(2:20, function(x) c(x,quantile(sapply(1:1000, function(z) collector.sim(x)),c(.025,.975))))
      [,1] [,2] [,3] [,4]   [,5]   [,6] [,7] [,8]   [,9] [,10] [,11] [,12]
         2    3    4    5  6.000  7.000    8    9 10.000    11    12 13.00
2.5%     2    3    4    5  7.000  8.000   10   12 14.000    16    19 20.00
97.5%    7   12   19   26 30.025 36.025   43   51 58.025    63    68 80.05
      [,13] [,14]   [,15]   [,16]   [,17]   [,18]   [,19]
         14 15.00  16.000  17.000  18.000  19.000  20.000
2.5%     24 26.00  29.000  32.000  33.975  36.000  39.000
97.5%    81 92.05 103.025 108.025 116.000 118.025 135.025

Dvs finns det 20 bilder i serien krävs det cirka 135 inköp.

(Ett tips är att börja med små simuleringar, t.ex. på vara 10 eller 100 loopar (i sapply) för att få en känsla för hur resultatet kommer att bli. Man kan då rätt snabbt se om det är några konstigheter. Större simuleringar ger större tillförlitlighet men de tar också längre (ibland mycket längre) tid att köra.)

Med 100 bilder:

> sapply(100, function(x) c(x,quantile(sapply(1:100, function(z) collector.sim(x)),c(.025,.975))))
       [,1]
      100.0
2.5%  303.2
97.5% 791.3

Om man har 99 samlarbilder av 100, hur många paket måste man köpa för att få den 100:e?

> sapply(99:100, function(x) c(x,quantile(sapply(1:1000, function(z) collector.sim(x)),c(.025,.975))))
        [,1]   [,2]
       99.00 100.00
2.5%  337.95 345.00
97.5% 787.05 817.05

Se även http://rec-puzzles.org/sol.pl/probability/coupon

Sultanens hemgift (Sultans Dowry Problem)

Se http://mathworld.wolfram.com/SultansDowryProblem.html:
A sultan has granted a commoner a chance to marry one of his n daughters. The 
commoner will be presented with the daughters one at a time and, when each 
daughter is presented, the commoner will be told the daughter's dowry (which 
is fixed in advance). Upon being presented with a daughter, the commoner must 
immediately decide whether to accept or reject her (he is not allowed to 
return to a previously rejected daughter). However, the sultan will allow the 
marriage to take place only if the commoner picks the daughter with the 
overall highest dowry. Then what is the commoner's best strategy, assuming he 
knows nothing about the distribution of dowries (B. Elbows)?

Since the commoner knows nothing about the distribution of the dowries, the 
best strategy is to wait until a certain number x of daughters have been 
presented, then pick the highest dowry thereafter. The exact number to skip is
determined by the condition that the odds that the highest dowry has already 
been seen is just greater than the odds that it remains to be seen and that if 
it is seen it will be picked. 
....
The problem is most commonly stated with n = 100 daughters, which gives the 
result that the commoner should wait until he has seen 37 of the daughters, 
then pick the first daughter with a dowry that is bigger than any preceding 
one. With this strategy, his odds of choosing the daughter with the highest 
dowry are surprisingly high: about 37% 
(B. Elbows; Honsberger 1979, pp. 104-110, Mosteller 1987).

Jag följer denna strategi. Notera att man ska välja det tal (efter skip.numret) som är större än samtliga föregående tal.

sultans.dowry.sim <- function(skip=.50, n=100) {
   s<-sample(1:(n*10),n,replace=F)
   max.s <- max(s)
   skip.pos<-floor(n*skip)
   for (i in skip.pos:n) {
     if (s[i] > max(s[1:(i-1)])) {
       return(s[i] == max.s)
     } 
   }
   return(s[n] == max.s)
}

Här går man igenom .1:.9 för att se vilket värde som är högst:

> sapply(seq(.1,1,length=10),function(z) c(z,mean(sapply(1:1000,function(x) sultans.dowry.sim(z)))))
      [,1]  [,2]  [,3] [,4]  [,5]  [,6] [,7]  [,8]  [,9] [,10]
[1,] 0.100 0.200 0.300 0.40 0.500 0.600 0.70 0.800 0.900 1.000
[2,] 0.227 0.326 0.372 0.37 0.359 0.327 0.25 0.198 0.092 0.006

# Finjustera några gånger (.2 - .5)
> sapply(seq(.2,.5,,length=10),function(z) c(z,mean(sapply(1:1000,function(x) sultans.dowry.sim(z)))))
      [,1]      [,2]      [,3]  [,4]      [,5]      [,6]  [,7]      [,8]
[1,] 0.200 0.2333333 0.2666667 0.300 0.3333333 0.3666667 0.400 0.4333333
[2,] 0.331 0.3250000 0.3460000 0.358 0.4000000 0.3760000 0.385 0.3620000
          [,9] [,10]
[1,] 0.4666667 0.500
[2,] 0.3560000 0.357

# intervallet (.3 - .4)
> sapply(seq(.3,.4,,length=10),function(z) c(z,mean(sapply(1:1000,function(x) sultans.dowry.sim(z)))))
      [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
[1,] 0.300 0.3111111 0.3222222 0.3333333 0.3444444 0.3555556 0.3666667
[2,] 0.365 0.3470000 0.3670000 0.3600000 0.3740000 0.4100000 0.3740000
          [,8]      [,9] [,10]
[1,] 0.3777778 0.3888889 0.400
[2,] 0.3730000 0.3640000 0.371

# intervallet (.35 - .4)
> sapply(seq(.35,.4,,length=10),function(z) c(z,mean(sapply(1:1000,function(x) sultans.dowry.sim(z)))))
      [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
[1,] 0.350 0.3555556 0.3611111 0.3666667 0.3722222 0.3777778 0.3833333
[2,] 0.371 0.3860000 0.3960000 0.3720000 0.3680000 0.3770000 0.3650000
          [,8]      [,9] [,10]
[1,] 0.3888889 0.3944444 0.400
[2,] 0.3480000 0.3610000 0.368

Och nu orkar jag inte finjustera längre, utan beslutar mig för att det verkar som om gränsen .377 är den bästa strategin.

Se även http://www.juliansimon.com/puzzles/50chall3.txt för en diskussion om detta problem.

Monty Hall

Det finns en massa bra sajter på nätet om detta kända problem: Här kollar jag de två strategierna:
monty1 <- function(change=0) {
  doors  <- 1:3
  prize  <- sample(doors, 1) # the price door
  picked <- sample(doors, 1) # pick this door

  # Monty opens one door that is not my picked door and not the prize door
  dif     <- setdiff(doors, c(picked, prize))

  # NOTE: a single value as first argument to sample is taken as size!
  if (length(dif) == 1) {
     opens   <- dif
  } else {
     opens   <- sample(dif, 1)
  }
  left    <- setdiff(doors, c(picked, opens)) # it is the one left
  if (change == 1) { # should we change?
    picked <- left
  }
  win <- sum(picked==prize)
  win
}

# testar strategin att inte byta
> sum(sapply(1:10000, function(x) monty1(0)))
[1] 3347

# testar strategin att bya
> sum(sapply(1:10000, function(x) monty1(1)))
[1] 6626

Dvs om man inte byter dörr vinner man cirka 1/3 av gånger, men byter man så vinner man 2/3 av gångerna.

Några kö-problem

Se http://mathforum.org/epigone/comp.soft-sys.math.mathematica/frekorprox/6hpeka$dd8@smc.vnet.net.
Meeting Under the Clock (This problem is posed by Julian Simon(1994))

Two persons agree to arrive at the two clock sometime between 1 pm and 2 pm 
and to stay for 20 minutes. What is the probability that they will be there
at the same time?

Simulering av detta:

under.the.clock.sim <- function() {
  p1 <- sample(60,1)
  p1 <- p1:(p1+20) # if >= 60 then the time is past 2 pm
  p2 <- sample(60,1)
  p2 <- p2:(p2+20)
  length(intersect(p1,p2))>0
}
under.the.clock.sim()

> sum(resamp.sim(10000,under.the.clock.sim))/10000
[1] 0.5686

Functionen resamp.sim som används definieras i all sin enkelhet så här:
resamp.sim <- function(trials=100, f, par=NULL) {
   if (!is.null(par)) {
     sapply(1:trials, function(x) f(par))
   } else {
     sapply(1:trials, function(x) f())
   }
}

En enklare variant: Om differansen mellan p1 och p2 <= 20 så kommer de att träffas.
under.the.clock2.sim <- function() {
  p1 <- sample(60,1)
  p2 <- sample(60,1)
  abs(p1-p2)<=20
}

> describe(sapply(1:100,function(x) sum(resamp.sim(1000,under.the.clock2.sim))/1000))
sapply(1:100, function(x) sum(resamp.sim(1000, under.the.clock2.sim))/1000) 
      n missing  unique    Mean     .05     .10     .25     .50     .75     .90 
    100       0      51  0.5676  0.5440  0.5489  0.5570  0.5670  0.5790  0.5891 
    .95 
 0.5910 

lowest : 0.526 0.534 0.541 0.542 0.543, highest: 0.590 0.591 0.592 0.598 0.601 
Nästa problem är från http://groups.google.com/groups?dq=&hl=en&lr=&ie=UTF-8&c2coff=1&safe=off&th=f0a097f2a90a0280
Lets say that 4 people agree to meet between 3:00 P.M. and 4:00 P.M.. 
One man can wait for 10 minutes. Another can wait for 15 minutes. 
Another can wait for 20 Minutes, and another still can wait for 30 
minutes. What is the probability that all 4 will meet?
Jag förutsätter här:
waiting4.sim <- function() {
    a <- sample(0:60,1)
    b <- sample(0:60,1)
    c <- sample(0:60,1)
    d <- sample(0:60,1)

    a <- a:(a+10)
    b <- b:(b+15)
    c <- c:(c+20)
    d <- d:(d+30)

    length(intersect(a,intersect(b,intersect(c,d))))>0
    
}
waiting4.sim()
> sum(sapply(1:10000, function(x) sum(waiting4.sim())))/10000
[1] 0.073

Dvs ungefär var 14:e gång kommer de att träffas:
> 1/0.073
[1] 13.69863
Ett tredje exempel: http://exn.ca/Stories/1997/04/09/01.asphttp://exn.ca/Stories/1997/04/09/01.asp. Om Robert Matthews:
Murphy's law of queues - The line next to you will usually finish first - is 
also rooted in simple math.

If you're waiting behind a person with a two-months supply of groceries, it's 
hardly a surprise if you get through more slowly than your neighbours. But 
what about joining a line that's identical in length from the ones on either 
side of you. Is it all in your mind when you are the last one to get through?

Actually...NO.
Over time, it's true that on average the lines will move more or less at the 
same rate, says Matthews. Each will suffer random delays when for example the 
cashier runs out of change or a customer's cheque or bank card is declined. 
But when we're lined up at any one particular time, we don't care about 
averages, we care only about that one instance. In these cases, chances of you
picking the fastest moving line is one out of the number of lines there are in
the store. If there are 10 lines, you have a 1 in 10 chance (or 10%) of 
picking the fastest line. Even if you're only concerned about beating the 
lines on either side of you, you only have a 1 in 3 chance that you will. 
You're more likely to be in a slower line.
En simulering av detta. Några förutsättningar: Chansen att jag ska ställa mig i den snabbaste kön är:
queue2.sim <- function(queues=5,persons=10) {
  tot.time <- rep(0,queues)
  # for each queue
  for (i in 1:queues) {
     # time for each of the persons in the queue
     for (j in 1: persons) {
        tot.time[i] <- tot.time[i] + sample(1:5,1)
     }
  }
  tot.time
}
Min kö är #1. Hur ofta kommer denna kö att:

Vi slumpar nu antalet kunder till att vara mellan 4 och 15. Låt oss säga att det motsvarar min kalkyl hur snabb/långsam kön är (dvs med 10 personer kan jag - rätt eller fel - göra en bedömning att den är snabb som 4 eller långsam som 15). Fortfarande gäller betjäningen på mellan 1:5 minuter, dock.

queue3.sim <- function(queues=5) {
  tot.time <- rep(0,queues)

  # for each queue
  for (i in 1:queues) {
     # time for each of the persons in the queue
     persons <- sample(4:15,1)
     for (j in 1: persons) {
        tot.time[i] <- tot.time[i] + sample(1:5,1)
     }
  }
  tot.time
}
Min kö är #1. Hur ofta kommer denna kö att:

Dvs, det verkar ingen roll vad jag tror om vilken som är den snabbaste/bästa kön. Det blir i stort sett samma sannolikheter som med exakt 10 personer.

OK, nästa steg är att jag faktiskt väljer kön med minst antal personer. Jag slumpar antalet för min kö och antar sedan att de andra köerna har fler personer.

queue4.sim <- function(queues=5) {
  tot.time <- rep(0,queues)

  persons <- rep(0,queues)
  persons[1] <- sample(3:10,1)
  for (q in 2:queues) {
    persons[q] <- persons[1] + sample(1:7,1)
  }

  # for each queue
  for (q in 1:queues) {
     # time for each of the persons in the queue
     for (j in 1:persons[q]) {
        tot.time[q] <- tot.time[q] + sample(1:5,1)
     }
  }
  tot.time
}
Min kö är #1. Hur ofta kommer denna kö att Dvs det är rätt stor chans att om jag jag kommer först om jag väljer den kortaste kön.

Litet problem i (engelsk) Lotto

Se http://groups.google.com/groups?hl=en&lr=&ie=UTF-8&c2coff=1&safe=off&th=64dbcce10d08e875&rnum=47:
One day I was asked 

What is the probability that no two numbers are consecutive in
a lottery draw?

Information:

National Lottery in UK is to draw 6 different numbers from 1 to 49.
so, for example,

(a)   1  4  7   12  19  44    - no consecutive numbers
(b)   3  6  17  18  44  46    - 17 and 18 are consecutive 
(c)   1  2  3   17  29  49    - 1, 2 and 3 are consecutive

We are asking the probability that class (a) occurs. Hope this is clear.

Observation shows that it is NOT a small number (actually near half-and-half).
......

A Monte Carlo simulation experiment produced 50,558 sets of six numbers 
between 1 and 49 with none consecutive out of 100,000 trials giving an 
estimated probability of .50558.
Simulering av detta i R:
> sum(sapply(1:10000, function(x) {ss<-sort(sample(1:49,6));ss2<-ss[2:6]-ss[1:5];sum(ss2!=1)==5}))/10000
[1] 0.5011

> describe(sapply(1:100, function(z) sum(sapply(1:1000, function(x) {ss<-sort(sample(1:49,6));ss2<-ss[2:6]-ss[1:5];sum(ss2!=1)==5}))/1000)) 
sapply(1:100, function(z) sum(sapply(1:1000, function(x) {     ss <- sort(sample(1:49, 6))     ss2 <- ss[2:6] - ss[1:5]     sum(ss2 != 1) == 5 }))/1000) 
      n missing  unique    Mean     .05     .10     .25     .50     .75     .90 
    100       0      51  0.5052  0.4809  0.4840  0.4948  0.5060  0.5160  0.5261 
    .95 
 0.5320 

lowest : 0.449 0.474 0.475 0.479 0.481, highest: 0.531 0.532 0.533 0.535 0.540 


> sum(sapply(1:100000, function(x) {ss<-sort(sample(1:49,6));ss2<-ss[2:6]-ss[1:5];sum(ss2!=1)==5}))/100000
[1] 0.50441

Frustration patience

Från Grinstead & Snells "Introduction of Probability" http://www.dartmouth.edu/~Echance/teaching_aids/books_articles/probability_book/book.pdf, sid 86: Frustraton patience.

Problemet är följande: Blanda en kortlek. För varje kort räkna 1,2,3, osv till 13, sedan 1,2,3 osv. Om det någon gång kommer upp samma valör som det man räknar så har man förlorat.

Hur ofta kommer ett sådant sammanfallande upp? Svaret ska bli 4.

> mean(sapply(1:10000,function(x) sum(sample(rep(1:13,4),replace=F)==rep(1:13,4))))
[1] 3.9812

Hur stor är chansen att man vinner, dvs att det inte kommer upp något sammanfallande?

> sum(sapply(1:10000, function(x) sum(sample(rep(1:13,4),replace=F)==rep(1:13,4))==0))/10000
[1] 0.0179
Korrekt svar är:
> 1/exp(1)**4
[1] 0.01831564

Knutar på ett snöre

Robert Matthews har skrivit en artikel om knutar på ett snöre Knotted rope: a topological example of Murphy's Law ("If rope can become knotted, it will do").

Simulering:
Matthews pratar om en random walk i 3 dimensioner.

Mitt försök till simulering är följande:

Här tillåter jag endast att gå ett steg i någon riktning, dvs antingen -1 eller 1 i någon av riktningarna x,y,z. Jag ignorerar problemet att det kanske krävs N steg innan man verkligen kan få en knut. (Vi kan ju tänka oss att varje punkt är sammanbunden på ett rep, eller är en orm etc...) Jag ignorerar också - och det är kanske allvarligare - att det krävs någon form av "rörelse" runt punkten för att det ska bli en knut (dvs att det inte räcker att man kommer tillbaka i samma punkt, man måste gå på ett speciellt sätt "in" och "ut" till/från punkten. Men i alla fall... Och möjligen skulle man lägga till en möjlighet att man inte rör sig alls.

# Notera listhanteringen.
murphy.knot.sim <- function() {
  pos <- c(0,0,0)
  positions <- list(c(0,0,0))
  step <- 0
  while (step < 300) {
    axis <- sample(1:3,1)
    dir  <- sample(c(-1,1), 1)
    pos[axis] <- pos[axis]+dir
    if (is.element(list(pos), positions)) {
      return(step) 
    }
    positions <- c(positions, list(pos))
    step <- step + 1
  }

  step
}
Medelantalet steg är cirka 5.3. Notera att det högsta antalet är bara 34!
> describe(ss<-sapply(1:10000, function(x) murphy.knot.sim()))
ss <- sapply(1:10000, function(x) murphy.knot.sim()) 
      n missing  unique    Mean     .05     .10     .25     .50     .75     .90 
  10000       0      32   5.268       1       1       2       4       7      11 
    .95 
     14 

lowest :  1  2  3  4  5, highest: 28 29 32 33 34 


# Standardavvikelsen:
> sd(ss)
[1] 4.363381
Fråga: Hur påverkas medellängden av hur många dimensioner det är?
murphy.knot2.sim <- function(n=3) {
  pos <- rep(0,n) # c(0,0,0)
  positions <- list(rep(0,n)) # list(c(0,0,0))
  step <- 0
  while (step < 300) {
    axis <- sample(1:n,1)
    dir  <- sample(c(-1,1), 1) # with a possibility of not "moving around"
    pos[axis] <- pos[axis]+dir
    if (is.element(list(pos), positions)) {
      return(step) 
    }
    positions <- c(positions, list(pos))
    step <- step + 1
  }

  step
}
Testar med dimensionerna 2:100 (bara med 100 simuleringar per dimension):
> sapply(2:10, function(z) c(z,mean(sapply(1:100, function(x) murphy.knot2.sim(z)))))
     [,1] [,2] [,3] [,4]  [,5]  [,6]  [,7]  [,8]  [,9]
[1,] 2.00 3.00 4.00 5.00  6.00  7.00  8.00  9.00 10.00
[2,] 3.24 5.18 7.89 8.95 11.25 13.51 14.42 19.24 20.66
Se för övrigt avsnittet om Random Walks i Grinstead & Snell. Där förklaras en hel del av ovanstående.

Cluster watching

Från http://www.numberwatch.co.uk/of_birthdays_and_clusters.htm:

As we gratefully witness the dying fall of yet another silly season, it is 
interesting to note how often birthdays provide the basis of the most popular 
fallacies. In August 2001 a classic of the genre appeared when the New 
Scientist announced that researchers in Scotland had established that 
anorexics are likely to be have been June babies. The team studied 446 women 
who had been diagnosed as anorexic and observed that 30% more than average 
were born in June. As the monthly average of births is about 37, we deduce 
that the June number must have been 48. At first sight this looks like a 
significant result (at least by epidemiological standards) since the 
probability of getting 48 or more in a random month is about 3%. But that is 
not what they are doing! They are making twelve such selections and then 
picking the biggest. Application of the theory of the statistics of extremes 
tells us that the probability of the largest of twelve such selections being 
48 or greater is 30%, which is not significant even by epidemiological 
standards.

Medelvärdet födelser per månad är c:a 37.

> 446/12
[1] 37.16667

Vi simulerar 446 födelser (12 olika månader) och ser om det är någon månad som är 30% över 37.

Värdet som vi ska jämföra med (som är 30% över medelvärdet) är:

> 446/12*1.3
[1] 48.31667

> sum(sapply(1:1000, function(x) sum(table(sample(1:12,446,replace=T))>=37.16667*1.3)>0))/1000
[1] 0.307

Dvs det finns - som det mycket riktigt står i artikeln - cirka 30% chans att en sådan fördelning uppstår av ren slump.

Om vi tittar på max() värdet så fördelar det sig enligt följande:

> describe(sapply(1:10000, function(x) max(table(sample(1:12,446,replace=T)))))
sapply(1:10000, function(x) max(table(sample(1:12, 446, replace = T)))) 
      n missing  unique    Mean     .05     .10     .25     .50     .75     .90 
  10000       0      25   47.36      43      44      45      47      49      52 
    .95 
     53 

lowest : 40 41 42 43 44, highest: 60 61 62 64 67 

Dvs: medelvärdet av max-värdena är cirka 47.3 (endast cirka 1 från det värde vi kontrollerar mot).

Hur många skulle vara födda i samma månad för att det skulle anses som signifikant (oavsett problement med post hoc-hypoteser)?

I describe() ovan ser vi att för .95 (dvs 5% one tailed) hittar vi värdet 53, vilket alltså är det sökta värdet för denna siginifikansnivå.

Om vi tar 99%-signifikansnivå blir det:

> quantile(sapply(1:10000, function(x) max(table(sample(1:12,446,replace=T)))),c(.01,.99))
 1% 99% 
 42  57 
dvs det krävs minst 57 födda i samma månad för att vi ska börja höja på ögonbrynen.

Lite senare i ovanstående artikel:

...
At the time, ninety people had died of vCJD. The probability of two of them 
having the same birthday is 0.999993848, i.e. a certainty. In fact there was 
an evens chance of three of them having the same birthday. Likewise, if we 
divided the UK up into 1000 areas of equal population, we would find by the 
same calculation that the probability of two of the ninety coming from the 
same area is about 0.98. Yet the original Queniborough two were claimed to be 
statistically significant and the hunt was begun for more 'linked pairs'.

För 90 personer, hur stor är chansen att minst två har samma födelsedag?

> sum(sapply(1:1000, function(x) sum(table(sample(1:365, 90,replace=T))>1)>0))/1000
[1] 1
Av 90 personer som bor i ett av 1000 områden, hur stor är chansen är att två eller flera bor i samma?
> sum(sapply(1:10000, function(x) sum(table(sample(1:1000, 90,replace=T))>1)>0))/10000
[1] 0.9833
Once the micro-cluster was identified, of course, the search was widened to 
include people who 'came from the same part of Leicestershire' to expand it to 
five. This is known in the trade as the method of the Texas sharp-shooter, who
sprays the side of a barn with bullets and then draws a target round the most
prominent cluster.

Om rekord (extreme values)

Från http://www.numberwatch.co.uk/record.htm. Se även http://www.numberwatch.co.uk/extreme_value_fallacy.htm samt ovan nämnda anorexi-fall-fall.
Let us invent a stationary process and create a little allegory. In the 
eighteenth century an eccentric amateur astronomer looks through his low-power
telescope each night, pointing roughly in the same direction, and counts the 
number of stars in its field of vision. He is keen to establish a record so 
that he can apply for entry in the contemporary equivalent of the Guinness 
Book of Records. The conditions are such that the average number of stars he 
sees is 100 and the average scatter (standard deviation) is plus or minus ten 
stars.

We can model this process by generating numbers from a normal distribution 
(for the cognoscenti the function rnorm(N,100,10) in Mathcad).

Our astronomer takes one reading each night and only writes down the current 
record, the highest value so far, so after a hundred nights he gets something 
like....
....
By now we have reached the tenth generation. It is the late twentieth century 
and the age of rationality has come to an end. The new scion of the family, 
realising that it has obtained no benefit from all these years of effort, 
publishes a paper warning of the increase in star density that is occurring 
and calculating all sorts of disasters that arise from the increased 
gravitation. He is appointed to a professorial chair in a new university and 
founds a new Department of Celestial Change. The cause is taken up by the 
Deputy Prime Minister and large sums of money are diverted from other research 
areas to a new programme for the study of celestial change. Disaster stories 
are regularly fed to the media, who sell more papers and advertising time.

A few old-fashioned scientists point out that the results are completely 
explained by the well-known statistics of extremes, which predict that the 
largest value will increase as the logarithm of the number of observations. 
There is no money in that, so they are ignored and the bandwagon rolls on.
Här simuleras 20 nätters observationer.
record.sim <- function(n=20) {
  record <- 0
  records <- c()
  for (i in 1:n) {
    this.night<-rnorm(100,100,10)  # mean: 100, sd:10
    record <- max(record, this.night)
    records <- c(records,record)
  }
  records
}

> record.sim(20)
 [1] 126.0962 126.0962 132.8600 132.8600 132.8600 132.8600 132.8600 132.8600
 [9] 132.8600 132.8600 132.8600 132.8600 132.8600 133.7859 133.7859 133.7859
[17] 133.7859 133.7859 133.7859 133.7859


> describe(record.sim(100))
record.sim() 
      n missing  unique    Mean 
    100       0       7   135.7 

          127.4 128.3 128.9 131.6 134.2 137.9 140.8
Frequency     3    15     1    12    15    22    32
%             3    15     1    12    15    22    32

Om PK experiment

Från http://www.fourmilab.ch/rpkp/experiments/statistics.html:
The data are right before our eyes, and the probability we calculated is 
correct, but we asked the wrong question, and in doing so fell into a trap 
littered with the bones of many a naïve researcher. Note the order in which 
we did things. We ran the experiment, examined the data, found something 
seemingly odd in it, then calculated the probability of that particular oddity
appearing by chance. We asked the question, "What is the probability of 
`999999' appearing in a 1000 digit random sequence?" and got the answer 
"less than one in a thousand", a result most people would consider significant.
But since we calculated the probability after seeing the data, in fact we were
asking the question "What is the chance that `999999' appears in a 1000 digit 
random sequence which contains one occurrence of `999999'?. The answer to that 
question is, of course, "certainty".
Tänkte nu simulera en sådan test, dvs bland 1000 slumpmässiga tal finna en speciell sekvens.
check.seq <- function(seq, pattern) {
  seq.len <- length(seq)
  pat.len <- length(pattern)
  sum(sapply(1:(seq.len-pat.len-1), function(x) all(seq[x:(x+pat.len-1)]==pattern)))
}

> sapply(1:100, function(x) check.seq(sample(0:9,1000,replace=T),1:4))
  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [38] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
 [75] 0 0 0 3 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

> sapply(1:100, function(x) check.seq(sample(0:9,1000,replace=T),1:4))
  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [38] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
 [75] 0 0 0 3 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Kollar nu förekomsten av 999999 (rep(9,6)) i en 1000-sekvens med 100-sim.
> sum(sapply(1:100, function(x) check.seq(sample(0:9,1000,replace=T),rep(9,6)))>0)
[1] 1
Med en simulering av 1000 körningar:
> sum(sapply(1:1000, function(x) check.seq(sample(0:9,1000,replace=T),rep(9,6)))>0)
[1] 0
Testar samma mönster på en 10000-sträng 100-sim.
> sapply(1:100, function(x) check.seq(sample(0:9,10000,replace=T),rep(9,6)))
  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
 [38] 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [75] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Simuleringar av några av Mostellers sannolikhetsproblem

Från http://www.juliansimon.com/puzzles/50chall1.txt:
Problem 2: Successive Wins
        
  To encourage Elmer's promising tennis career, his 
  father offers him a prize if he wins (at least) two 
  tennis sets in a row in a three-set series to be played 
  with his father and the club champion alternately:  
  father-champion-father or champion-father-champion, 
  according to Elmer's choice.  The champion is a better 
  player than Elmer's father.  Which series should Elmer 
  choose?
Simulering i R:
 
succ.wins2 <- function() {
   win.loss <- c("W","L")
   # notera att man kan ha prob= i sample
   cc <- sample(win.loss, 2, prob=c(2,8),replace=T)
   ff <- sample(win.loss, 2, prob=c(4,6),replace=T) 

   res.cfc <- 0
   res.fcf <- 0

   # C - F - C
   if ( ((cc[1] == "W") & (ff[1] =="W")) | ((ff[1]=="W") & (cc[2]=="W"))) {
       res.cfc <- 1
   }

   # F - C - F
   if ( ((ff[1] == "W") & (cc[1] =="W")) | ((cc[1]=="W") & (ff[2]=="W"))) {
       res.fcf <- 1
   }

   c(res.cfc, res.fcf)

}

> ss<-sapply(1:10000, function(x) succ.wins2());sum(ss[1,]);sum(ss[2,])
[1] 1473
[1] 1299
Ett annat problem:
   Problem 1.5: The Red and the Black (from Gardner, 1983, p. 42))
  "Shuffle a packet of four cards - two red, two black - 
   and deal them face down in a row.  Two cards are picked 
   at random, say by placing a penny on each.  What is the 
   probability that those two cards are the same color?"
Simulering:
> sum(sapply(1:1000, function(x) {a<-sample(rep(0:1,each=2),2);sum(a)==0|sum(a)==2}))/1000
[1] 0.347


cards1.sim <- function(x) {
   # Red: 1 Black: 2
   array1<- c(1,1,2,2)
   array2<- sample(array1)
   a<-array2[1]
   b<-array2[2]
   c <- a-b
   c
}

> sum(resamp.sim(1000,cards1.sim)==0)/1000
[1] 0.332

> describe(sapply(1:100, function(x) sum(resamp.sim(1000,cards1.sim)==0)/1000))
sapply(1:100, function(x) sum(resamp.sim(1000, cards1.sim) == 0)/1000) 
      n missing  unique    Mean     .05     .10     .25     .50     .75     .90 
    100       0      48   0.335  0.3119  0.3160  0.3220  0.3350  0.3462  0.3540 
    .95 
 0.3590 

lowest : 0.296 0.301 0.306 0.308 0.309, highest: 0.360 0.371 0.374 0.378 0.380 

En annan spelkortssimulering

Från en annan Julian Simons bok: http://63.126.124.88/content/teaching/philosophy/part3/chapIII-6.txt:
  You are standing in the warehouse of a playing-card factory
  that has been hit by a tornado.  Cards are scattered everywhere,
  some not yet wrapped and others ripped out of their packages.
  The factory makes a variety of decks - for poker without a joker,
  poker with a joker, and pinochle; magician's decks; decks made of
  paper and others of plastic; cards of various sizes; and so on.
    Two hours from now a friend will join you for a game of
  near-poker with these cards. Each hand will be chosen as randomly
  as possible from the huge heap of cards, and then burned. What
  odds should you attach to getting the combination two-of-a-kind -
  two cards of different or the same suit but of the same number or
  picture - in a five-card draw?
Dvs det är frågan om en kortlek med återläggning. (Frågan är om det är exakt ett par eller åtminstone ett par, jag antar det förra.) sannolikheten för exakt ett par:
> sum(sapply(1:10000, function(x) sum(table(sample(1:13,5,replace=T))==2)==1))/10000 
[1] 0.4637
Triss:
> sum(sapply(1:10000, function(x) sum(table(sample(1:13,5,replace=T))==3)==1))/10000 
[1] 0.0506

Back to my homepage
Created by Hakan Kjellerstrand hakank@bonetmail.com