#! /usr/bin/R
sink(file="/dev/null") # to stop getting the Hmisc announcement text
# 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)
}


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]
}

# p.derange.score(6,8)

sink() # switch output back to stdout

tag(HTML)
    tag(HEAD)
        tag(TITLE)
            cat("Probability of a derangement score r from n")
        untag(TITLE)
    untag(HEAD)

    lf(2)

    tag(BODY, bgcolor = "lime")
        lf(2)
        tag(center)
        cat("<H1>Probability of a derangement score r from n</H1>")

        comments("Let's start with some testing")
        prob <- 0
        score <- as.numeric(scanText(formData$score))
        n <- as.numeric(scanText(formData$n))
        round <- as.numeric(scanText(formData$round))
        if ((trunc(score) != score) || (score < 0)) {
        cat("<H2>Score must be a positive integer, go back and try again!</H2>")
           prob <- 1
        }
        if ((trunc(n) != n) || (n < 0) || (n > 290)){
           cat("<H2>n must be a positive integer &lt;201, go back and try again!</H2>")
           prob <- 1
         }
        if (score > n) {
           cat("<H2>Score must be smaller than or equal to n, go back and try again!</H2>")
           prob <- 1
         }
        if (trunc(round) != round) {
           cat("<H2>Score must be an integer, go back and try again!</H2>")
           prob <- 1
        }
        if ((round < 0) | (round > 8)) {
           cat("<H2>Program only offers rounding to between zero and eight significant digits.  Go back and try again!</H2>")
           prob <- 1
         }
        if (prob) {
           cat("<H1>Go back and try again!</H1>")
           cat("CGI script written by Chris Evans using David Firth's excellent GGIwithR package.  ")
           linkto("Contact me if something isn't working right ...", "http://www.psyctc.org/cgi-bin/mailto.pl?chris") ; br()
         }
        lf(2)
        comments("end of testing")
        if (!prob) {
           comments("Got into results section")
           cat("<TABLE BORDER='1' CELLSPACING='1' CELLPADDING='3' WIDTH='600' BGCOLOR='lime')")
           cat("<TR>")
           cat("<TH BGCOLOR='white' COLSPAN='2'>")
           cat("<B>Your input</B>")
           cat("</TH>")
           cat("</TR>")
           cat("<tr>")
           cat("<td>")
              now <- system('date +%Y-%m-%d,%T', intern=TRUE)
              host <- system("echo $REMOTE_ADDR", intern=TRUE)
              cat("Request from: ")
           cat("</td>")
           cat("<td>")
              cat(host, " at", now,"<BR>")
           cat("</td>")
       cat("</tr>")
       cat("<tr>")
           cat("<td>")
              cat("Score: ")
           cat("</td>")
           cat("<td>")
              cat(score)
           cat("</td>")
       cat("</tr>")
       cat("<tr>")
           cat("<td>")
              cat("n: ")
           cat("</td>")
           cat("<td>")
              cat(n)
           cat("</td>")
       cat("</tr>")
       cat("<tr>")
           cat("<td>")
              cat("Significant digits wanted: ")
           cat("</td>")
           cat("<td>")
              cat(round)
            cat("</td>")
       cat("</tr>")
       cat("<tr>")
            cat("<th bgcolor='white' COLSPAN='2'>")
            cat("<b>Results</b>")
            cat("</th>")
            cat("</tr>")
            cat("<tr>")
            cat("<td COLSPAN='2'><b><center>")
            p <- p.derange.score(score,n)
            p <- signif(p,round)
            cat("<I>p</I> = ",p)
            cat("</center></b>")
            cat("</td>")
            cat("</tr>")
            cat("<tr>")
            cat("<th bgcolor='white' COLSPAN='2'>")
            cat("<b>Explanation</b>")
            cat("</th>")
            cat("</tr>")
            cat("<tr>")
            cat("<td COLSPAN='2'>That's the probability of scoring",score,"<I>or more</I> from ",n,"possible matches is",p,".  If this is &lt;0.05 then you are entitled on the conventional logic of significance testing, to say that the test has shown statistically significant evidence that something better than chance is at work here.  You can set a more stringent criterion if you like, of course.")
           cat("<P>It's as simple as that: if you achieved a score with a probability lower than your preset criterion probability (alpha) then you have statistical support that something non-random is at work in your matching the objects back to their origins.  Whether that will impress people will depend on whet her you have eliminated trite ways in which you might have achieved this.")
        cat("<P>If you are using this method and publishing the results or finding them of use to you, then ")
       linkto("I'd love to hear about what you're doing.", "http://www.psyctc.org/cgi-bin/mailto.pl?chris") ; br()


           cat("</td>")
        cat("</tr>")
        cat("<tr>")
        cat("<th bgcolor='white' COLSPAN='2'>")
        cat("<b>Technicalities</b>")
        cat("</th>")
        cat("</tr>")
        cat("<tr>")
        cat("<td COLSPAN='2'>")
        cat("Calculation done in R by a CGI script written by Chris Evans using David Firth's excellent GGIwithR package.")
        cat("  Neither myself, David Firth or anyone in the R team accept responsibility for the results or consequences of their use, the maths of this is fairly easy and spelled out in my publication on the topic and R is very reliable but it's always possible for things to go wrong so think through what you see here carefully before doing anything serious with it.")
     linkto("Contact me if something isn't working right ...", "http://www.psyctc.org/cgi-bin/mailto.pl?chris") ; br()

     lf()

     cat("<p>")
        cat("Output produced at ", date())
     cat("</p>")

     lf()
     cat("</td>")
     cat("</tr>")
     cat("</table>")
     }
  untag(BODY)
untag(HTML)
lf()
sink(paste("/home/xychris0/R/CGIwithR-log/derange.p.R.",as.character(now),sep=""))
cat("program = derange.p.R ; host = ",host, "; ", now, "; score = ",score,"; n = ",n,"; p = ",p,"; round = ",round,"\n")
sink()