Här är lite spridda saker om simulering av sannolikhetsproblem, statistik etc.
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:-)
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, 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.
På 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:
Amazonlänk till denna bok.
(Mosteller är en av de stora inom detta område och skrev 1989 tillsammans med Persi Diaconis "Methods for studying coincidences" som fortfarande citeras när man pratar om slump, tillfälligheter och skrock. Tyvärr har jag inte hittat detta paper någonstans. [Att Persi Diaconis tillika är en duktig magiker, liksom för övrigt en annan av mina skeptiska husgudar Martin Gardner, är extra skoj.]).
Detta är ett inte helt färdigt manuskript om resampling på lite mer abstrakt (filosofisk) nivå, men har en hel del exempel på simuleringar.
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.
Denna sajt har många intressanta småartiklar om hur statistik (miss)brukas i media och annorstädes. Det finns även länkar till introduktionsmaterial om sannolikhetslära, t.ex. en mycket trevlig bok (PDF) http://www.dartmouth.edu/~chance/teaching_aids/books_articles/probability_book/book.html som (också) betonar simuleringens roll när man ska lära sig och arbeta sannolikhetslära. Det finns exempel i Mathematica, True BASIC, Maple och samt Java applets (inklusive kod).
Speciellt kan rekommenderas "Statistical Issues in ESP Research" av Ray Hyman (en av de stora inom debunkningsbranchen).
John Allen Paulos sajt. Han är en av de mest aktiva inom området att försöka få folk (oss!) att förstå statistik och är mycket kritisk till hur media använder (missbrukar) statistik. Speciellt boken "Innumeracy" är.' mycket populär. Den rekommenderas varmt. Länkar till boken:
Paulos har också en kolumn Who's counting? i där han skriver om relaterade ämnen.
[Lite skoj är att Paulos tidigare böcker "I think, therefore I laugh" och "Mathematics and Humor" köpte jag i början/mitten av 80-talet när jag var som mest intresserad av vitsar. I dessa böcker gör Paulos intressanta och skojiga analyser av vitsar och humor med matematiska övertoner. De påverkade mig rätt mycket...]
Robert Matthews har gjort en massa roliga (i betydelsen intressanta) analyser om vardagsproblem (Murphys lagar) t.ex. om och varför en tappad smörgås tenderar att komma med smörsidan ner, varför köerna vid sidan om "alltid" går snabbare, varför ett snöre tenderar att få en knut och har analyserat begreppet coincidence.
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.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.
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
Några korta förklaringar av viktiga funktioner i R:
> 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
(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.369668cirka 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
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
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.
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.
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.69863Ett 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:
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:
> sum(sapply(1:1000, function(x) order(queue2.sim())[1]==1))/1000 [1] 0.224
> sum(sapply(1:1000, function(x) 1 %in% order(queue2.sim())[1:3]))/1000 [1] 0.642
> sum(sapply(1:1000, function(x) 1 %in% order(queue2.sim())[3:5]))/1000 [1] 0.541
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:
> sum(sapply(1:1000, function(x) order(queue3.sim())[1]==1))/1000 [1] 0.21
> sum(sapply(1:1000, function(x) 1 %in% order(queue3.sim())[1:3]))/1000 [1] 0.622
> sum(sapply(1:1000, function(x) 1 %in% order(queue3.sim())[3:5]))/1000 [1] 0.578
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
> sum(sapply(1:1000, function(x) order(queue4.sim())[1]==1))/1000 [1] 0.79
> sum(sapply(1:1000, function(x) 1 %in% order(queue4.sim())[1:3]))/1000 [1] 0.995
> sum(sapply(1:1000, function(x) 1 %in% order(queue4.sim())[3:5]))/1000 1] 0.039
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
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.0179Korrekt svar är:
> 1/exp(1)**4 [1] 0.01831564
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.
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 57dvs 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] 1Av 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.
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
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] 1Med 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] 0Testar 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
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
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.4637Triss:
> sum(sapply(1:10000, function(x) sum(table(sample(1:13,5,replace=T))==3)==1))/10000 [1] 0.0506