simulation.scm

(define time-since-last
  (lambda (α t)
    (- 1 (exp (/ (- t) α)))))

(define ρ 200)

(define wait-time
  (lambda (α β x)
    (* ρ (expt x (- α 1)) (expt (- 1 x) (- β 1)))))

(define average
  (lambda (values)
    (/ (fold-left + 0 values) (length values))))

(define variance
  (lambda (values)
    (let* ((average (average values))
	   (squared-sum (fold-left + 0 (map (lambda [x] (* (- x average) (- x average)))
					    values))))
      (/ squared-sum (length values)))))

(define standard-error  (lambda (variance num-samples) (sqrt (/ variance num-samples))))

(define run-monte-carlo
  (lambda (num-samples f)
    (let ((sample-values (map (lambda (i) (f))
			      (iota num-samples))))
      (cons (average sample-values) (apply max sample-values)))))


;; average waiting time 
(define average-waiting-time
  (lambda (customer)
    (run-monte-carlo 200000
		     (lambda () (wait-time (car customer) (cdr customer) (random 1.0))))))

(define yellow-customer (cons 2 5))
(define red-customer (cons 2 2))
(define blue-customer (cons 5 1))

(define display-customer-times
  (lambda (color customer)
    (let ((waiting-times (average-waiting-time customer)))
      (display "Customer type: ") (display color) (newline)
      (display "average waiting time: ") (display (car waiting-times)) (newline)
      (display "max waiting time: ") (display (cdr waiting-times)) (newline)
      (display "δ: ") (display (- (cdr waiting-times) (car waiting-times))) (newline))))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Printing customer times ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(display-customer-times 'yellow yellow-customer)
(display-customer-times 'red red-customer)
(display-customer-times 'blue blue-customer)


;; queue length calculations

;; average interarrival time

(define max-interarrival-time 600)

(define interarrival-time-probability
  (run-monte-carlo 200000
		   (lambda ()
		     (time-since-last 100 (random max-interarrival-time)))))

(define simulation-time 600)

;; sample simulation code
(define queue-length
  (lambda (average-waiting-time average-interarrival-time)
    (fold-left (lambda (queue-data t)
		 (let ((queue-length (car queue-data))
		       (time-since (cadr queue-data))
		       (processing-time (caddr queue-data)))
		   (let* ((data1 (if (> processing-time average-waiting-time)
				     (list (- queue-length 1) (+ 1 time-since) 0)
				     (list queue-length (+ 1 time-since) (+ 1 processing-time))))
			  (queue-length (car data1))
			  (time-since (cadr data1))
			  (processing-time (caddr data1)))
		     (if (> time-since average-interarrival-time)
			 (list (+ queue-length 1) 0 processing-time)
			 data1))))
	       ;; queue-length time-since processing time
	       (list 0 0 0)
	       (iota simulation-time))))


;; (display (queue-length (car average-waiting-time)
;; 		       (car average-interarrival-time)))
;; (newline)