Deming's funnel

From Rosetta Code
Revision as of 22:05, 16 April 2014 by rosettacode>Purple24 (Added Ruby)
Deming's funnel is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

W Edwards Deming was an American statistician and management guru who used physical demonstrations to illuminate his teachings. In one demonstration Deming repeatedly dropped marbles through a funnel at a target, marking where they landed, and observing the resulting pattern. He applied a sequence of "rules" to try to improve performance. In each case the experiment begins with the funnel positioned directly over the target.

  • Rule 1: The funnel remains directly above the target.
  • Rule 2: Adjust the funnel position by shifting the target to compensate after each drop. E.g. If the last drop missed 1 cm east, move the funnel 1 cm to the west of its current position.
  • Rule 3: As rule 2, but first move the funnel back over the target, before making the adjustment. E.g. If the funnel is 2 cm north, and the marble lands 3 cm north, move the funnel 3 cm south of the target.
  • Rule 4: The funnel is moved directly over the last place a marble landed.

Apply the four rules to the set of 50 pseudorandom displacements provided (e.g in the Racket solution) for the dxs and dys. Output: calculate the mean and standard-deviations of the resulting x and y values for each rule.

Note that rules 2, 3, and 4 give successively worse results. Trying to deterministically compensate for a random process is counter-productive, but -- according to Deming -- quite a popular pastime: see the Further Information, below for examples.

Stretch goal 1: Generate fresh pseudorandom data. The radial displacement of the drop from the funnel position is given by a Gaussian distribution (standard deviation is 1.0) and the angle of displacement is uniformly distributed.

Stretch goal 2: Show scatter plots of all four results.


Further information

D

Translation of: Python

<lang d>import std.stdio, std.math, std.algorithm, std.range, std.typecons;

auto mean(T)(in T[] xs) /*pure*/ {

   return xs.sum / xs.length;

}

auto stdDev(T)(in T[] xs) /*pure*/ {

   immutable m = xs.mean;
   return sqrt(xs.map!(x => (x - m) ^^ 2).sum / xs.length);

}

alias TF = double function(in double, in double) pure nothrow;

auto funnel(T)(in T[] dxs, in T[] dys, in TF rule) {

   T x = 0, y = 0;
   immutable(T)[] rxs, rys;
   foreach (const dx, const dy; zip(dxs, dys)) {
       immutable rx = x + dx;
       immutable ry = y + dy;
       x = rule(x, dx);
       y = rule(y, dy);
       rxs ~= rx;
       rys ~= ry;
   }
   return tuple(rxs, rys);

}

void experiment(T)(in string label,

                  in T[] dxs, in T[] dys, in TF rule) {
   //immutable (rxs, rys) = funnel(dxs, dys, rule);
   immutable rs = funnel(dxs, dys, rule);
   label.writeln;
   writefln("Mean x, y:    %.4f, %.4f", rs[0].mean, rs[1].mean);
   writefln("Std dev x, y: %.4f, %.4f", rs[0].stdDev, rs[1].stdDev);
   writeln;

}

void main() {

   immutable dxs = [
   -0.533,  0.270,  0.859, -0.043, -0.205, -0.127, -0.071,  0.275,
    1.251, -0.231, -0.401,  0.269,  0.491,  0.951,  1.150,  0.001,
   -0.382,  0.161,  0.915,  2.080, -2.337,  0.034, -0.126,  0.014,
    0.709,  0.129, -1.093, -0.483, -1.193,  0.020, -0.051,  0.047,
   -0.095,  0.695,  0.340, -0.182,  0.287,  0.213, -0.423, -0.021,
   -0.134,  1.798,  0.021, -1.099, -0.361,  1.636, -1.134,  1.315,
    0.201,  0.034,  0.097, -0.170,  0.054, -0.553, -0.024, -0.181,
   -0.700, -0.361, -0.789,  0.279, -0.174, -0.009, -0.323, -0.658,
    0.348, -0.528,  0.881,  0.021, -0.853,  0.157,  0.648,  1.774,
   -1.043,  0.051,  0.021,  0.247, -0.310,  0.171,  0.000,  0.106,
    0.024, -0.386,  0.962,  0.765, -0.125, -0.289,  0.521,  0.017,
    0.281, -0.749, -0.149, -2.436, -0.909,  0.394, -0.113, -0.598,
    0.443, -0.521, -0.799,  0.087];
   immutable dys = [
    0.136,  0.717,  0.459, -0.225,  1.392,  0.385,  0.121, -0.395,
    0.490, -0.682, -0.065,  0.242, -0.288,  0.658,  0.459,  0.000,
    0.426,  0.205, -0.765, -2.188, -0.742, -0.010,  0.089,  0.208,
    0.585,  0.633, -0.444, -0.351, -1.087,  0.199,  0.701,  0.096,
   -0.025, -0.868,  1.051,  0.157,  0.216,  0.162,  0.249, -0.007,
    0.009,  0.508, -0.790,  0.723,  0.881, -0.508,  0.393, -0.226,
    0.710,  0.038, -0.217,  0.831,  0.480,  0.407,  0.447, -0.295,
    1.126,  0.380,  0.549, -0.445, -0.046,  0.428, -0.074,  0.217,
   -0.822,  0.491,  1.347, -0.141,  1.230, -0.044,  0.079,  0.219,
    0.698,  0.275,  0.056,  0.031,  0.421,  0.064,  0.721,  0.104,
   -0.729,  0.650, -1.103,  0.154, -1.720,  0.051, -0.385,  0.477,
    1.537, -0.901,  0.939, -0.411,  0.341, -0.411,  0.106,  0.224,
   -0.947, -1.424, -0.542, -1.032];
   static assert(dxs.length == dys.length);
   experiment("Rule 1:", dxs, dys, (z, dz) => 0.0);
   experiment("Rule 2:", dxs, dys, (z, dz) => -dz);
   experiment("Rule 3:", dxs, dys, (z, dz) => -(z + dz));
   experiment("Rule 4:", dxs, dys, (z, dz) => z + dz);

}</lang>

Output:
Rule 1:
Mean x, y:    0.0004, 0.0702
Std dev x, y: 0.7153, 0.6462

Rule 2:
Mean x, y:    0.0008, -0.0103
Std dev x, y: 1.0371, 0.8999

Rule 3:
Mean x, y:    0.0438, -0.0063
Std dev x, y: 7.9871, 4.7784

Rule 4:
Mean x, y:    3.1341, 5.4210
Std dev x, y: 1.5874, 3.9304

PARI/GP

This is a work-in-progress.

<lang parigp>drop(drops, rule, rnd)={

 my(v=vector(drops),target=0);
 v[1]=rule(target, 0);
 for(i=2,drops,
   target=rule(target, v[i-1]);
   v[i]=rnd(n)+target
 );
 v

}; R=[-.533-.136*I,.27-.717*I,.859-.459*I,-.043+.225*I,-.205-1.39*I,-.127-.385*I,-.071-.121*I,.275+.395*I,1.25-.490*I,-.231+.682*I,-.401+.0650*I,.269-.242*I,.491+.288*I,.951-.658*I,1.15-.459*I,.001,-.382-.426*I,.161-.205*I,.915+.765*I,2.08+2.19*I,-2.34+.742*I,.034+.0100*I,-.126-.0890*I,.014-.208*I,.709-.585*I,.129-.633*I,-1.09+.444*I,-.483+.351*I,-1.19+1.09*I,.02-.199*I,-.051-.701*I,.047-.0960*I,-.095+.0250*I,.695+.868*I,.34-1.05*I,-.182-.157*I,.287-.216*I,.213-.162*I,-.423-.249*I,-.021+.00700*I,-0.134-.00900*I,1.8-.508*I,.021+.790*I,-1.1-.723*I,-.361-.881*I,1.64+.508*I,-1.13-.393*I,1.32+.226*I,.201-.710*I,.034-.0380*I,.097+.217*I,-.17-.831*I,.054-.480*I,-.553-.407*I,-.024-.447*I,-.181+.295*I,-.7-1.13*I,-.361-.380*I,-.789-.549*I,.279+.445*I,-.174+.0460*I,-.009-.428*I,-.323+.0740*I,-.658-.217*I,.348+.822*I,-.528-.491*I,.881-1.35*I,.021+.141*I,-.853-1.23*I,.157+.0440*I,.648-.0790*I,1.77-.219*I,-1.04-.698*I,.051-.275*I,.021-.0560*I,.247-.0310*I,-.31-.421*I,.171-.0640*I,-.721*I,.106-.104*I,.024+.729*I,-.386-.650*I,.962+1.10*I,.765-.154*I,-.125+1.72*I,-.289-.0510*I,.521+.385*I,.017-.477*I,.281-1.54*I,-.749+.901*I,-.149-.939*I,-2.44+.411*I,-.909-.341*I,.394+.411*I,-.113-.106*I,-.598-.224*I,.443+.947*I,-.521+1.42*I,-.799+.542*I,.087+1.03*I]; rule1(target, result)=0; rule2(target, result)=target-result; rule3(target, result)=-result; rule4(target, result)=result; mean(v)=sum(i=1,#v,v[i])/#v; stdev(v,mu=mean(v))=sqrt(sum(i=1,#v,(v[i]-mu)^2)/#v); main()={

 my(V);
 V=apply(f->drop(100,f,n->R[n]), [rule1, rule2, rule3, rule4]);
 for(i=1,4,
   print("Method #"i);
   print("Means: ", mean(real(V[i])), "\t", mean(imag(V[i])));
   print("StDev: ", stdev(real(V[i])), "\t", stdev(imag(V[i])));
   print()
 )

}</lang>

Python

Translation of: Racket

<lang python>import math

dxs = [-0.533, 0.27, 0.859, -0.043, -0.205, -0.127, -0.071, 0.275, 1.251,

      -0.231, -0.401, 0.269, 0.491, 0.951, 1.15, 0.001, -0.382, 0.161, 0.915,
      2.08, -2.337, 0.034, -0.126, 0.014, 0.709, 0.129, -1.093, -0.483, -1.193, 
      0.02, -0.051, 0.047, -0.095, 0.695, 0.34, -0.182, 0.287, 0.213, -0.423,
      -0.021, -0.134, 1.798, 0.021, -1.099, -0.361, 1.636, -1.134, 1.315, 0.201, 
      0.034, 0.097, -0.17, 0.054, -0.553, -0.024, -0.181, -0.7, -0.361, -0.789,
      0.279, -0.174, -0.009, -0.323, -0.658, 0.348, -0.528, 0.881, 0.021, -0.853,
      0.157, 0.648, 1.774, -1.043, 0.051, 0.021, 0.247, -0.31, 0.171, 0.0, 0.106,
      0.024, -0.386, 0.962, 0.765, -0.125, -0.289, 0.521, 0.017, 0.281, -0.749,
      -0.149, -2.436, -0.909, 0.394, -0.113, -0.598, 0.443, -0.521, -0.799, 
      0.087]

dys = [0.136, 0.717, 0.459, -0.225, 1.392, 0.385, 0.121, -0.395, 0.49, -0.682,

      -0.065, 0.242, -0.288, 0.658, 0.459, 0.0, 0.426, 0.205, -0.765, -2.188, 
      -0.742, -0.01, 0.089, 0.208, 0.585, 0.633, -0.444, -0.351, -1.087, 0.199,
      0.701, 0.096, -0.025, -0.868, 1.051, 0.157, 0.216, 0.162, 0.249, -0.007, 
      0.009, 0.508, -0.79, 0.723, 0.881, -0.508, 0.393, -0.226, 0.71, 0.038, 
      -0.217, 0.831, 0.48, 0.407, 0.447, -0.295, 1.126, 0.38, 0.549, -0.445, 
      -0.046, 0.428, -0.074, 0.217, -0.822, 0.491, 1.347, -0.141, 1.23, -0.044, 
      0.079, 0.219, 0.698, 0.275, 0.056, 0.031, 0.421, 0.064, 0.721, 0.104, 
      -0.729, 0.65, -1.103, 0.154, -1.72, 0.051, -0.385, 0.477, 1.537, -0.901, 
      0.939, -0.411, 0.341, -0.411, 0.106, 0.224, -0.947, -1.424, -0.542, -1.032]

def funnel(dxs, rule):

   x, rxs = 0, []
   for dx in dxs:
       rxs.append(x + dx)
       x = rule(x, dx)
   return rxs

def mean(xs): return sum(xs) / len(xs)

def stddev(xs):

   m = mean(xs)
   return math.sqrt(sum((x-m)**2 for x in xs) / len(xs))

def experiment(label, rule):

   rxs, rys = funnel(dxs, rule), funnel(dys, rule)
   print label
   print 'Mean x, y    : %.4f, %.4f' % (mean(rxs), mean(rys))
   print 'Std dev x, y : %.4f, %.4f' % (stddev(rxs), stddev(rys))
   print


experiment('Rule 1:', lambda z, dz: 0) experiment('Rule 2:', lambda z, dz: -dz) experiment('Rule 3:', lambda z, dz: -(z+dz)) experiment('Rule 4:', lambda z, dz: z+dz)</lang>

Output:
Rule 1:
Mean x, y    : 0.0004, 0.0702
Std dev x, y : 0.7153, 0.6462

Rule 2:
Mean x, y    : 0.0009, -0.0103
Std dev x, y : 1.0371, 0.8999

Rule 3:
Mean x, y    : 0.0439, -0.0063
Std dev x, y : 7.9871, 4.7784

Rule 4:
Mean x, y    : 3.1341, 5.4210
Std dev x, y : 1.5874, 3.9304

Alternative: [Generates pseudo-random data and gives some interpretation.] The funnel experiment is performed in one dimension. The other dimension would act similarly. <lang python>from random import gauss from math import sqrt from pprint import pprint as pp

NMAX=50

def statscreator():

   sum_ = sum2 = n = 0
   def stats(x):
       nonlocal sum_, sum2, n
       sum_ += x
       sum2 += x*x
       n    += 1.0
       return sum_/n, sqrt(sum2/n - sum_*sum_/n/n)
   return stats

def drop(target, sigma=1.0):

   'Drop ball at target'
   return gauss(target, sigma)

def deming(rule, nmax=NMAX):

    Simulate Demings funnel in 1D. 
   
   stats = statscreator()
   target = 0
   for i in range(nmax):
       value = drop(target)
       mean, sdev = stats(value)
       target = rule(target, value)
       if i == nmax - 1:
           return mean, sdev

def d1(target, value):

    Keep Funnel over target. 
   return target


def d2(target, value):

    The new target starts at the center, 0,0 then is adjusted to
   be the previous target _minus_ the offset of the new drop from the
   previous target. 
   
   return -value   # - (target - (target - value)) = - value

def d3(target, value):

    The new target starts at the center, 0,0 then is adjusted to
   be the previous target _minus_ the offset of the new drop from the
   center, 0.0. 
   
   return target - value

def d4(target, value):

    (Dumb). The new target is where it last dropped. 
   
   return value


def printit(rule, trials=5):

   print('\nDeming simulation. %i trials using rule %s:\n %s'
         % (trials, rule.__name__.upper(), rule.__doc__))
   for i in range(trials):
       print('    Mean: %7.3f, Sdev: %7.3f' % deming(rule))


if __name__ == '__main__':

   rcomments = [ (d1, 'Should have smallest deviations ~1.0, and be centered on 0.0'),
                 (d2, 'Should be centred on 0.0 with larger deviations than D1'),
                 (d3, 'Should be centred on 0.0 with larger deviations than D1'),
                 (d4, 'Center wanders all over the place, with deviations to match!'),
               ]
   for rule, comment in rcomments:
       printit(rule)
       print('  %s\n' % comment)</lang>
Output:
Deming simulation. 5 trials using rule D1:
  Keep Funnel over target. 
    Mean:  -0.161, Sdev:   0.942
    Mean:  -0.092, Sdev:   0.924
    Mean:  -0.199, Sdev:   1.079
    Mean:  -0.256, Sdev:   0.820
    Mean:  -0.211, Sdev:   0.971
  Should have smallest deviations ~1.0, and be centered on 0.0


Deming simulation. 5 trials using rule D2:
  The new target starts at the center, 0,0 then is adjusted to
    be the previous target _minus_ the offset of the new drop from the
    previous target. 
    Mean:  -0.067, Sdev:   4.930
    Mean:   0.035, Sdev:   4.859
    Mean:  -0.080, Sdev:   2.575
    Mean:   0.147, Sdev:   4.948
    Mean:   0.050, Sdev:   4.149
  Should be centred on 0.0 with larger deviations than D1


Deming simulation. 5 trials using rule D3:
  The new target starts at the center, 0,0 then is adjusted to
    be the previous target _minus_ the offset of the new drop from the
    center, 0.0. 
    Mean:   0.006, Sdev:   1.425
    Mean:  -0.039, Sdev:   1.436
    Mean:   0.030, Sdev:   1.305
    Mean:   0.009, Sdev:   1.419
    Mean:   0.001, Sdev:   1.479
  Should be centred on 0.0 with larger deviations than D1


Deming simulation. 5 trials using rule D4:
  (Dumb). The new target is where it last dropped. 
    Mean:   5.252, Sdev:   2.839
    Mean:   1.403, Sdev:   3.073
    Mean:  -1.525, Sdev:   3.650
    Mean:   3.844, Sdev:   2.715
    Mean:  -7.697, Sdev:   3.715
  Center wanders all over the place, with deviations to match!

Racket

The stretch solutions can be obtained by uncommenting radii etc. (delete the 4 semi-colons) to generate fresh data, and scatter-plots can be obtained by deleting the #; . <lang racket>#lang racket (require math/distributions math/statistics plot)

(define dxs '(-0.533 0.270 0.859 -0.043 -0.205 -0.127 -0.071 0.275 1.251 -0.231

             -0.401 0.269 0.491 0.951 1.150 0.001 -0.382 0.161 0.915 2.080 -2.337 
             0.034 -0.126 0.014 0.709 0.129 -1.093 -0.483 -1.193 0.020 -0.051
             0.047 -0.095 0.695 0.340 -0.182 0.287 0.213 -0.423 -0.021 -0.134 1.798
             0.021 -1.099 -0.361 1.636 -1.134 1.315 0.201 0.034 0.097 -0.170 0.054 
             -0.553 -0.024 -0.181 -0.700 -0.361 -0.789 0.279 -0.174 -0.009 -0.323
             -0.658 0.348 -0.528 0.881 0.021 -0.853 0.157 0.648 1.774 -1.043 0.051 
             0.021 0.247 -0.310 0.171 0.000 0.106 0.024 -0.386 0.962 0.765 -0.125 
             -0.289 0.521 0.017 0.281 -0.749 -0.149 -2.436 -0.909 0.394 -0.113 -0.598
             0.443 -0.521 -0.799 0.087))

(define dys '(0.136 0.717 0.459 -0.225 1.392 0.385 0.121 -0.395 0.490 -0.682 -0.065

             0.242 -0.288 0.658 0.459 0.000 0.426 0.205 -0.765 -2.188 -0.742 -0.010 
             0.089 0.208 0.585 0.633 -0.444 -0.351 -1.087 0.199 0.701 0.096 -0.025 
             -0.868 1.051 0.157 0.216 0.162 0.249 -0.007 0.009 0.508 -0.790 0.723
             0.881 -0.508 0.393 -0.226 0.710 0.038 -0.217 0.831 0.480 0.407 0.447
             -0.295 1.126 0.380 0.549 -0.445 -0.046 0.428 -0.074 0.217 -0.822 0.491 
             1.347 -0.141 1.230 -0.044 0.079 0.219 0.698 0.275 0.056 0.031 0.421 0.064
             0.721 0.104 -0.729 0.650 -1.103 0.154 -1.720 0.051 -0.385 0.477 1.537 
             -0.901 0.939 -0.411 0.341 -0.411 0.106 0.224 -0.947 -1.424 -0.542 -1.032))
(define radii (map abs (sample (normal-dist 0 1) 100)))
(define angles (sample (uniform-dist (- pi) pi) 100))
(define dxs (map (λ (r theta) (* r (cos theta))) radii angles))
(define dys (map (λ (r theta) (* r (sin theta))) radii angles))

(define (funnel dxs rule)

 (let ([x 0])
   (for/fold ([rxs null])
     ([dx dxs])
     (let ([rx (+ x dx)])
       (set! x (rule x dx))
       (cons rx rxs)))))

(define (experiment label rule)

 (define (p s) (real->decimal-string s 4))
 (let ([rxs (funnel dxs rule)]
       [rys (funnel dys rule)])
   (displayln label)
   (printf "Mean x, y   : ~a, ~a\n" (p (mean rxs)) (p (mean rys)))
   (printf "Std dev x, y: ~a, ~a\n\n" (p (stddev rxs)) (p (stddev rys)))
   #;(plot (points (map vector rxs rys)
         #:x-min -15 #:x-max 15 #:y-min -15 #:y-max 15))))

(experiment "Rule 1:" (λ (z dz) 0)) (experiment "Rule 2:" (λ (z dz) (- dz))) (experiment "Rule 3:" (λ (z dz) (- (+ z dz)))) (experiment "Rule 4:" (λ (z dz) (+ z dz))) </lang>

Output:
Rule 1:
Mean x, y   : 0.0004, 0.0702
Std dev x, y: 0.7153, 0.6462

Rule 2:
Mean x, y   : 0.0009, -0.0103
Std dev x, y: 1.0371, 0.8999

Rule 3:
Mean x, y   : 0.0439, -0.0063
Std dev x, y: 7.9871, 4.7784

Rule 4:
Mean x, y   : 3.1341, 5.4210
Std dev x, y: 1.5874, 3.9304

Ruby

Translation of: Python

<lang ruby>def funnel(dxs, &rule)

 x, rxs = 0, []
 for dx in dxs
   rxs << (x + dx)
   x = rule[x, dx]
 end
 rxs

end

def mean(xs) xs.inject(:+) / xs.size end

def stddev(xs)

 m = mean(xs)
 Math.sqrt(xs.inject(0.0){|sum,x| sum + (x-m)**2} / xs.size)

end

def experiment(label, dxs, dys, &rule)

 rxs, rys = funnel(dxs, &rule), funnel(dys, &rule)
 puts label
 puts 'Mean x, y    : %7.4f, %7.4f' % [mean(rxs), mean(rys)]
 puts 'Std dev x, y : %7.4f, %7.4f' % [stddev(rxs), stddev(rys)]
 puts

end

dxs = [ -0.533, 0.270, 0.859, -0.043, -0.205, -0.127, -0.071, 0.275,

        1.251, -0.231, -0.401,  0.269,  0.491,  0.951,  1.150,  0.001,
       -0.382,  0.161,  0.915,  2.080, -2.337,  0.034, -0.126,  0.014,
        0.709,  0.129, -1.093, -0.483, -1.193,  0.020, -0.051,  0.047,
       -0.095,  0.695,  0.340, -0.182,  0.287,  0.213, -0.423, -0.021,
       -0.134,  1.798,  0.021, -1.099, -0.361,  1.636, -1.134,  1.315,
        0.201,  0.034,  0.097, -0.170,  0.054, -0.553, -0.024, -0.181,
       -0.700, -0.361, -0.789,  0.279, -0.174, -0.009, -0.323, -0.658,
        0.348, -0.528,  0.881,  0.021, -0.853,  0.157,  0.648,  1.774,
       -1.043,  0.051,  0.021,  0.247, -0.310,  0.171,  0.000,  0.106,
        0.024, -0.386,  0.962,  0.765, -0.125, -0.289,  0.521,  0.017,
        0.281, -0.749, -0.149, -2.436, -0.909,  0.394, -0.113, -0.598,
        0.443, -0.521, -0.799,  0.087]

dys = [ 0.136, 0.717, 0.459, -0.225, 1.392, 0.385, 0.121, -0.395,

        0.490, -0.682, -0.065,  0.242, -0.288,  0.658,  0.459,  0.000,
        0.426,  0.205, -0.765, -2.188, -0.742, -0.010,  0.089,  0.208,
        0.585,  0.633, -0.444, -0.351, -1.087,  0.199,  0.701,  0.096,
       -0.025, -0.868,  1.051,  0.157,  0.216,  0.162,  0.249, -0.007,
        0.009,  0.508, -0.790,  0.723,  0.881, -0.508,  0.393, -0.226,
        0.710,  0.038, -0.217,  0.831,  0.480,  0.407,  0.447, -0.295,
        1.126,  0.380,  0.549, -0.445, -0.046,  0.428, -0.074,  0.217,
       -0.822,  0.491,  1.347, -0.141,  1.230, -0.044,  0.079,  0.219,
        0.698,  0.275,  0.056,  0.031,  0.421,  0.064,  0.721,  0.104,
       -0.729,  0.650, -1.103,  0.154, -1.720,  0.051, -0.385,  0.477,
        1.537, -0.901,  0.939, -0.411,  0.341, -0.411,  0.106,  0.224,
       -0.947, -1.424, -0.542, -1.032]

experiment('Rule 1:', dxs, dys) {|z, dz| 0} experiment('Rule 2:', dxs, dys) {|z, dz| -dz} experiment('Rule 3:', dxs, dys) {|z, dz| -(z+dz)} experiment('Rule 4:', dxs, dys) {|z, dz| z+dz}</lang>

Output:
Rule 1:
Mean x, y    :  0.0004,  0.0702
Std dev x, y :  0.7153,  0.6462

Rule 2:
Mean x, y    :  0.0009, -0.0103
Std dev x, y :  1.0371,  0.8999

Rule 3:
Mean x, y    :  0.0439, -0.0063
Std dev x, y :  7.9871,  4.7784

Rule 4:
Mean x, y    :  3.1341,  5.4210
Std dev x, y :  1.5874,  3.9304