R : Copyright 2003, The R Development Core Team Version 1.8.0 (2003-10-08) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. [Previously saved workspace restored] > # derangements.R is a trivial R program written by Chris Evans in 2001 > # copyright is asserted by me, Chris Evans > # Rampton Hospital, Retford, Notts. DN22 0PD Britain > # by telephone at [+44|0] 1777 247242 > # and 'fax at: [+44|0] 1777 247213 > # or http://www.psyctc.org/cgi-bin/mailto.pl?webmaster > > # like all my shabby programming that I put on the www, this is released under the > # "atrribution, share-alike" creative commons licence http://creativecommons.org/licenses/by-sa/1.0/ > # what that says is: > # You are free: > # > # * to copy, distribute, display, and perform the work > # * to make derivative works > # * to make commercial use of the work > # > # Under the following conditions: > # > # Attribution. You must give the original author credit. > # > # Share Alike. If you alter, transform, or build upon this work, you may distribute the resulting work > # only under a license identical to this one. > # > # * For any reuse or distribution, you must make clear to others the license terms of this work. > # * Any of these conditions can be waived if you get permission from the author. > # If you do change or improve on this (shouldn't be hard as I'm no programmer!) I would hugely appreciate > # receiving a copy and your permission to replace this with any improved version with full attribution to you > # Ditto if you translate this into another language, human or computer form. > # > # The theory behind this is fully described in: > # Evans, C., Hughes, J. & Houston, J. (2002) Significance testing the validity of ideographic methods: > # a little derangement goes a long way. British Journal of Mathematical and Statistical Psychology, 55, 385-390. > # and if you contact me, I will endeavour to send you a copy of that if it looks unlikely that you'd find other ways > # of getting it. > # > # This program differs from a program for S+ only in having to declare a function, factorial() which comes with S+ but not > # the version of R on which I'm testing this (1.7.1) and in explicitly declaring tmp at the end of all.derangements() since > # R won't return it to the console (does return it for assignment) if you just end the function with the assignment to tmp > > > factorial <- function(n) { + gamma(n+1) + } > > all.derangements <- function(n){ + cumprob <- prob <- number <- term <- score <- rev(0:n) + for (m in 1:n) { + i <- m+1 + s <- n-m + term[i] <- ((-1)^(m))/(factorial(m)) + } + term[1] <- 1 + for (i in 0:n) { + s <- i+1 + prob[s] <- (sum(term[1:s]))/factorial(n-i) + } + number <- factorial(n)*prob + for (s in 0:n) { + m <- n-s + i <- m+1 + cumprob[i] <- sum(prob[1:i]) + } + tmp <- cbind(score,number,prob,cumprob) + tmp + } > > p.derange.score <- function(score,n){ + if (score > n) stop("Score cannot be greater than n") + if (score == (n-1)) stop ("Score cannot be n-1") + cumprob <- prob <- term <- rev(0:n) + for (m in 1:n) { + i <- m+1 + s <- n-m + term[i] <- ((-1)^(m))/(factorial(m)) + } + term[1] <- 1 + for (i in 0:n) { + s <- i+1 + prob[s] <- (sum(term[1:s]))/factorial(n-i) + } + for (s in 0:n) { + m <- n-s + i <- m+1 + cumprob[i] <- sum(prob[1:i]) + } + cumprob[n+1-score] + } > all.derangements(8) score number prob cumprob [1,] 8 1 2.480159e-05 2.480159e-05 [2,] 7 0 0.000000e+00 2.480159e-05 [3,] 6 28 6.944444e-04 7.192460e-04 [4,] 5 112 2.777778e-03 3.497024e-03 [5,] 4 630 1.562500e-02 1.912202e-02 [6,] 3 2464 6.111111e-02 8.023313e-02 [7,] 2 7420 1.840278e-01 2.642609e-01 [8,] 1 14832 3.678571e-01 6.321181e-01 [9,] 0 14833 3.678819e-01 1.000000e+00 > > p.derange.score(6,8) [1] 0.000719246 > > >