Abs Puzzle

If someone asked you to implement abs function (say, for ints), how would you do that? A simple experiment shows that almost everyone comes up with something like that, whether it’s F#, C#, Scala or anything else:

1: let abs x = if x < 0 then -x else x

The obvious question - do the answers match your expectations? That depends…

absx

[Originally posted here.]

It’s quite obvious when you think about the corner cases - here it’s the minimal integer value:

-(-2147483648) = -2147483648

So, in Java/Scala the absolute value can be negative:

Note that if the argument is equal to the value of Integer.MIN_VALUE, the most negative representable int value, the result is that same value, which is negative.

In .NET you’ll get an exception - personally, I think that’s a better behaviour for this function:
OverflowException - value equals Int32.MinValue

When I started to think about different examples for our future meetup (collecting all kinds of crashes on the way - from ‘CLR detected an invalid program’ to ‘unexpected exception’), I couldn’t resist but look at what happens when you move some checks to the type level.

Let’s define a refinement type for the natural numbers in F* and make sure it works (you can try the examples in your browser):

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
    
    
module Nums

type nat = x:int{x >= 0}

let good_x: nat = 42
let bad_x: nat = -42
input(6,17-6,20) : Error : Expected an expression of type: x_3:int{GTE x_3 0} but got
(op_Minus(42)): x_357_1:int{Eq2 int int x_357_1 (Minus (42))} Type checking failed: Nums

Now it’s turn for our abs function:

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
val abs: int -> nat
let abs x = if x < 0 then -x

let _ = abs 2147483647  
let _ = abs (-2147483647) 

// both following statements give Syntax error:                                
let _ = abs 2147483648
let _ = abs (-2147483648)

The first error is kind of expected, because the literal 2147483648 is too large for int. However, I’d think the second one should be ok, and after a quick look at the F* lexer not sure why they both fail ("Allow &lt;max_int+1&gt; to parse as min_int"). Anyway, there’s another interesting case:

 1: 
 2:   
      
let _ = abs (2147483647 + 1)
// And the result is...                                                         
Verified module: Nums

A bit strange, isn’t it? If you ask Scala or F# what 2147483647 + 1 is, you get -2147483648, that means abs (2147483647 + 1) can’t be of nat type! Here’s a bunch of checks for ‘+1’ with different refinements:

 1: 
 2: 
 3: 
 4: 
    
    
    
val check: x:nat -> y:int{y > x}
let check x = x + 1 

let _ = check (2147483647 + 1)
Error : Let-bound variable op_Addition(2147483647, 1) escapes its scope in type
x_8_1:int{GT x_8_1 x_11_3}; insert an explicit type annotation. (expected type none)
Type checking failed: Nums    
 1: 
    
    
    
val check: x:nat -> y:int{y = (x+1)} 
Error : Let-bound variable op_Addition(2147483647, 1) escapes its scope in type
x_6_1:int{Eq2 int int x_6_1 (Add (x_9_5) (1))}; insert an explicit type annotation.
(expected type none) Type checking failed: Nums
 1: 
    
val check: x:nat -> y:int{y > 0} 
Verified module: Nums                                                            
 1: 
    
    
    
val check: x:nat -> y:int{y < 0 \/ y > x}
Error : Let-bound variable op_Addition(2147483647, 1) escapes its scope in type
x_8_1:int{((LT x_8_1 0) || (GT x_8_1 x_11_3))}; insert an explicit type annotation.
(expected type none) Type checking failed: Nums

All the cases above succeeded for the max value:

let _ = check 2147483647

I also tried these examples on my old-slightly-working F* 0.7 alpha Mono build, where the results seem to be more reasonable:

 1: 
    
    
let _ = abs 2147483648 
Error : Expected an expression of type: int but got (2147483648.000000): float
Type checking failed: Nums
 1: 
    
    
    
    
    
let _ = abs (-2147483648) 
 starting backend
 starting derefinement
 starting expanding types
 finished with expanding types
Verified module: Nums                                                    

Maybe one day I’ll finally setup a VM with Windows to take a look at IL, but for now the conclusion is to always be on the lookout.

P.S. Pex can find the failing case (for some reason, at least in online version this example works only with C# and VB):

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
    
    
    
    
using System;
using System.Diagnostics.Contracts;

public class Program {
  public static void Puzzle(int x) {
    var y = x < 0 ? -x : x;
    Contract.Assert(y >= 0);
  }
}
x            | Output/Exception  | Error Message
-----------------------------------------------------------
0            |                   |
int.MinValue | ContractException | Assertion failed: y >= 0    

 
 

FsCheck and co: Testing Excel Financial Functions

Excel Financial Functions library is a .NET library written by Luca Bolognese that provides the full set of financial functions from Excel. It can be used from both F# and C# as well as from other .NET languages. The main goal for the library is compatibility with Excel, by providing the same functions, with the same behaviour.

Excel-Funcs-Logo

[Originally posted here.]

The library was recently moved to GitHub and its new home is here. The goal was to turn an archive with sources into a community project, including the build scripts, tests, documentation and a fancy homepage. And thanks to FAKE and FSharp.Formatting that was much more fun than I had expected.

So, what did we want to see in the new version?
- the library should be as close as possible to the original;
- the tests, running automatically during the build;
- mono compatibility;
- a couple of pages with nicely formatted docs;
- (ideally) no new unnecessary dependencies.

The library itself doesn’t depend on the Excel interop libraries (available only in .NET), that means it can be used on mono right out of the box. The only problem was in checking the results against Excel: there was an amazing test suite in form of console app, which I wanted to leave. But how to modify it to match the new requirements?

Choosing the Framework

The first step, of course, is to ask the community ^_^ - what frameworks people like most and use, in what cases, how that worked for them. Thanks again to @mausch, @kurt2001, @brandewinder, @foxyjackfox, @sergey_tihon, @davefancher, @c4fsharp, @kitlovesfsharp and @_____c for the feedback!

But then I realized that there’s no particular need in any additional frameworks: the functions are rather simple, you can compare floats or dates with NUnit, why add new dependencies? A single FsUnit-like shouldEqual is more than enough (yes, couldn’t resist):

 1: 
 2: 
 3: 
 4: 
 5: 
[<Literal>]
let PRECISION = 1e-6

let inline shouldEqual msg exp act =
    Assert.AreEqual(exp, float act, PRECISION, msg)

The only thing left was to implement the tests itself.

If you’re choosing the framework too, I’d recommend to check the thread Question about Testing Frameworks in the F# OpenSource list.

The Search for Values

And this is the point I stuck at: there’re almost 200,000 of tests! If the values match Excel results, the test is passed… but we don’t want to use Excel interop because of mono. Ok, let’s just store the arguments and expected values. There’re 52 functions to test, in general case the different functions have different number and/or types of parameters. What is the best way to read them from the files, parse, pass to tests, compare the actual and expected results?

Here’s the most concise solution I came up with, it’s based on F#’s inlining (source is here):

 1: 
 2: 
 3: 
 4: 
 5: 
[<Test>]  
let disc() = runTests "disc" parse6 Financial.Disc   

[<Test>]
let price() = runTests "price" parse8 Financial.Price

The test data is read from the “price” file. The function has 8 arguments (luckily, it’s a tuple for C# users convenience) and the compiler is able to infer all the types - in our case they’re standard (float, DateTime, int…) and have TryParse static method.

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
let inline parse str =
    let mutable res = Unchecked.defaultof<_>
    let _ = (^a: (static member TryParse: string * byref< ^a > -> bool) (str, &res))
    res

let inline parse3 [| a; b; c |] =
    (parse a, parse b), parse c 

As you may notice, there’re some duplications - we still need several parseN functions.

Back to our minimalistic approach, the function runTests should be inlined too, so the compiler can infer the types. The error messages are pretty detailed: they contain index, function name, arguments, expected and actual results - everything you hopefully won’t ever need to know :)

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
let inline runTests fname parsef f =
    readTestData fname    
    |> Seq.iteri (fun i data ->
        let param, expected = parsef data
        let actual = f param
        shouldEqual (sprintf "%d - %s(%A)" i fname param) expected actual)

Original console tests are in the repo too - all tests pass (checked for Excel 2010).

Testing the Properties

Something new this version got is the tests of some properties, e.g. the modified duration shouldn’t be greater than maturity, the cumulative accrint should be equal to the sum of accrints in corresponding coupon periods etc.

The best way to do that is through randomized testing with FsCheck.

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
[<Test>]
let ``tbill price is less than 100``() =
    fsCheck (fun (sd: DateTime) t disc' ->
        let md, disc = sd.AddDays (toFloat t * 365.), toFloat disc'

        tryTBillPrice sd md disc    
        ==>
        lazy (Financial.TBillPrice(sd, md, disc) - 100. < PRECISION))

The most challenging part is expressing the properties in the way FsCheck can understand easily. For example, each function has tryX method to check the parameters. That’s exactly what we need to be sure the random values make sense. However, when the logic behind such a function is not trivial, the arguments just get exhausted.

First, the dates - when there’re several dates and you know which one comes first, how many days can be between them etc, it makes sense to keep as date the first one only and define the others as date + n days (as in the example above).

Second, converting ints to floats usually works better than straightforward generating of floats. Probably, there’re other (and better) ways - for example, I was thinking about writing a smart generator both for dates and floats - prices, rates and so on, but the current approach worked quite well and solved the problem with arguments.

Several community projects using FsCheck:

Bonus Reading:

Why Does Anyone Care?

This is the only library replicating Excel behavior, with a few differences, though, as explained here. Anyone who needs this functionality can use the library on any platform.

What was also interesting is to compare the results with other apps. For example, spotted case when Excel for Mac didn’t match Windows Excel results, that was a bit surprising. But in Libre Office almost every function is more or less different!

Some OO functions just returned errors for any input…

Seems like they can be fixed now ;)

 

 

It almost looks like a property

Based on a true story. Only the names, places, and events have been changed.

Character

Many people today forget that the tools are just tools. Doesn’t matter if it’s OO or functional programming, take any language: there’re plenty of ways to write awesome and even more - to write terrifying code. lets . compose . these . functions . until . your . type-checker
. explodes . btw . did . you . save . everything?

[Originally posted here.]

Today I’ll share a simple code-reducing trick, which I came up with a couple of years ago, but still like it for some reason (if you don’t - just deal with it). So what the problem is?

Now answer the question - Do you believe there might be anything except the classes?
(No) No problem. At all. There’s a class, with a bunch of fields. The end.
(Yes) Ok, let’s say the objects of this class suppose to describe the weather conditions. Moving on,

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
class Weather
{
    public double Temp { get; set; }
    public double Humidity { get; set; }
    public double Pressure { get; set; } 
    public double Rainmm { get; set; }
}

You may say, the temperature is usually reported as a range. Here we go:

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
class Weather
{
    public double MinTemp { get; set; }
    public double MaxTemp { get; set; }
    public double WindChill { get; set; } 
    ... // other properties
}

Nice, but there’re also different observations for the night/morning/afternoon/evening forecasts, which one may want to include into the same forecast:

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
class Weather
{
    public double MinNightTemp { get; set; }
    public double MinDayTemp { get; set; }
    ... // other properties
}

Now add several predictions for future dates, whatever. On a rainy day you may find out the number of fields somehow approached fifty, *sigh*, add another one and forget about them.

And then the Universe decides that it was not fun enough and you should make the predictions more consistent with reality (or at least pretend to do so), for example, by including the information from other sources like <insert your favorite weather forecast website here, I’d better take an umbrella anyway>. You quickly run an experiment with the min/max temperatures and expect to get a 10% improvement in accuracy by weighting the observations from different sources… But there’re more than 50 properties already, remember?

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
var w1 = 0.8;
var w2 = 1 - w1;
var result = new Weather 
{ 
    MinNightTemp = w1 * orig.MinNightTemp + w2 * ext.MinNightTemp, 
    MaxNightTemp = w1 * orig.MaxNightTemp + w2 * ext.MaxNightTemp,  
    ... // oh, why is that intern on vacation now? who will fill in all this stuff?
}

It could be anything else, the key is that there’s a set of operations, similar for all the fields, and the fields are somewhat similar too - in our case, they’re all of type double.

Let’s create an array instead, where the fields are the elements of this array:

 1: 
 2: 
1: for(int i = 0; i < n; ++i)
2:     x[i] = w1 * y[i] + w2 * z[i];

But… but now we don’t know what is what. Ok, here’s the answer:

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
enum P 
{
    MinNightTemp,
    MaxNightTemp,
    ...
}
...   
var res = x[(int)P.MinNightTemp];

What about a kind of vector? It looks even better:

 1: 
 2: 
var x = w1 * y + w2 * z; var res = x[P.MinNightTemp]; // compare to x.MinNightTemp
  • Saves a lot of typing, decreses the probability of potential mistakes;
  • You don’t need even a dictionary - an array is already enough;
  • The enum part can help to preserve the order (in original case that was mandatory requirement), nothing breaks if you add a new feature somewhere in the middle;
  • Every field automatically has its own underlying value, you avoid the mistake of assigning the same value to the different constants (e.g. MinTempInd = 42 and MaxTempInd = 42);
  • There’re good chances no one will ever notice, that Weather is not a separate class any more. It looks almost like a property.

 

Instead of conclusion

There’re simple yet annoying problems, and to solve them you don’t need a dataframe or any complicated data structure or even a separate class. C# is not that bad in the end, no one forces you to write only overly-OOP code.

That doesn’t mean you can’t come up with more concise solutions using a different language or solving a different problem. Only that there’re ways - in everything - to do the familiar things differently, the ways not everyone notices.

 

Log(music) with Undertone

Whether you’re a physicist, biologist, economist, data scientist, whatever, the chances are you’ll meet the Markov chains sooner or later. And today we’ll try them on… music!

Christmas

If you’re not up to the math part, feel free to skip it.

[Originally posted here].

Math

The system states changes are represented with a help of transition matrix, describing the probabilities of possible transitions. This matrix is a key to determine the state of the model after n steps.

The math behind that is very straightforward - we need just a matrix power. Oh, wait - what if it is fractional, for example, when representing time? There’re many ways to compute it, let’s take a look at some of them.

1. Eigen vectors

This method is based on properties of eigen values decomposition:

That’s actually how Matlab mpower function looks like:

 1: [V,D] = eig(M)
 2: An = V * (D .^ n) * inv(V) 

We can take advantage of Math.NET matrix decompositions to implement this method:

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
open namespaces, define constants

/// Matrix power: eigen values decomposition method
/// M^p = V * (D .^ p) / V
let mpowEig (m: SparseMatrix) p =
    let evd = m.Evd()
    let v, d = evd.EigenVectors(), evd.D()
    mapInPlace (fun x -> x ** p) d
    v * d * v.Inverse()

What are the pros/cons?

+ the method is rather efficient
- inverse matrix might not exist 1

- some matrices cannot be diagonalised:

2. Generator matrices

If we can find a matrix , such as , then (remember series?)

We are interested in finding a generator with zero row sums and nonnegative off-diagonal entries. So we’ll use Mercator’s series:

In matrix case it’s

The common concern is that the series methods are not very efficient and may be not accurate enough. It’s easy to compare the methods by computing the squared error between the actual and expected values (say, compute ).

Another problem is that a generator matrix (and power too) might have negative elements - not really what you want when dealing with probabilities. We can make them non-negative and keep the rows sums equal to zero… but the power properties won’t be preserved. For example, you can find out that .

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
17: 
18: 
19: 
20: 
21: 
22: 
23: 
24: 
25: 
26: 
27: 
28: 
29: 
30: 
31: 
32: 
33: 
34: 
35: 
/// generator matrix
let genm (a: SparseMatrix) =
    let size = a.RowCount
    let mid = constDiag size 1.0
    let a_i = a - mid

    let rec logm (i, qpow: SparseMatrix) res =
        if i = MAX_ITER then res
        else
           let qdiv = qpow.Divide (float i)
           if sumSq qdiv < PREC then res + qdiv
           else 
               logm (i+1, -qpow * a_i) (res + qdiv)
    // find log(A - I)
    // update negative off-diagonal entries, the row-sum = 0
    logm (1, a_i) (zeroCreate size size) |> updateneg

/// exp series approximation
let expSeries (m: SparseMatrix)  =
    let size = m.RowCount
    
    let rec sum (i, k, mpow: SparseMatrix) res =
        if i = MAX_ITER then res
        else
            let mpowi, ki = mpow * m, k * float i
            let mdiv = mpowi.Divide ki
            if sumSq mdiv < PREC then res + mdiv
            else
                sum (i+1, ki, mpowi) (res + mdiv) 

    let mid = constDiag size 1.0
    sum (1, 1.0, mid) mid

/// Matrix power: series method
let mpowSeries m (p: float) = SparseMatrix.OfMatrix (genm m * p) |> expSeries

+ Works for non-diagonizable matrices
- Not very efficient
- May be not accurate enough or even fail to converge

These are just the most simple algorithms - would be interesting to check the others, but let’s move on.

References

  1. 19 Dubious Ways to Compute the Exponential of a Matrix
  2. Finding Generators for Markov Chains via empirical transition matrices, with applications to credit ratings

Music

It’s not a secret, that Markov models are often used for generative music, here’s one of the examples: Programming Instrumental Music From Scratch (there’re also some interesting references in the comments). We’ll do a different thing, though: given notes, we can calculate the transition probabilities and… generate new melodies.

But everything in its time.

We need to somehow define the sequences of notes and play them. And thanks to Undertone that’s pretty easy! First of all, let’s download the piano samples - University of Iowa has a collection of recordings here. You can get them manually or using this script. It’s also ok to generate the notes programmatically - but somewhat a bit more realistic is more fun ^_^

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
17: 
18: 
19: 
20: 
21: 
22: 
23: 
24: 
25: 
26: 
27: 
28: 
29: 
30: 
31: 
32: 
33: 
34: 
35: 
36: 
37: 
38: 
39: 
40: 
41: 
42: 
43: 
44: 
45: 
46: 
47: 
48: 
49: 
50: 
open namespaces

// Size: 2/4

/// 1/16 (semiquaver)
[<Literal>] 
let Sq = 0.0625
/// 1/8 (quaver)
[<Literal>] 
let Q1 = 0.125
/// 3/8
[<Literal>] 
let Q3 = 0.375

...

/// Mapping between note enum and file name convention 
let noteToPianoName note =
    match note with
    | Note.C         -> "C" 
    | Note.Csharp    -> "Db"
    | Note.D         -> "D"
    | Note.Dsharp    -> "Eb"
    | Note.E         -> "E"
    | Note.F         -> "F"
    | Note.Fsharp    -> "Gb"
    | Note.G         -> "G"
    | Note.Gsharp    -> "Ab"
    | Note.A         -> "A"
    | Note.Asharp    -> "Bb"
    | Note.B         -> "B"
    | _ -> failwith "invalid note"

/// Read piano note from .aiff file
let readPianoNote() =
    let cache = Dictionary<_,_> HashIdentity.Structural
    let inline path (noteName, octave) = 
        Path.Combine(dir, "Piano.ff." + noteName + string octave + ".aiff")

    fun (note, octave) ->
        let noteKey = noteToPianoName note, octave
        match cache.TryGetValue noteKey with
        | true, wave -> wave
        | _ -> 
            let wave = IO.read (path noteKey) |> Seq.toArray
            cache.Add(noteKey, wave)
            wave

let makeNote = readPianoNote() 
let getMelody = Seq.collect (fun (noct, time) -> makeNote noct |> Seq.take (noteValue time))

The note is defined with its symbol (C = Do, F# = Fis = Fa-diesis = F sharp), octave and time (♪ - 1/8, ♬ - 2 notes * 1/16). The sequence of notes gives a melody, which can be played and even saved later. Can you guess the following one?

Original

Sounds recognizable, but awkwardly stumbling - makes me feel much better about my piano skills. I’d recommend to watch this performace, it’s amazing! (and that pianist-feeling suddenly goes down again)). I believe there’s a way to make the generated sound smoother and so on, but that’s not the goal now.

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
Define commonly used notes

/// Sample melody
let ent = [
    D3,Sq; Ds3,Sq
    E3,Sq; C4,Q1; E3,Sq; C4,Q1; E3,Sq; C4,Q3
    C4,Sq; D4,Sq; Ds4,Sq
    E4,Sq; C4,Sq; D4,Sq; E4,Q1; C4,Sq; D4,Q1
    C4,Sq; D3,Sq; Ds3,Sq ...
    C4,Q3
]

Player.Play (getMelody ent)

(Note, the example above has 3 times faster tempo that the one you’ll get with downloaded notes)

The next thing to do is to calculate the transition probabilities. Consider the following sequence: E3 - 1/16, C4 - 1/8, E3 - 1/16, C4 - 3/8, C4 - 1/16.

The first option is to simply count transitions between notes and then normalize the rows, so the sum are equal to 1:


But wait, what about time? There’re different notes, we can take that into account too: if the ‘smallest’ length is 1/16, then we can count, how many such intervals each note takes, multiplying by 162:

E3 - 1/16, C4 - 1/8, E3 - 1/16, C4 - 3/8, C4 - 1/16

So, we know the transitions and times when the next note should be played. A uniform random number generator will choose this next note: e.g. if transitions are 0.4 0.2 0.1 0.3 and generated number is 0.65, that’s the 3rd one (this operation is called choose below).
The note at time t, given transition matrix P is

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
17: 
18: 
19: 
20: 
21: 
22: 
23: 
24: 
25: 
26: 
27: 
28: 
29: 
30: 
31: 
32: 
33: 
34: 
35: 
36: 
37: 
38: 
39: 
40: 
41: 
42: 
43: 
44: 
45: 
46: 
47: 
/// Compute the times and probabilities of transitions between notes
let transitions xs =
    let noteToId = Seq.distinctBy fst xs |> Seq.mapi (fun i (note, _) -> note, i) |> dict

    // # of unique notes and # of notes in melody
    let n, m = noteToId.Count, List.length xs

    let matrix = Array.init n (fun _ -> Array.zeroCreate n)
    let times = Array.create m 0.

    // fill the matrix
    Seq.pairwise xs
    |> Seq.iteri (fun i ((note1, time1), (note2, _)) ->
        let n1, n2 = noteToId.[note1], noteToId.[note2]
        // option #1 - ignore times
        //   matrix.[n1].[n2] <- matrix.[n1].[n2] + 1.
        // option #2 - add extra probability for transition to the same note
        //    matrix.[n1].[n1] <- matrix.[n1].[n1] + time1 * 16.0
        //    matrix.[n1].[n2] <- matrix.[n1].[n2] + 1.0
        // option #3 - take time into account (scale by 16)
        let p = time1*16.-1. in if p > 0. then matrix.[n1].[n1] <- matrix.[n1].[n1] + p
        matrix.[n1].[n2] <- matrix.[n1].[n2] + 1.0
        
        times.[i+1] <- times.[i] + time1)
    
    normalizeRows matrix, times, Seq.toArray noteToId.Keys

let matrix, times, idToNote = transitions ent

/// Find index of the next note given probability
let transTo (ps: _[]) = ...

/// Generate music from transition matrices in a file
let genMusic matricesFileName (fstNote, idToNote: _[]) =
    // transition matrices and times
    let ps, ts = readMatrices matricesFileName
    let n = ps.Count
    let rand = System.Random()
    
    let rec gen i prevInd res =
        if i >= n-1 then List.rev res
        else
            let j = transTo ps.[i].[prevInd] (rand.NextDouble())
            let next = idToNote.[j], ts.[i]
            gen (i+1) j (next::res)

    gen 0 0 [fstNote] 

I tried to compute the transition matrix - and both method failed! Turns out this is the real-life ‘unlucky’ matrix. To be sure that doesn’t work, checked R’s expm package:

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
logm(x, method = "Eigen")
Error in logm(x, method = "Eigen") : non diagonalisable matrix
logm(x, method = "Higham08")
Error in solve.default(X[ii, ii] + X[ij, ij], S[ii, ij] - sumU) :
system is computationally singular: reciprocal condition number = 0
In addition: Warning messages:
1: In sqrt(S[ij, ij]) : NaNs produced
2: In sqrt(S[ij, ij]) : NaNs produced

Fortunately, we don’t need fraction power here: if we count time in 1/16th, it becomes integer! 1/16 is 1, 3/8 is 6 and so on. Good old multiplication saves the situation.

Here’s several examples of generated melodies - there’re only 10 notes are in the game, but the results are already pretty different from the original:

Generated-1
Generated-2
Generated-3

What if we add an ‘extra-probability’ for note to stay in the same state?

E3 - 1/16, C4 - 1/8, E3 - 1/16, C4 - 3/8, C4 - 1/16

For this matrix, the generator does exist and we can programmatically compose something like that:

Generated-series-1
Generated-series-2

Why stop here? We could define the similar transitions for note length; or add chords - say, generated according to the current tonality; or switch tonalities/modes… Well, there’s a whole ‘musical math’, solfeggio, full of different functions - plenty of space to experiment!

Happy Holidays!

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
let jingleBells = [
    E4,Sq; E4,Sq; E4,Q1
    E4,Sq; E4,Sq; E4,Q1
    E4,Sq; G4,Sq; C4,Sq; D4,Sq
    E4,Q1*2.0
    F4,Sq; F4,Sq; F4,Sq; F4,Sq
    F4,Sq; E4,Sq; E4,Sq; E4,Sq
    E4,Sq; D4,Sq; D4,Sq; E4,Sq
    D4,Q1; G4,Q1
]

Player.Play (getMelody jingleBells)

     

open MathNet.Numerics.LinearAlgebra.Double
open MathNet.Numerics.LinearAlgebra.Generic
open Matrix
open SparseMatrix

[<Literal>]
let MAX_ITER = 10000
[<Literal>]
let PREC = 1e-60

val mpowEig : m:SparseMatrix -> p:float -> Matrix<float>

Full name: LogMusic.mpowEig

 Matrix power: eigen values decomposition method
 M^p = V * (D .^ p) / V

val m : SparseMatrix
Multiple items
type SparseMatrix =
  inherit Matrix
  new : storage:SparseCompressedRowMatrixStorage<float> -> SparseMatrix + 6 overloads
  member CreateMatrix : numberOfRows:int * numberOfColumns:int * ?fullyMutable:bool -> Matrix<float>
  member CreateVector : size:int * ?fullyMutable:bool -> Vector<float>
  member FrobeniusNorm : unit -> float
  member IndexedEnumerator : unit -> IEnumerable<Tuple<int, int, float>>
  member InfinityNorm : unit -> float
  member IsSymmetric : bool
  member KroneckerProduct : other:Matrix<float> * result:Matrix<float> -> unit
  member LowerTriangle : unit -> Matrix<float> + 1 overload
  member NonZerosCount : int
  ...

Full name: MathNet.Numerics.LinearAlgebra.Double.SparseMatrix

--------------------
SparseMatrix(storage: MathNet.Numerics.LinearAlgebra.Storage.SparseCompressedRowMatrixStorage<float>) : unit
SparseMatrix(order: int) : unit
SparseMatrix(rows: int, columns: int) : unit

val p : float
val evd : Factorization.Evd
Matrix.Evd() : Factorization.Evd
val v : Matrix<float>
val d : Matrix<float>
Factorization.Evd.EigenVectors() : Matrix<float>
Factorization.Evd.D() : Matrix<float>
val mapInPlace : f:(float -> float) -> A:#Matrix<float> -> unit

Full name: MathNet.Numerics.LinearAlgebra.Double.Matrix.mapInPlace

val x : float
Matrix.Inverse() : Matrix<float>
val genm : a:SparseMatrix -> Matrix<float>

Full name: LogMusic.genm

 generator matrix

val a : SparseMatrix
val size : int
property Matrix.RowCount: int
val mid : SparseMatrix
val constDiag : n:int -> f:float -> SparseMatrix

Full name: MathNet.Numerics.LinearAlgebra.Double.SparseMatrix.constDiag

val a_i : SparseMatrix
val logm : (int * SparseMatrix -> Matrix<float> -> Matrix<float>)
val i : int
val qpow : SparseMatrix
val res : Matrix<float>
val MAX_ITER : int

Full name: LogMusic.MAX_ITER

val qdiv : Matrix<float>
Matrix.Divide(scalar: float) : Matrix<float>
Matrix.Divide(scalar: float, result: Matrix<float>) : unit
Multiple items
val float : value:'T -> float (requires member op_Explicit)

Full name: Microsoft.FSharp.Core.Operators.float

--------------------
type float = System.Double

Full name: Microsoft.FSharp.Core.float

--------------------
type float<'Measure> = float

Full name: Microsoft.FSharp.Core.float<_>

val sumSq : m:Matrix<float> -> float

Full name: LogMusic.sumSq

 sum of squares

val PREC : float

Full name: LogMusic.PREC

val zeroCreate : rows:int -> cols:int -> SparseMatrix

Full name: MathNet.Numerics.LinearAlgebra.Double.SparseMatrix.zeroCreate

val updateneg : m:Matrix<float> -> Matrix<float>

Full name: LogMusic.updateneg

 update negative

val expSeries : m:SparseMatrix -> Matrix<float>

Full name: LogMusic.expSeries

 exp series approximation

val sum : (int * float * SparseMatrix -> Matrix<float> -> Matrix<float>)
val k : float
val mpow : SparseMatrix
val mpowi : SparseMatrix
val ki : float
val mdiv : Matrix<float>
val mpowSeries : m:SparseMatrix -> p:float -> Matrix<float>

Full name: LogMusic.mpowSeries

 Matrix power: series method

SparseMatrix.OfMatrix(matrix: Matrix<float>) : SparseMatrix
open System.IO
open System.Collections.Generic
open Undertone
open Undertone.Waves
Multiple items
type LiteralAttribute =
  inherit Attribute
  new : unit -> LiteralAttribute

Full name: Microsoft.FSharp.Core.LiteralAttribute

--------------------
new : unit -> LiteralAttribute

val Sq : float

Full name: LogMusic.Sq

 1/16 (semiquaver)

val Q1 : float

Full name: LogMusic.Q1

 1/8 (quaver)

val Q3 : float

Full name: LogMusic.Q3

 3/8

let inline noteValue time = Time.noteValue 15. time

/// Note samples directory
[<Literal>]
let dir = @"C:\samples"

val noteToPianoName : note:Note -> string

Full name: LogMusic.noteToPianoName

 Mapping between note enum and file name convention

val note : Note
type Note =
  | Cflat = -1
  | C = 0
  | Csharp = 1
  | Dflat = 1
  | D = 2
  | Dsharp = 3
  | Eflat = 3
  | E = 4
  | Fflat = 4
  | Esharp = 5
  | F = 5
  | Fsharp = 6
  | Gflat = 6
  | G = 7
  | Gsharp = 8
  | Aflat = 8
  | A = 9
  | Asharp = 10
  | Bflat = 10
  | B = 11
  | Bsharp = 12

Full name: Undertone.Note

Note.C: Note = 0
Note.Csharp: Note = 1
Note.D: Note = 2
Note.Dsharp: Note = 3
Note.E: Note = 4
Note.F: Note = 5
Note.Fsharp: Note = 6
Note.G: Note = 7
Note.Gsharp: Note = 8
Note.A: Note = 9
Note.Asharp: Note = 10
Note.B: Note = 11
val failwith : message:string -> 'T

Full name: Microsoft.FSharp.Core.Operators.failwith

val readPianoNote : unit -> (Note * int -> float [])

Full name: LogMusic.readPianoNote

 Read piano note from .aiff file

val cache : Dictionary<(string * int),float []>
Multiple items
type Dictionary<'TKey,'TValue> =
  new : unit -> Dictionary<'TKey, 'TValue> + 5 overloads
  member Add : key:'TKey * value:'TValue -> unit
  member Clear : unit -> unit
  member Comparer : IEqualityComparer<'TKey>
  member ContainsKey : key:'TKey -> bool
  member ContainsValue : value:'TValue -> bool
  member Count : int
  member GetEnumerator : unit -> Enumerator<'TKey, 'TValue>
  member GetObjectData : info:SerializationInfo * context:StreamingContext -> unit
  member Item : 'TKey -> 'TValue with get, set
  ...
  nested type Enumerator
  nested type KeyCollection
  nested type ValueCollection

Full name: System.Collections.Generic.Dictionary<_,_>

--------------------
Dictionary() : unit
Dictionary(capacity: int) : unit
Dictionary(comparer: IEqualityComparer<'TKey>) : unit
Dictionary(dictionary: IDictionary<'TKey,'TValue>) : unit
Dictionary(capacity: int, comparer: IEqualityComparer<'TKey>) : unit
Dictionary(dictionary: IDictionary<'TKey,'TValue>, comparer: IEqualityComparer<'TKey>) : unit

module HashIdentity

from Microsoft.FSharp.Collections

val Structural<'T (requires equality)> : IEqualityComparer<'T> (requires equality)

Full name: Microsoft.FSharp.Collections.HashIdentity.Structural

val path : (string * 'a -> string)
val noteName : string
val octave : 'a
type Path =
  static val DirectorySeparatorChar : char
  static val AltDirectorySeparatorChar : char
  static val VolumeSeparatorChar : char
  static val InvalidPathChars : char[]
  static val PathSeparator : char
  static member ChangeExtension : path:string * extension:string -> string
  static member Combine : params paths:string[] -> string + 3 overloads
  static member GetDirectoryName : path:string -> string
  static member GetExtension : path:string -> string
  static member GetFileName : path:string -> string
  ...

Full name: System.IO.Path

Path.Combine(params paths: string []) : string
Path.Combine(path1: string, path2: string) : string
Path.Combine(path1: string, path2: string, path3: string) : string
Path.Combine(path1: string, path2: string, path3: string, path4: string) : string
val dir : string

Full name: LogMusic.dir

 Note samples directory

Multiple items
val string : value:'T -> string

Full name: Microsoft.FSharp.Core.Operators.string

--------------------
type string = System.String

Full name: Microsoft.FSharp.Core.string

val octave : int
val noteKey : string * int
Dictionary.TryGetValue(key: string * int, value: byref<float []>) : bool
val wave : float []
module IO

from Undertone

val read : file:string -> seq<float>

Full name: Undertone.IO.read

module Seq

from Microsoft.FSharp.Collections

val toArray : source:seq<'T> -> 'T []

Full name: Microsoft.FSharp.Collections.Seq.toArray

Dictionary.Add(key: string * int, value: float []) : unit
val makeNote : (Note * int -> float [])

Full name: LogMusic.makeNote

val getMelody : (((Note * int) * float) list -> seq<float>)

Full name: LogMusic.getMelody

val collect : mapping:('T -> #seq<'U>) -> source:seq<'T> -> seq<'U>

Full name: Microsoft.FSharp.Collections.Seq.collect

val noct : Note * int
val time : float
val take : count:int -> source:seq<'T> -> seq<'T>

Full name: Microsoft.FSharp.Collections.Seq.take

val noteValue : time:float -> int

Full name: LogMusic.noteValue

let D3 = Note.D, 3
let Ds3 = Note.Dsharp, 3
let E3 = Note.E, 3
let Fs3 = Note.Fsharp, 3
let G3 = Note.G, 3
let A3 = Note.A, 3
let C4 = Note.C, 4
let D4 = Note.D, 4
let Ds4 = Note.Dsharp, 4
let E4 = Note.E, 4

let G4 = Note.G, 4
let F4 = Note.F, 4

val ent : ((Note * int) * float) list

Full name: LogMusic.ent

 Sample melody

val D3 : Note * int

Full name: LogMusic.D3

val Ds3 : Note * int

Full name: LogMusic.Ds3

val E3 : Note * int

Full name: LogMusic.E3

val C4 : Note * int

Full name: LogMusic.C4

val D4 : Note * int

Full name: LogMusic.D4

val Ds4 : Note * int

Full name: LogMusic.Ds4

val E4 : Note * int

Full name: LogMusic.E4

E3,Sq; C4,Q1; E3,Sq; C4,Q1; E3,Sq; C4,Q3
    A3,Sq; G3,Sq
    Fs3,Sq; A3,Sq; C4,Sq; E4,Q1; D4,Sq; C4,Sq; A3,Sq
    D4,Q3; D3,Sq; Ds3,Sq
    E3,Sq; C4,Q1; E3,Sq; C4,Q1; E3,Sq; C4,Q3
    C4,Sq; D4,Sq; Ds4,Sq
    E4,Sq; C4,Sq; D4,Sq; E4,Q1; C4,Sq; D4,Q1
    C4,Q3; D4,Sq; Ds4,Sq
    E4,Sq; C4,Sq; D4,Sq; E4,Q1; C4,Sq; D4,Sq; C4,Sq
    E4,Sq; C4,Sq; D4,Sq; E4,Q1; C4,Sq; D4,Sq; C4,Sq
    E4,Sq; C4,Sq; D4,Sq; E4,Q1; C4,Sq; D4,Q1
type Player =
  private new : unit -> Player
  static member Play : sampleSource:seq<float> -> IPlayer

Full name: Undertone.Player

static member Player.Play : sampleSource:seq<float> -> IPlayer
val transitions : xs:('a * float) list -> float [] [] * float [] * 'a [] (requires equality)

Full name: LogMusic.transitions

 Compute the times and probabilities of transitions between notes

val xs : ('a * float) list (requires equality)
val noteToId : IDictionary<'a,int> (requires equality)
val distinctBy : projection:('T -> 'Key) -> source:seq<'T> -> seq<'T> (requires equality)

Full name: Microsoft.FSharp.Collections.Seq.distinctBy

val fst : tuple:('T1 * 'T2) -> 'T1

Full name: Microsoft.FSharp.Core.Operators.fst

val mapi : mapping:(int -> 'T -> 'U) -> source:seq<'T> -> seq<'U>

Full name: Microsoft.FSharp.Collections.Seq.mapi

val note : 'a (requires equality)
val dict : keyValuePairs:seq<'Key * 'Value> -> IDictionary<'Key,'Value> (requires equality)

Full name: Microsoft.FSharp.Core.ExtraTopLevelOperators.dict

val n : int
val m : int
property ICollection.Count: int
Multiple items
type List<'T> =
  new : unit -> List<'T> + 2 overloads
  member Add : item:'T -> unit
  member AddRange : collection:IEnumerable<'T> -> unit
  member AsReadOnly : unit -> ReadOnlyCollection<'T>
  member BinarySearch : item:'T -> int + 2 overloads
  member Capacity : int with get, set
  member Clear : unit -> unit
  member Contains : item:'T -> bool
  member ConvertAll<'TOutput> : converter:Converter<'T, 'TOutput> -> List<'TOutput>
  member CopyTo : array:'T[] -> unit + 2 overloads
  ...
  nested type Enumerator

Full name: System.Collections.Generic.List<_>

--------------------
List() : unit
List(capacity: int) : unit
List(collection: IEnumerable<'T>) : unit

val length : list:'T list -> int

Full name: Microsoft.FSharp.Collections.List.length

val matrix : float [] []
module Array

from Microsoft.FSharp.Collections

val init : count:int -> initializer:(int -> 'T) -> 'T []

Full name: Microsoft.FSharp.Collections.Array.init

val zeroCreate : count:int -> 'T []

Full name: Microsoft.FSharp.Collections.Array.zeroCreate

val times : float []
val create : count:int -> value:'T -> 'T []

Full name: Microsoft.FSharp.Collections.Array.create

val pairwise : source:seq<'T> -> seq<'T * 'T>

Full name: Microsoft.FSharp.Collections.Seq.pairwise

val iteri : action:(int -> 'T -> unit) -> source:seq<'T> -> unit

Full name: Microsoft.FSharp.Collections.Seq.iteri

val note1 : 'a (requires equality)
val time1 : float
val note2 : 'a (requires equality)
val n1 : int
val n2 : int
val normalizeRows : m:float [] [] -> float [] []

Full name: LogMusic.normalizeRows

 Normalize rows, so the row-sums are equal to 1

property IDictionary.Keys: ICollection<'a>
val matrix : float [] []

Full name: LogMusic.matrix

val times : float []

Full name: LogMusic.times

val idToNote : (Note * int) []

Full name: LogMusic.idToNote

val transTo : ps:float [] -> (float -> int)

Full name: LogMusic.transTo

 Find index of the next note given probability

val ps : float []
let last = ps.Length-1
    let rec findj j v =
        if ps.[j] >= v || j = last then j
        else findj (j+1) (v - ps.[j])
    findj 0
val genMusic : matricesFileName:string -> fstNote:('a * float) * idToNote:'a [] -> ('a * float) list

Full name: LogMusic.genMusic

 Generate music from transition matrices in a file

val matricesFileName : string
val fstNote : 'a * float
val idToNote : 'a []
val ps : List<ResizeArray<float []>>
val ts : float []
val readMatrices : path:string -> List<ResizeArray<float []>> * float []

Full name: LogMusic.readMatrices

 Read transition matrices and times

property List.Count: int
val rand : System.Random
namespace System
Multiple items
type Random =
  new : unit -> Random + 1 overload
  member Next : unit -> int + 2 overloads
  member NextBytes : buffer:byte[] -> unit
  member NextDouble : unit -> float

Full name: System.Random

--------------------
System.Random() : unit
System.Random(Seed: int) : unit

val gen : (int -> int -> ('a * float) list -> ('a * float) list)
val prevInd : int
val res : ('a * float) list
val rev : list:'T list -> 'T list

Full name: Microsoft.FSharp.Collections.List.rev

val j : int
System.Random.NextDouble() : float
val next : 'a * float
val jingleBells : ((Note * int) * float) list

Full name: LogMusic.jingleBells

val G4 : Note * int

Full name: LogMusic.G4

val F4 : Note * int

Full name: LogMusic.F4

  1. the example matrix isn’t actually a transition matrix - the row sums should be equal to 1 - it’s here just to demonstrate a possible problem.

  2. the sequence C4 - 1/8, E3 - 1/16 is almost the same as C4 - 1/16, C4 - 1/16, E3 - 1/16. That gives us 50/50 transitions: after each period of time C4 can either stay as C4 or transform to E3.

Trying out Deedle with Bones and Regression

I usually don’t need to run a regression anywhere, but it’s kind of chasing me recently, starting with the Asset Pricing class and several variations of returns regressions (signed up to look at the familiar things from a different point of view… well, I definitely succeeded: have you ever thought about drawing the returns, prices and discount factors in space, all at once? 1. But I ‘cheated’ and completed the assignments with R.

[Originally posted here.]

Though that was only the beginning - my cousin, MD student, was measuring the deflection of bones and other samples with different loads. And this time I decided to try out Deedle and help her to explore the data.

Hint for Mono users: if XS doesn’t load the main Deedle project, you can manually update the fsproj file (delete the reference to FSharp.Core), reload it, add references to Math.NET and FSharp.Data libs - it’ll work nicely.

Let’s start with loading experimental data. No more string splits and manual parsing!

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
#r "FSharp.Data.dll"
#load "Deedle.fsx"

open System
open Deedle
// load data from a csv file
let ds = Frame.ReadCsv(__SOURCE_DIRECTORY__ + "/deflection.csv")
val ds : Frame<int,string> = 
        p         sample type length<mm> width<mm> height<mm> deflection<mm> 
  0  -> 0.98      bone        47         7.3       4.4        0.06 
  1  -> 1.96      bone        47         7.3       4.4        0.12 

Then we checked out some properties of different samples groups - say, the average sample deformation. For simplicity we’ll use only the bones group, load (“p”) and deflection columns.

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16:   
17: 
18: 
19: 
20: 
21: 
22: 
23: 
24: 
25: 
26: 
27: 
// choose several columns and group the data by sample type
let bySample = 
    ds.Columns.[["sample type"; "p"; "deflection<mm>"]] 
    |> Frame.groupRowsByString "sample type"
val bySample : Frame<(string * int),string> =
                  sample type p         deflection<mm> 
  bone      0  -> bone        0.98      0.06           
            1  -> bone        1.96      0.12           
  ...      ...    ...         ...                      
  duralumin 6  -> duralumin   0.98      0.03           
  ...      ...    ...         ...                      

// average deflection by sample type
bySample.Columns.[["sample type"; "deflection<mm>"]] |> Frame.meanLevel Pair.get1Of2
val it : Frame<string,string> =
               deflection<mm>    
  bone      -> 0.228333333333333 
  duralumin -> 0.116666666666667 
  ...       -> ...               

// select the data for bones
let bones = (Frame.nest bySample).["bone"]
val bones : Frame<int,string> =
       sample type p         deflection<mm> 
 0  -> bone        0.98      0.06           
 1  -> bone        1.96      0.12           
 ...-> ...         ...       ...            

You may notice that some of the values in the table are missing (the handwriting can be completely unparsable!), by default they are omited, but we can always specify how we want this data to be filled using Direction or a custom function.

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
// note that there're missing values in this dataset
let deflections = bones?``deflection<mm>``
val deflections : Series<int,float> =
 0  -> 0.06      
 1  -> 0.12      
 ...-> ...       
 12 -> <missing> 
Series.mean deflections
 val it : float = 0.2283333333 
// omit missing values
deflections |> Series.dropMissing |> Series.mean
 val it : float = 0.2283333333 
// fill missing values by copying forward
deflections |> Series.fillMissing Direction.Forward |> Series.mean
 val it : float = 0.2542857143 

Now let’s check if there’s any relation between the deflection and load. In theory, it’s supposed to be linear and we’re going to test that with a linear regression.

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
/// Find slope and intercept with linear regression
let linearRegression xs ys = (...)
// drop rows with missing values
let bonesreg = Frame.dropSparseRows bones
val bonesreg : Frame<int,string> =
        sample type p    deflection<mm> 
   0 -> bone        0.98 0.06           
 ... -> ...         ...  ...            
   5 -> bone        5.88 0.41           

let load = Series.values bonesreg?p
let defl = Series.values bonesreg?``deflection<mm>``
let slope, intercept = linearRegression load defl   
val slope : float = 0.07142857143 
val intercept : float = -0.01666666667 

Chart

Does this line make a good fit? Here is a chart with a couple of samples from this dataset.

On the other hand a classical metric like R^2 can help to answer this question too, especially when it’s extremely simple to add a new column to the dataframe and perform some operations.

The new library is tried out, the lab is completed - everyone is happy ^_^

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16:   
17: 
18: 
19: 
20: 
21: 
22: 
23: 
24: 
25: 
26: 
27: 
bonesreg?prediction <- intercept + slope * bonesreg?p 
bonesreg?residualsq <- bonesreg?prediction - bonesreg?``deflection<mm>`` 
                       |> Series.mapValues (fun x -> x*x)
bonesreg
val it : Frame<int,string> =
       sample type p    deflection<mm> prediction         residualsq           
  0 -> bone        0.98 0.06           0.0533333333333334 4.44444444444441E-05 
  1 -> bone        1.96 0.12           0.123333333333333  1.11111111111111E-05 
 ...-> bone        ...  ...            ...                ...                  

let sdvs = Frame.sdv bonesreg
val sdvs : Series<string,float> =
  sample type    -> <missing>           
  p              -> 1.83341211951923    
  deflection<mm> -> 0.131059782796503   
  prediction     -> 0.130958008537088   
  residualsq     -> 1.72132593164778E-05

// compute the metrics:
let rsquare = let x = sdvs.["prediction"] / sdvs.["deflection<mm>"] in x * x
val rsquare : float = 0.9984475063 
let df = Frame.countRows bonesreg - 2 |> float
val df : float = 4.0 
let tvalue = sqrt (rsquare / (1. - rsquare) * df)
val tvalue : float = 50.71981861 
let se = (Series.sum bonesreg?residualsq) / df |> sqrt
val se : float = 0.005773502692 

Yes, it’s that simple.





namespace System
namespace Deedle
val ds : Frame<int,string>

Full name: Regression.ds

Multiple items
module Frame

from Deedle

--------------------
type Frame<'TRowKey,'TColumnKey (requires equality and equality)> =
  interface IDynamicMetaObjectProvider
  interface INotifyCollectionChanged
  interface IFsiFormattable
  interface IFrame
  new : names:seq<'TColumnKey> * columns:seq<ISeries<'TRowKey>> -> Frame<'TRowKey,'TColumnKey>
  private new : rowIndex:IIndex<'TRowKey> * columnIndex:IIndex<'TColumnKey> * data:IVector<IVector> -> Frame<'TRowKey,'TColumnKey>
  member AddSeries : column:'TColumnKey * series:ISeries<'TRowKey> -> unit
  member AddSeries : column:'TColumnKey * series:seq<'V> -> unit
  member AddSeries : column:'TColumnKey * series:ISeries<'TRowKey> * lookup:Lookup -> unit
  member AddSeries : column:'TColumnKey * series:seq<'V> * lookup:Lookup -> unit
  ...

Full name: Deedle.Frame<_,_>

--------------------
type Frame =
  static member CreateEmpty : unit -> Frame<'R,'C> (requires equality and equality)
  static member FromColumns : cols:Series<'TColKey,Series<'TRowKey,'V>> -> Frame<'TRowKey,'TColKey> (requires equality and equality)
  static member FromColumns : cols:Series<'TColKey,ObjectSeries<'TRowKey>> -> Frame<'TRowKey,'TColKey> (requires equality and equality)
  static member FromColumns : columns:seq<KeyValuePair<'ColKey,ObjectSeries<'RowKey>>> -> Frame<'RowKey,'ColKey> (requires equality and equality)
  static member FromColumns : columns:seq<KeyValuePair<'ColKey,Series<'RowKey,'V>>> -> Frame<'RowKey,'ColKey> (requires equality and equality)
  static member FromColumns : rows:seq<Series<'ColKey,'V>> -> Frame<'ColKey,int> (requires equality)
  static member FromRecords : values:seq<'T> -> Frame<int,string>
  static member FromRecords : series:Series<'K,'R> -> Frame<'K,string> (requires equality)
  static member FromRowKeys : keys:seq<'K> -> Frame<'K,string> (requires equality)
  static member FromRows : rows:Series<'TColKey,Series<'TRowKey,'V>> -> Frame<'TColKey,'TRowKey> (requires equality and equality)
  ...

Full name: Deedle.Frame

--------------------
new : names:seq<'TColumnKey> * columns:seq<ISeries<'TRowKey>> -> Frame<'TRowKey,'TColumnKey>

static member Frame.ReadCsv : path:string * ?hasHeaders:bool * ?inferTypes:bool * ?inferRows:int * ?schema:string * ?separators:string * ?culture:string -> Frame<int,string>
static member Frame.ReadCsv : stream:IO.Stream * ?hasHeaders:bool * ?inferTypes:bool * ?inferRows:int * ?schema:string * ?separators:string * ?culture:string -> Frame<int,string>
val bySample : Frame<(string * int),string>

Full name: Regression.bySample

property Frame.Columns: ColumnSeries<int,string>
val groupRowsByString : column:'a -> frame:Frame<'b,'a> -> Frame<(string * 'b),'a> (requires equality and equality)

Full name: Deedle.Frame.groupRowsByString

property Frame.Columns: ColumnSeries<(string * int),string>
val meanLevel : keySelector:('R -> 'a) -> frame:Frame<'R,'C> -> Frame<'a,'C> (requires equality and equality and equality)

Full name: Deedle.Frame.meanLevel

module Pair

from Deedle

val get1Of2 : v:'a * 'b -> 'a

Full name: Deedle.Pair.get1Of2

val bones : Frame<int,string>

Full name: Regression.bones

val nest : frame:Frame<('R1 * 'R2),'C> -> Series<'R1,Frame<'R2,'C>> (requires equality and equality and equality)

Full name: Deedle.Frame.nest

val deflections : Series<int,float>

Full name: Regression.deflections

Multiple items
module Series

from Deedle

--------------------
type Series<'K,'V (requires equality)> =
  interface IFsiFormattable
  interface ISeries<'K>
  new : pairs:seq<KeyValuePair<'K,'V>> -> Series<'K,'V>
  new : keys:seq<'K> * values:seq<'V> -> Series<'K,'V>
  new : index:IIndex<'K> * vector:IVector<'V> * vectorBuilder:IVectorBuilder * indexBuilder:IIndexBuilder -> Series<'K,'V>
  member Aggregate : aggregation:Aggregation<'K> * observationSelector:Func<DataSegment<Series<'K,'V>>,KeyValuePair<'TNewKey,OptionalValue<'R>>> -> Series<'TNewKey,'R> (requires equality)
  member Aggregate : aggregation:Aggregation<'K> * keySelector:Func<DataSegment<Series<'K,'V>>,'TNewKey> * valueSelector:Func<DataSegment<Series<'K,'V>>,'R> -> Series<'TNewKey,'R> (requires equality)
  member Append : otherSeries:Series<'K,'V> -> Series<'K,'V>
  member AsyncMaterialize : unit -> Async<Series<'K,'V>>
  override Equals : another:obj -> bool
  ...

Full name: Deedle.Series<_,_>

--------------------
type Series =
  static member ofNullables : values:seq<Nullable<'a0>> -> Series<int,'a0> (requires default constructor and value type and 'a0 :> ValueType)
  static member ofObservations : observations:seq<'a0 * 'a1> -> Series<'a0,'a1> (requires equality)
  static member ofOptionalObservations : observations:seq<'K * OptionalValue<'a1>> -> Series<'K,'a1> (requires equality)
  static member ofValues : values:seq<'a0> -> Series<int,'a0>

Full name: Deedle.FSharpSeriesExtensions.Series

--------------------
new : pairs:seq<Collections.Generic.KeyValuePair<'K,'V>> -> Series<'K,'V>
new : keys:seq<'K> * values:seq<'V> -> Series<'K,'V>
new : index:Indices.IIndex<'K> * vector:IVector<'V> * vectorBuilder:Vectors.IVectorBuilder * indexBuilder:Indices.IIndexBuilder -> Series<'K,'V>

val mean : series:Series<'K,'V> -> 'V (requires equality and member ( + ) and member DivideByInt and member get_Zero)

Full name: Deedle.Series.mean

val dropMissing : series:Series<'K,'T> -> Series<'K,'T> (requires equality)

Full name: Deedle.Series.dropMissing

val fillMissing : direction:Direction -> series:Series<'K,'T> -> Series<'K,'T> (requires equality)

Full name: Deedle.Series.fillMissing

type Direction =
  | Backward = 0
  | Forward = 1

Full name: Deedle.Direction

Direction.Forward: Direction = 1
val linearRegression : xs:seq<float> -> ys:seq<float> -> float * float

Full name: Regression.linearRegression

 Find slope and intercept with linear regression

val xs : seq<float>
val ys : seq<float>
let x, y, xx, xy, n =
        Seq.zip xs ys
        |> Seq.fold (fun (xsum, ysum, xxsum, xysum, n) (x, y) ->
            xsum+x, ysum+y, xxsum+x*x, xysum+x*y, n+1.) (0.,0.,0.,0.,0.)
    let slope = (n * xy - x * y) / (n * xx - x * x)
    let intercept = (y - slope * x) / n
    slope, intercept
val bonesreg : Frame<int,string>

Full name: Regression.bonesreg

val dropSparseRows : frame:Frame<'R,'C> -> Frame<'R,'C> (requires equality and equality)

Full name: Deedle.Frame.dropSparseRows

val load : seq<float>

Full name: Regression.load

val values : series:Series<'K,'T> -> seq<'T> (requires equality)

Full name: Deedle.Series.values

val defl : seq<float>

Full name: Regression.defl

val slope : float

Full name: Regression.slope

val intercept : float

Full name: Regression.intercept

val mapValues : f:('T -> 'R) -> series:Series<'K,'T> -> Series<'K,'R> (requires equality)

Full name: Deedle.Series.mapValues

val x : float
val sdvs : Series<string,float>

Full name: Regression.sdvs

val sdv : frame:Frame<'R,'C> -> Series<'C,float> (requires equality and equality)

Full name: Deedle.Frame.sdv

val rsquare : float

Full name: Regression.rsquare

val df : float

Full name: Regression.df

val countRows : frame:Frame<'R,'C> -> int (requires equality and equality)

Full name: Deedle.Frame.countRows

Multiple items
val float : value:'T -> float (requires member op_Explicit)

Full name: Microsoft.FSharp.Core.Operators.float

--------------------
type float<'Measure> = float

Full name: Microsoft.FSharp.Core.float<_>

--------------------
type float = Double

Full name: Microsoft.FSharp.Core.float

val tvalue : float

Full name: Regression.tvalue

val sqrt : value:'T -> 'U (requires member Sqrt)

Full name: Microsoft.FSharp.Core.Operators.sqrt

val se : float

Full name: Regression.se

val sum : series:Series<'K,'V> -> 'V (requires equality and member ( + ) and member get_Zero)

Full name: Deedle.Series.sum

  1. to be honest, still don’t get why anyone would need that, maybe that’s the point where being a PhD helps?