F# | Jamie Dixon's Home

Consuming Azure ML With F#

September 16, 2014 2 Comments

(This post is a continuation of this one)

So with a model that works well enough, I selected only that model and saved it

Created a new experiment and used that model with the base data. I then marked the project columns as the input and the score as the output (green and blue circle respectively)

After running it, I published it as a web service

And voila, an endpoint ready to go. I then took the auto generated script and opened up a new Visual Studio F# project to use it. The problem was that this is the data structure that the model needs

FeatureVector = new Dictionary<string, string>() 
    {
        { "Precinct", "0" },
        { "VRN", "0" },
        { "VRstatus", "0" },
        { "VRlastname", "0" },
        { "VRfirstname", "0" },
        { "VRmiddlename", "0" },
        { "VRnamesufx", "0" },
        { "VRstreetnum", "0" },
        { "VRstreethalfcode", "0" },
        { "VRstreetdir", "0" },
        { "VRstreetname", "0" },
        { "VRstreettype", "0" },
        { "VRstreetsuff", "0" },
        { "VRstreetunit", "0" },
        { "VRrescity", "0" },
        { "VRstate", "0" },
        { "Zip Code", "0" },
        { "VRfullresstreet", "0" },
        { "VRrescsz", "0" },
        { "VRmail1", "0" },
        { "VRmail2", "0" },
        { "VRmail3", "0" },
        { "VRmail4", "0" },
        { "VRmailcsz", "0" },
        { "Race", "0" },
        { "Party", "0" },
        { "Gender", "0" },
        { "Age", "0" },
        { "VRregdate", "0" },
        { "VRmuni", "0" },
        { "VRmunidistrict", "0" },
        { "VRcongressional", "0" },
        { "VRsuperiorct", "0" },
        { "VRjudicialdistrict", "0" },
        { "VRncsenate", "0" },
        { "VRnchouse", "0" },
        { "VRcountycomm", "0" },
        { "VRschooldistrict", "0" },
        { "11/6/2012", "0" },
        { "Voted Ind", "0" },
    },
    GlobalParameters = new Dictionary<string, string>() 
    {
    }
};

And since I am only using 6 of the columns, it made sense to reload the Wake County Voter Data with just the needed columns. I went back to the original CSV and did that. Interestingly, I could not set the original dataset as the publish input so I added a project column module that does nothing

With that in place, I republished the service and opened Visual Studio. I decided to start with a script. I was struggling though the async when Tomas P helped me on Stack Overflow here. I’ll say it again, the F# community is tops. In any event, here is the initial script:


#r @"C:\Program Files (x86)\Reference Assemblies\Microsoft\Framework\.NETFramework\v4.5\System.Net.Http.dll"
#r @"..\packages\Microsoft.AspNet.WebApi.Client.5.2.2\lib\net45\System.Net.Http.Formatting.dll"

open System
open System.Net.Http
open System.Net.Http.Headers
open System.Net.Http.Formatting
open System.Collections.Generic

type scoreData = {FeatureVector:Dictionary<string,string>;GlobalParameters:Dictionary<string,string>}
type scoreRequest = {Id:string; Instance:scoreData}

let invokeService () = async {
    let apiKey = ""
    let uri = "https://ussouthcentral.services.azureml.net/workspaces/19a2e623b6a944a3a7f07c74b31c3b6d/services/f51945a42efa42a49f563a59561f5014/score"
    use client = new HttpClient()
    client.DefaultRequestHeaders.Authorization <- new AuthenticationHeaderValue("Bearer",apiKey)
    client.BaseAddress <- new Uri(uri)

    let input = new Dictionary<string,string>()
    input.Add("Zip Code","27519")
    input.Add("Race","W")
    input.Add("Party","UNA")
    input.Add("Gender","M")
    input.Add("Age","45")
    input.Add("Voted Ind","1")

    let instance = {FeatureVector=input; GlobalParameters=new Dictionary<string,string>()}
    let scoreRequest = {Id="score00001";Instance=instance}

    let! response = client.PostAsJsonAsync("",scoreRequest) |> Async.AwaitTask
    let! result = response.Content.ReadAsStringAsync() |> Async.AwaitTask

    if response.IsSuccessStatusCode then
        printfn "%s" result
    else
        printfn "FAILED: %s" result
    response |> ignore
}

invokeService() |> Async.RunSynchronously

Unfortunately, when I run it, it fails. Below is the Fiddler trace:

So it looks like the Json Serializer is postpending the “@” symbol. I changed the records to types and voila:

You can see the final script here.

So then throwing in some different numbers.

A millennial: ["27519","W","D","F","25","1","1","0.62500011920929"]
A senior citizen: ["27519","W","D","F","75","1","1","0.879632294178009"]

I wonder why social security never gets cut?

In any event, just to check the model:

A 15 year old: ["27519","W","D","F","15","1","0","0.00147285079583526"]

Filed under Analytics, Azure ML, F#

Sql Saturday and MVP Monday

September 9, 2014 Leave a comment

Thanks to everyone who came to my session on F# Type Providers. The code is found here.

Also, my article on the Eject-A-Bed was selected for MVP Mondays. You can see a link here.

Filed under F#, Netduino

Fun with Statistics and Charts

September 2, 2014 1 Comment

I am preparing my Raleigh Code Camp submission ‘Nerd Dinner With Brains” this weekend. If you are not familiar, Nerd Dinner is the canonical example of a MVC application and is very familiar to Web Devs who want to learn MVC the Microsoft way. You can see the walkthrough here. For everything that Nerd Dinner is, it is not … smart. There is no business rules outside of some basic input validation, which is pretty representative of many “Boring Line Of Business Applications (BLOBAs according to Scott Waschlan). Not coincidently, the lack of business logic is the biggest reason many BLOBAs don’t have many unit tests –> if all you are doing is wire framing a database, what business logic needs to be tested?

The talk is going to take the Nerd Diner wireframe and inject some analytics to the application. To that end, I first considered the person who is attending the dinner. All we know about them is their name and possibly their location. So what can a name tell you? Turns out, plenty.

As I showed in this post, there is a great source of the number of names given by gender, yearOfBrith, and stateOfBirth from the US census. Picking up where that post left off, I loaded in the entire data set into memory.

My first question was, “given a name, can I tell what gender the person is?” This is very straight forward to calculate.

 1 let genderSearch name = 
 2     let nameFilter = usaData
 3                         |> Seq.filter(fun r -> r.Mary = name)
 4                         |> Seq.groupBy(fun r -> r.F)
 5                         |> Seq.map(fun (n,a) -> n,a |> Seq.sumBy(fun (r) -> r.``14``)) 
 6  
 7     let nameSum = nameFilter |> Seq.sumBy(fun (n,c) -> c)
 8     nameFilter 
 9         |> Seq.map(fun (n,c) -> n, c, float c/float nameSum) 
10         |> Seq.toArray
11 
12 genderSearch "James"
13

And the REPL shows me that is is very likely that “James” is a male:

I can then set up in the web.config file a confidence point where there name is a male/female, I am thinking 75%. Once we have that, the app can respond differently. Perhaps we have a product-placement advertisement that becomes a male-focused if we are reasonably certain that the user is a male. Perhaps we can be more subtle and change the theme of the site, or the page navigation, to induce the person to do additional things on the site.

In any event, I then wanted to tackle age. I spun up some code to isolate a person’s age

 1 let ageSearch name =
 2     let nameFilter = usaData
 3                         |> Seq.filter(fun r -> r.Mary = name)
 4                         |> Seq.groupBy(fun r -> r.``1910``)
 5                         |> Seq.map(fun (n,a) -> n,a |> Seq.sumBy(fun (r) -> r.``14``)) 
 6                         |> Seq.toArray
 7     let nameSum = nameFilter |> Seq.sumBy(fun (n,c) -> c)
 8     nameFilter 
 9         |> Seq.map(fun (n,c) -> n, c, float c/float nameSum) 
10         |> Seq.toArray

I had no idea if names have a certain age connotation so I decided to do some basic charting. Isaac Abraham pointed me to FSharp.Chart which is a great way to do some basic charting for discovery.

1 let chartData = ageSearch "James"
2                     |> Seq.map(fun (y,c,p) -> y, c)
3                     |> Seq.sortBy(fun (y,c) -> y)
4     
5 Chart.Line(chartData).ShowChart()

And sure enough, the name “James” has a real ebb and flow for its popularity.

so if the user has a name of “James”, you can make a reasonable assumption they are male and probably born before 1975. Cue up the Van Halen!

And yes, because I had to:

1 let chartData = ageSearch "Britney"
2                     |> Seq.map(fun (y,c,p) -> y, c)
3                     |> Seq.sortBy(fun (y,c) -> y)

Kinda does match her career, no?

Anyway, back to the task at hand. In terms of analytics, I want to be a bit more precise then eyeballing a chart. I started with the following code:

 1 ageSearch "James"
 2     |> Seq.map(fun (y,c,p) -> float c)
 3     |> Seq.average
 4 
 5 ageSearch "James"
 6     |> Seq.map(fun (y,c,p) -> float c)
 7     |> Seq.min
 8 
 9 ageSearch "James"
10     |> Seq.map(fun (y,c,p) -> float c)
11     |> Seq.max
12

With these basic statistics out of the way, I then wanted to look at when the name was no longer popular. I decided to use 1 standard deviation away from the average to determine an outlier. First the standard deviation:

 1 let variance (source:float seq) =
 2     let mean = Seq.average source
 3     let deltas = Seq.map(fun x -> pown(x-mean) 2) source
 4     Seq.average deltas
 5 
 6 let standardDeviation(values:float seq) =
 7     sqrt(variance(values))
 8 
 9 ageSearch "James"
10     |> Seq.map(fun (y,c,p) -> float c)
11     |> standardDeviation
12 
13 let standardDeviation' = ageSearch "James"
14                             |> Seq.map(fun (y,c,p) -> float c)
15                             |> standardDeviation
16 
17 let average = ageSearch "James"
18                 |> Seq.map(fun (y,c,p) -> float c)
19                 |> Seq.average
20 
21 let attachmentPoint = average+standardDeviation'

And then I can get the last year that the name was within 1 standard deviation above the average (greater than 71,180 names given):

1 
2 let popularYears = ageSearch "James"
3                         |> Seq.map(fun (y,c,p) -> y, float c)
4                         |> Seq.filter(fun (y,c) -> c > attachmentPoint)
5                         |> Seq.sortBy(fun (y,c) -> y)
6                         |> Seq.last

So “James” is very likely a male and likely born before 1964. Cue up the Pink Floyd!

The last piece was the state of birth –> can I guess the state of birth for a user? I first looked at the states on a plot

1 let chartData' = stateSearch "James"
2                     |> Seq.map(fun (s,c,p) -> s,c)
3     
4 Chart.Column(chartData').ShowChart()
5

Nothing really stands out at me –> states with the most births have the most names. I could do an academic exercise of seeing what states favor certain names, but that does not help me with Nerd Dinner in guessing the state of birth when given a name.

I pressed on to look at the top 10 states:

 1 let topTenStates = stateSearch "James"
 2                     |> Seq.sortBy(fun (s,c,p) -> -c-1)
 3                     |> Seq.take 10
 4 
 5 let topTenTotal = topTenStates 
 6                     |> Seq.sumBy(fun (s,c,p) -> c)
 7 let total = stateSearch "James"
 8                 |> Seq.sumBy(fun (s,c,p) -> c)
 9 
10 float topTenTotal/float total

So 50% of “James” were born in 10 states. Again, I am not sure there is any actionable information here. For example, if a majority of “James” were born in MI, I might have something (cue up the Bob Seger).

Interestingly, there are certain number of names where the state of birth does matter. For example, consider “Jose”:

Unsurprisingly, the two states are CA and TX. Just using James and Jose as an example:

James is a male born before 1964
Jose is a male born before 2008 in either TX or CA

As an academic exercise, we could construct a random forest to find the names with the greatest state affinity. However, that won’t help us on Nerd Dinner so I am leaving that out for another day.

This analysis does not account for a host of factors (person not born in the USA, nicknames, etc..), but it is still better than the nothing that Nerd Dinner currently has. This analysis is not particular sophisticated but I often find that even the most basic statistics can be very powerful if used correctly. That will be the next part of the talk…

Filed under Analytics, F#

Consuming and Analyzing Census Data Using F#

August 19, 2014 3 Comments

As part of my Nerd Dinner refactoring, I wanted to add the ability to guess a person’s age and gender based on their name. I did a quick search on the internet and the only place that I found that has an API is here and it doesn’t have everything I am looking for. Fortunately, the US Census website has some flat files with the kind of data I am looking for here.

I grabbed the data and pumped it into Azure Blob Storage here. You can swap out the state code to get each dataset. I then loaded in a list of State Codes found here that match to the file names.

I then fired up Visual Studio and created a new FSharp project. I added FSharp.Data to use a Type Provider to access the data. I don’t need to install the Azure Storage .dlls b/c the blobs are public and I just have to read the file

Once Nuget was done with its magic, I opened up the script file, pointed to the newly-installed FSharp.Data, and added a reference to the datasets on blob storage:

#r "../packages/FSharp.Data.2.0.9/lib/portable-net40+sl5+wp8+win8/FSharp.Data.dll"
open FSharp.Data


type censusDataContext = CsvProvider<"https://portalvhdspgzl51prtcpfj.blob.core.windows.net/censuschicken/AK.TXT">
type stateCodeContext = CsvProvider<"https://portalvhdspgzl51prtcpfj.blob.core.windows.net/censuschicken/states.csv">

(Note that I am going add FSharp as a language to my Live Writer code snippet add-in at a later date)

In any event, I then printed out all of the codes to see what it looks like:

let stateCodes =  stateCodeContext.Load("https://portalvhdspgzl51prtcpfj.blob.core.windows.net/censuschicken/states.csv");
stateCodes.Rows |> Seq.iter(fun r -> printfn "%A" r)

And by changing the lambda slightly like so,

stateCodes.Rows |> Seq.iter(fun r -> printfn "%A" r.Abbreviation)

I get all of the state codes

I then tested the census data with code and results are expected

let arkansasData = censusDataContext.Load("https://portalvhdspgzl51prtcpfj.blob.core.windows.net/censuschicken/AK.TXT");
arkansasData.Rows |> Seq.iter(fun r -> printfn "%A" r)

So then I created a method to load all of the state census data and giving me the length of the total:

let stateCodes =  stateCodeContext.Load("https://portalvhdspgzl51prtcpfj.blob.core.windows.net/censuschicken/states.csv");
let usaData = stateCodes.Rows 
                |> Seq.collect(fun r -> censusDataContext.Load(System.String.Format("https://portalvhdspgzl51prtcpfj.blob.core.windows.net/censuschicken/{0}.TXT",r.Abbreviation)).Rows)
                |> Seq.length

Since this is a I/O bound operation, it made sense to load the data asynchronously, which speeded things up considerably. You can see my question over on Stack Overflow here and the resulting code takes about 50% of the time on a my dual-processor machine:

stopwatch.Start()
let fetchStateDataAsync(stateCode:string)=
    async{
        let uri = System.String.Format("https://portalvhdspgzl51prtcpfj.blob.core.windows.net/censuschicken/{0}.TXT",stateCode)
        let! stateData =  censusDataContext.AsyncLoad(uri)
        return stateData.Rows
    }


let usaData' = stateCodes.Rows
                    |> Seq.map(fun r -> fetchStateDataAsync(r.Abbreviation))
                    |> Async.Parallel
                    |> Async.RunSynchronously
                    |> Seq.collect id
                    |> Seq.length
stopwatch.Stop()
printfn "Parallel: %A" stopwatch.Elapsed.Seconds

With the data in hand, it was time to analyze the data to see if there is anything we can do. Since 23 seconds is a bit too long to wait for a page load ( Smile ), I will need to put the 5.5 million records into a format that can be easily searched. Thinking what we want is:

Given a name, what is the gender?

Given a name, what is the age?

Given a name, what is their state of birth?

Also, since we have their current location, we can also input the name and location and answer those questions. If we make the assumption that their location is the same as their birth state, we can narrow down the list even further.

In any event, I first added a GroupBy to the name:

let nameSum = usaData' 
                |> Seq.groupBy(fun r -> r.Mary)
                |> Seq.toArray

And then I summed up the counts of the names

let nameSum = usaData' 
                |> Seq.groupBy(fun r -> r.Mary)
                |> Seq.map(fun (n,a) -> n,a |> Seq.sumBy(fun (r) -> r.``14``)) 
                |> Seq.toArray

And then the total in the set:

let totalNames = nameSum |> Seq.sumBy(fun (n,c) -> c)

And then applied a simple average and sorted it descending

let nameAverage = nameSum 
                    |> Seq.map(fun (n,c) -> n,c,float c/ float totalNames)
                    |> Seq.sortBy(fun (n,c,a) -> -a - 1.)
                    |> Seq.toArray

So I feel really special that my parents gave me the most popular name in the US ever…

And focusing back to the task on hand, I want to determine the probability that a person is male or female based on their name:

let nameSearch = usaData'
                    |> Seq.filter(fun r -> r.Mary = "James")
                    |> Seq.groupBy(fun r -> r.F)
                    |> Seq.map(fun (n,a) -> n,a |> Seq.sumBy(fun (r) -> r.``14``)) 
                    |> Seq.toArray

So 18196 parents thought is would be a good idea to name their daughter ‘James’. I created a quick function like so:

let nameSearch' name = 
    let nameFilter = usaData'
                        |> Seq.filter(fun r -> r.Mary = name)
                        |> Seq.groupBy(fun r -> r.F)
                        |> Seq.map(fun (n,a) -> n,a |> Seq.sumBy(fun (r) -> r.``14``)) 

    let nameSum = nameFilter |> Seq.sumBy(fun (n,c) -> c)
    nameFilter 
        |> Seq.map(fun (n,c) -> n, c, float c/float nameSum) 
        |> Seq.toArray

nameSearch' "James"

So if I see the name “James”, there is a 99% chance it is a male. This can lead to a whole host of questions like variance of names, names that are closest to gender neutral, etc…. Leaving those questions to another day, I now have something I can put into Nerd Dinner. Now, if there was only a way to handle nicknames and friendly names….

You can see the full code here.

Filed under F#, Statistics

Controlling Servos Using Netdunio and Phidgets

August 12, 2014 1 Comment

As part of the Terminator program I am creating, I need a way of controlling servos to point the laser (and then gun) and different targets. I decided to create a POC project and evaluate two different ways of controlling the servos. As step one, I purchased a pan and tilt chassis from here.

After playing with the servos from the kit, I decided to use my old stand-by servos that had a much higher quality and whose PWM signals I already know how to use. With the chassis done, I needed a laser pointer so I figured why not get a shark with fricken laser?

I found one here.

So with the servos and laser ready to go, it was time to code. I started with Netduninos:

    public class Program
    {
        private const uint TILT_SERVO_STRAIGHT = 1500;
        private const uint TILT_SERVO_MAX_UP = 2000;
        private const uint TILT_SERVO_MAX_DOWN = 1000;
        private const uint PAN_SERVO_STRAIGHT = 1500;
        private const uint PAN_SERVO_MAX_LEFT = 1000;
        private const uint PAN_SERVO_MAX_RIGHT = 2000;

        private static PWM _tiltServo = null;
        private static PWM _panServo = null;
        private static uint _tiltServoCurrentPosition = 0;
        private static uint _panServoCurrentPosition = 0;

        public static void Main()
        {
            SetUpServos();
            InputPort button = new InputPort(Pins.ONBOARD_BTN, false, Port.ResistorMode.Disabled);
            while (true)
            {
                if (button.Read())
                {
                    MoveServo();
                }
            }

        }

        private static void SetUpServos()
        {
            uint period = 20000;
            _tiltServoCurrentPosition = TILT_SERVO_STRAIGHT;
            _panServoCurrentPosition = PAN_SERVO_STRAIGHT;

            _tiltServo = new PWM(PWMChannels.PWM_PIN_D3, period, _tiltServoCurrentPosition, PWM.ScaleFactor.Microseconds, false);
            _tiltServo.Start();
            _panServo = new PWM(PWMChannels.PWM_PIN_D5, period, _panServoCurrentPosition, PWM.ScaleFactor.Microseconds, false);
            _panServo.Start();
        }

        private static void MoveServo()
        {
            _panServo.Duration = PAN_SERVO_MAX_LEFT;
            Thread.Sleep(2000);
            _panServo.Duration = PAN_SERVO_MAX_RIGHT;
            Thread.Sleep(2000);
            _panServo.Duration = PAN_SERVO_STRAIGHT;
            Thread.Sleep(2000);
            _tiltServo.Duration = TILT_SERVO_MAX_UP;
            Thread.Sleep(2000);
            _tiltServo.Duration = TILT_SERVO_MAX_DOWN;
            Thread.Sleep(2000);
            _tiltServo.Duration = TILT_SERVO_STRAIGHT;
        }

    }

And sure enough the servos are behaving as expected

I then implemented a similar app using Phidgets. Because the code is being executed on the PC, I could use F# to code (It does not look like the Netdunino/Microframework supports F#?)

open System
open Phidgets

let _servoController = new AdvancedServo()
let mutable _isServoControllerReady = false

let servoController_Attached(args:Events.AttachEventArgs) =
    let servoController =  args.Device  :?> AdvancedServo
    servoController.servos.[0].Engaged <- true
    servoController.servos.[7].Engaged <- true
    _isServoControllerReady <- true

[<EntryPoint>]
let main argv = 
    _servoController.Attach.Add(servoController_Attached)
    _servoController.``open``()

    while true do
        if _isServoControllerReady = true then
            _servoController.servos.[0].Position<- 100.
            _servoController.servos.[7].Position<- 100.

    Console.ReadKey() |> ignore
    printfn "%A" argv
    0

The choice then becomes using the Netduino or the Phidgets with my Kinect program. I decided to defer the decision and use an interface for now.

type IWeaponsSystem =
    abstract member Activate: unit -> unit
    abstract member AquireTarget : float*float -> bool
    abstract member Fire: int -> bool

My decision about using Phidgets or Netduino is a series of trade-offs. I can code Phidgets in C# or F# but I have to code Netduino in C#. I would prefer to do this in F# so that makes me learn towards Phidgets. I can put the Netduino anywhere and have it communicate via an Ethernet signal but I have to have the Phidgets wired to the PC. Since the targeting system needs to be near the Kinect and the Kinect has to be tethered to the PC also, there is no real advantage of using the mobile Netduino. Finally, the Phidgets API handles all communication to the servo control board for me, with the Netduino I would have to hook up a router to the Netduino and write the Ethernet communication code. So I am leaning towards Phidgets, but since I am not sure, the interface allows me to swap in the Netduino at a later point without changing any code. Love me some O in SOLID…

Up next, integrating the targeting system into the Terminator program.

Filed under F#, Netduino, Phidgets

Neural Network Part 3: Perceptrons

August 5, 2014 1 Comment

I went back to my code for building a Perceptron and I made some changes. I realized that although McCaffrey combines the code together, there are actually two actions for the perceptron: training and predicting. I created a diagram to help me keep the functions that I need for each in mind:

I also skeletoned out some data structures that I think I need:

With the base diagrams out of the way, I created different data structures that were tailored to each action. These are a bit different than the diagrams –> I didn’t go back and update the diagrams because the code is where you would look to see how the system works:

type observation = {xValues:float List}
type weightedObservation = {xws:(float*float) List}
type confirmedObservation = {observation:observation;yExpected:float}
type weightedConfirmedObservation = {weightedObservation:weightedObservation;yExpected:float}
type neuronInput = {weightedObservation:weightedObservation;bias:float}

type cycleTrainingInput = {weightedConfirmedObservation:weightedConfirmedObservation;bias:float;alpha:float}
type adjustmentInput = {weightedConfirmedObservation:weightedConfirmedObservation;bias:float;alpha:float;yActual:float}
type adjustmentOutput = {weights:float List; bias:float}
type rotationTrainingInput = {confirmedObservations:confirmedObservation List;weights:float List;bias:float;alpha:float}
type trainInput = {confirmedObservations:confirmedObservation List; weightSeedValue:float;biasSeedValue:float;alpha:float; maxEpoches:int}

type cyclePredictionInput = {weightedObservation:weightedObservation;bias:float}
type rotationPredictionInput = {observations:observation List;weights:float List;bias:float}
type predictInput = {observations:observation List;weights:float List;bias:float}

Note that I am composing data structures with the base being an observation. The observation is a list of different xValues for a given, well, observation. The weighted observation is the XValue paired with the perceptron weights. The confirmedObservation is for training –> given an observation, what was the actual output?

With the data structures out of the way, I went to the Perceptron and added in the basic functions for creating seed values:

    member this.initializeWeights(xValues, randomSeedValue) =
        let lo = -0.01
        let hi = 0.01
        let xWeight = (hi-lo) * randomSeedValue + lo
        xValues |> Seq.map(fun w -> xWeight)

    member this.initializeBias(randomSeedValue) = 
        let lo = -0.01
        let hi = 0.01
        (hi-lo) * randomSeedValue + lo

Since I was doing TDD, here are the unit tests I used for these functions:

        [TestMethod]
        public void initializeWeightsUsingHalfSeedValue_ReturnsExpected()
        {
            var weights = _perceptron.initializeWeights(_observation.xValues, .5);
            var weightsList = new List<double>(weights);

            var expected = 0.0;
            var actual = weightsList[0];
            Assert.AreEqual(expected, actual);
        }

        [TestMethod]
        public void initializeWeightsUsingLessThanHalfSeedValue_ReturnsExpected()
        {
            var weights = _perceptron.initializeWeights(_observation.xValues, .4699021627);
            var weightsList = new List<double>(weights);

            var expected = -0.00060;
            var actual = Math.Round(weightsList[0],5);
            Assert.AreEqual(expected, actual);
        }

        [TestMethod]
        public void initializeBiasesUsingHalfSeedValue_ReturnsExpected()
        {
            var expected = 0.0;
            var actual = _perceptron.initializeBias(.5);
            Assert.AreEqual(expected, actual);
        }

        [TestMethod]
        public void initializeBiasesUsingLessThanHalfSeedValue_ReturnsExpected()
        {
            var expected = -0.00060;
            var bias = _perceptron.initializeBias(.4699021627);
            var actual = Math.Round(bias, 5);
            Assert.AreEqual(expected, actual);
        }

        [TestMethod]
        public void initializeBiasesUsingGreaterThanHalfSeedValue_ReturnsExpected()
        {
            var expected = 0.00364;
            var bias = _perceptron.initializeBias(.6820621978);
            var actual = Math.Round(bias,5);
            Assert.AreEqual(expected, actual);
        }

I then created a base neuron and activation function that would work for both training and predicting:

    member this.runNeuron(input:neuronInput) = 
        let xws = input.weightedObservation.xws
        let output = xws 
                    |> Seq.map(fun (xValue,xWeight) -> xValue*xWeight) 
                    |> Seq.sumBy(fun x -> x)
        output + input.bias

    member this.runActivation(input) = 
        if input < 0.0 then -1.0 else 1.0

        [TestMethod]
        public void runNeuronUsingNormalInput_ReturnsExpected()
        {
            var expected = -0.0219;
            var perceptronOutput = _perceptron.runNeuron(_neuronInput);
            var actual = Math.Round(perceptronOutput, 4);
            Assert.AreEqual(expected, actual);
        }

        [TestMethod]
        public void runActivationUsingNormalInput_ReturnsExpected()
        {
            var expected = -1;
            var actual = _perceptron.runActivation(-0.0219);
            Assert.AreEqual(expected, actual);
        }

I then created the functions for training –> specifically to return adjusted weights and biases based on the result of the activation function

    member this.calculateWeightAdjustment(xValue, xWeight, alpha, delta) =
        match delta > 0.0, xValue >= 0.0 with
            | true,true -> xWeight - (alpha * abs(delta) * xValue)
            | false,true -> xWeight + (alpha * abs(delta) * xValue)
            | true,false -> xWeight - (alpha * abs(delta) * xValue)
            | false,false -> xWeight + (alpha * abs(delta) * xValue)

    member this.calculateBiasAdjustment(bias, alpha, delta) =
        match delta > 0.0 with
            | true -> bias - (alpha * abs(delta))
            | false -> bias + (alpha * abs(delta))

    member this.runAdjustment (input:adjustmentInput) =
        match input.weightedConfirmedObservation.yExpected = input.yActual with
            | true -> 
                let weights = input.weightedConfirmedObservation.weightedObservation.xws |> Seq.map(fun (x,w) -> w)
                let weights' = new List<float>(weights)
                {adjustmentOutput.weights=weights';adjustmentOutput.bias=input.bias}
            | false ->
                let delta = input.yActual - input.weightedConfirmedObservation.yExpected
                let weights' =  input.weightedConfirmedObservation.weightedObservation.xws 
                                    |> Seq.map(fun (xValue, xWeight) -> 
                                        this.calculateWeightAdjustment(xValue,xWeight,input.alpha,delta))
                                    |> Seq.toList
                let weights'' = new List<float>(weights')                                         
                let bias' = this.calculateBiasAdjustment(input.bias,input.alpha,delta)
                {adjustmentOutput.weights=weights'';adjustmentOutput.bias=bias'}

        [TestMethod]
        public void calculateWeightAdjustmentUsingPositiveDelta_ReturnsExpected()
        {
            var xValue = 1.5;
            var xWeight = .00060;
            var delta = 2;

            var weightAdjustment = _perceptron.calculateWeightAdjustment(xValue, xWeight, _alpha, delta);
            var actual = Math.Round(weightAdjustment, 4);
            var expected = -.0024;
            Assert.AreEqual(expected, actual);
        }

        [TestMethod]
        public void calculateWeightAdjustmentUsingNegativeDelta_ReturnsExpected()
        {
            var xValue = 1.5;
            var xWeight = .00060;
            var delta = -2;

            var weightAdjustment = _perceptron.calculateWeightAdjustment(xValue, xWeight, _alpha, delta);
            var actual = Math.Round(weightAdjustment, 5);
            var expected = .0036;
            Assert.AreEqual(expected, actual);
        }

        [TestMethod]
        public void calculateBiasAdjustmentUsingPositiveDelta_ReturnsExpected()
        {
            var bias = 0.00364;
            var delta = 2;
            var expected = .00164;
            var actual = _perceptron.calculateBiasAdjustment(bias, _alpha, delta);
            Assert.AreEqual(expected, actual);
        }

        [TestMethod]
        public void calculateBiasAdjustmentUsingNegativeDelta_ReturnsExpected()
        {
            var bias = 0.00364;
            var delta = -2;
            var expected = .00564;
            var actual = _perceptron.calculateBiasAdjustment(bias, _alpha, delta);
            Assert.AreEqual(expected, actual);
        }

        [TestMethod]
        public void runAdjustmentUsingMatchingData_ReturnsExpected()
        {

            var adjustmentInput = new adjustmentInput(_weightedConfirmedObservation, _bias, _alpha, -1.0);
            var adjustedWeights = _perceptron.runAdjustment(adjustmentInput);
            var expected = .0065;
            var actual = Math.Round(adjustedWeights.weights[0],4);
            Assert.AreEqual(expected, actual);
        }

        [TestMethod]
        public void runAdjustmentUsingNegativeData_ReturnsExpected()
        {
            weightedConfirmedObservation weightedConfirmedObservation = 
                new NeuralNetworks.weightedConfirmedObservation(_weightedObservation, 1.0);
            var adjustmentInput = new adjustmentInput(weightedConfirmedObservation, _bias, _alpha, -1.0);
            var adjustedWeights = _perceptron.runAdjustment(adjustmentInput);
            var expected = .0125;
            var actual = Math.Round(adjustedWeights.weights[0], 4);
            Assert.AreEqual(expected, actual);
        }

        [TestMethod]
        public void runAdjustmentUsingPositiveData_ReturnsExpected()
        {
            var adjustmentInput = new adjustmentInput(_weightedConfirmedObservation, _bias, _alpha, 1.0);
            var adjustedWeights = _perceptron.runAdjustment(adjustmentInput);
            var expected = .0005;
            var actual = Math.Round(adjustedWeights.weights[0], 4);
            Assert.AreEqual(expected, actual);
        }

With these functions ready, I could run a training cycle for a given observation

    member this.runTrainingCycle (cycleTrainingInput:cycleTrainingInput) =
        let neuronTrainingInput = 
            {neuronInput.weightedObservation=cycleTrainingInput.weightedConfirmedObservation.weightedObservation;
            neuronInput.bias=cycleTrainingInput.bias}
        let neuronResult = this.runNeuron(neuronTrainingInput)
        let activationResult = this.runActivation(neuronResult)
        let adjustmentInput = {weightedConfirmedObservation=cycleTrainingInput.weightedConfirmedObservation;
                                bias=cycleTrainingInput.bias;alpha=cycleTrainingInput.alpha;
                                yActual=activationResult}
        this.runAdjustment(adjustmentInput)

        [TestMethod]
        public void runTrainingCycleUsingNegativeData_ReturnsExpected()
        {
            var cycleTrainingInput = new cycleTrainingInput(_weightedConfirmedObservation, _bias, _alpha);
            var adjustmentOutput = _perceptron.runTrainingCycle(cycleTrainingInput);
            var expected = .0125;
            var actual = Math.Round(adjustmentOutput.weights[0], 4);
            Assert.AreEqual(expected, actual);
        }

        [TestMethod]
        public void runTrainingCycleUsingPositiveData_ReturnsExpected()
        {
            var cycleTrainingInput = new cycleTrainingInput(_weightedConfirmedObservation, _bias, _alpha);
            var adjustmentOutput = _perceptron.runTrainingCycle(cycleTrainingInput);
            var expected = .0065;
            var actual = Math.Round(adjustmentOutput.weights[0], 4);
            Assert.AreEqual(expected, actual);
        }

And then I could run a cycle for each of the observations in the training set, a rotation. I am not happy that I am mutating the weights and biases here, though I am not sure how to fix that. I looked for a Seq.Scan function where the results of a function applied to the 1st element of a Seq is used in the input of the next –> all I could see were examples of threading a collector of int (like Seq.mapi). This will be something I will ask the functional ninjas when I see them again.

    member this.runTrainingRotation(rotationTrainingInput: rotationTrainingInput)=
        let mutable weights = rotationTrainingInput.weights
        let mutable bias = rotationTrainingInput.bias
        let alpha = rotationTrainingInput.alpha
        for i=0 to rotationTrainingInput.confirmedObservations.Count-1 do
            let currentConfirmedObservation = rotationTrainingInput.confirmedObservations.[i]
            let xws = Seq.zip currentConfirmedObservation.observation.xValues weights
            let xws' = new List<(float*float)>(xws)
            let weightedObservation = {xws=xws'}
            let weightedTrainingObservation = {weightedObservation=weightedObservation;yExpected=currentConfirmedObservation.yExpected}
            let cycleTrainingInput = {
                                cycleTrainingInput.weightedConfirmedObservation=weightedTrainingObservation;
                                cycleTrainingInput.bias=bias;
                                cycleTrainingInput.alpha=alpha}
            let cycleOutput = this.runTrainingCycle(cycleTrainingInput)
            weights <- cycleOutput.weights
            bias <- cycleOutput.bias
        {adjustmentOutput.weights=weights; adjustmentOutput.bias=bias}

        [TestMethod]
        public void runTrainingRotationUsingNegativeData_ReturnsExpected()
        {

            var xValues = new List<double>();
            xValues.Add(3.0);
            xValues.Add(4.0);
            var observation = new observation(xValues);
            var yExpected = -1.0;
            var confirmedObservation0 = new confirmedObservation(observation, yExpected);

            xValues = new List<double>();
            xValues.Add(1.5);
            xValues.Add(2.0);
            yExpected = -1.0;
            var confirmedObservation1 = new confirmedObservation(observation, yExpected);

            var trainingObservations = new List<confirmedObservation>();
            trainingObservations.Add(confirmedObservation0);
            trainingObservations.Add(confirmedObservation1);
            
            var weights = new List<double>();
            weights.Add(.0065);
            weights.Add(.0123);

            var rotationTrainingInput = new rotationTrainingInput(trainingObservations, weights, _bias, _alpha);
            var trainingRotationOutput = _perceptron.runTrainingRotation(rotationTrainingInput);

            var expected = -0.09606;
            var actual = Math.Round(trainingRotationOutput.bias, 5);
            Assert.AreEqual(expected, actual);

        }

        [TestMethod]
        public void runTrainingRotationUsingPositiveData_ReturnsExpected()
        {
            var xValues = new List<double>();
            xValues.Add(3.0);
            xValues.Add(4.0);
            var observation = new observation(xValues);
            var yExpected = 1.0;
            var confirmedObservation0 = new confirmedObservation(observation, yExpected);

            xValues = new List<double>();
            xValues.Add(1.5);
            xValues.Add(2.0);
            yExpected = 1.0;
            var confirmedObservation1 = new confirmedObservation(observation, yExpected);

            var trainingObservations = new List<confirmedObservation>();
            trainingObservations.Add(confirmedObservation0);
            trainingObservations.Add(confirmedObservation1);


            var weights = new List<double>();
            weights.Add(.0065);
            weights.Add(.0123);

            var rotationTrainingInput = new rotationTrainingInput(trainingObservations, weights, _bias, _alpha);
            var trainingRotationOutput = _perceptron.runTrainingRotation(rotationTrainingInput);

            var expected = -.09206;
            var actual = Math.Round(trainingRotationOutput.bias, 5);
            Assert.AreEqual(expected, actual);

        }

With the rotation done, I could write the train function which runs rotations for N number of times to tune the weights and biases:

    member this.train(trainInput:trainInput) = 
        let currentObservation = trainInput.confirmedObservations.[0].observation
        let weights = this.initializeWeights(currentObservation.xValues,trainInput.weightSeedValue)
        let weights' = new List<float>(weights)
        let mutable bias = this.initializeBias(trainInput.biasSeedValue)
        let alpha = trainInput.alpha
        for i=0 to trainInput.maxEpoches do
            let rotationTrainingInput={rotationTrainingInput.confirmedObservations=trainInput.confirmedObservations;
                rotationTrainingInput.weights = weights'; 
                rotationTrainingInput.bias=bias;
                rotationTrainingInput.alpha=trainInput.alpha}
            this.runTrainingRotation(rotationTrainingInput) |> ignore
        {adjustmentOutput.weights=weights'; adjustmentOutput.bias=bias}

        [TestMethod]
        public void trainUsingTestData_RetunsExpected()
        {
            var xValues = new List<double>();
            xValues.Add(1.5);
            xValues.Add(2.0);
            var observation = new observation(xValues);
            var yExpected = -1.0;
            var confirmedObservation0 = new confirmedObservation(observation, yExpected);

            xValues = new List<double>();
            xValues.Add(2.0);
            xValues.Add(3.5);
            observation = new observation(xValues);
            yExpected = -1.0;
            var confirmedObservation1 = new confirmedObservation(observation, yExpected);

            xValues = new List<double>();
            xValues.Add(3.0);
            xValues.Add(5.0);
            observation = new observation(xValues);
            yExpected = -1.0;
            var confirmedObservation2 = new confirmedObservation(observation, yExpected);

            xValues = new List<double>();
            xValues.Add(3.5);
            xValues.Add(2.5);
            observation = new observation(xValues);
            yExpected = -1.0;
            var confirmedObservation3 = new confirmedObservation(observation, yExpected);

            xValues = new List<double>();
            xValues.Add(4.5);
            xValues.Add(5.0);
            observation = new observation(xValues);
            yExpected = 1.0;
            var confirmedObservation4 = new confirmedObservation(observation, yExpected);

            xValues = new List<double>();
            xValues.Add(5.0);
            xValues.Add(7.5);
            observation = new observation(xValues);
            yExpected = 1.0;
            var confirmedObservation5 = new confirmedObservation(observation, yExpected);

            xValues = new List<double>();
            xValues.Add(5.5);
            xValues.Add(8.0);
            observation = new observation(xValues);
            yExpected = 1.0;
            var confirmedObservation6 = new confirmedObservation(observation, yExpected);

            xValues = new List<double>();
            xValues.Add(6.0);
            xValues.Add(6.0);
            observation = new observation(xValues);
            yExpected = 1.0;
            var confirmedObservation7 = new confirmedObservation(observation, yExpected);


            var trainingObservations = new List<confirmedObservation>();
            trainingObservations.Add(confirmedObservation0);
            trainingObservations.Add(confirmedObservation1);
            trainingObservations.Add(confirmedObservation2);
            trainingObservations.Add(confirmedObservation3);
            trainingObservations.Add(confirmedObservation4);
            trainingObservations.Add(confirmedObservation5);
            trainingObservations.Add(confirmedObservation6);
            trainingObservations.Add(confirmedObservation7);

            var random = new Random();
            var weightSeedValue = random.NextDouble();
            var biasSeedValue = random.NextDouble();
            var alpha = .001;
            var maxEpoches = 100;

            var trainInput = new trainInput(trainingObservations, weightSeedValue, biasSeedValue, alpha, maxEpoches);
            var trainOutput = _perceptron.train(trainInput);

            Assert.IsNotNull(trainOutput);

        }

With the training out of the way, I could concentrate on the prediction. The prediction was much easier because there are no adjustments and the rotation is run once. The data structures are also simpler because I don’t have to pass in the knownY values. I also only have 1 covering (all be it long) unit test that looks that the results of the prediction.

    member this.runPredictionCycle (cyclePredictionInput:cyclePredictionInput) = 
        let neuronInput = 
            {neuronInput.weightedObservation=cyclePredictionInput.weightedObservation;
            neuronInput.bias=cyclePredictionInput.bias}
        let neuronResult = this.runNeuron(neuronInput)
        this.runActivation(neuronResult)

    member this.runPredictionRotation (rotationPredictionInput:rotationPredictionInput) = 
        let output = new List<List<float>*float>();
        let weights = rotationPredictionInput.weights
        for i=0 to rotationPredictionInput.observations.Count-1 do
            let currentObservation = rotationPredictionInput.observations.[i];
            let xws = Seq.zip currentObservation.xValues weights
            let xws' = new List<(float*float)>(xws)
            let weightedObservation = {xws=xws'}
            let cyclePredictionInput = {
                            cyclePredictionInput.weightedObservation = weightedObservation;
                            cyclePredictionInput.bias = rotationPredictionInput.bias}
            let cycleOutput = this.runPredictionCycle(cyclePredictionInput)
            output.Add(currentObservation.xValues, cycleOutput)
        output

    member this.predict(predictInput:predictInput) = 
        let rotationPredictionInput = {
                        rotationPredictionInput.observations = predictInput.observations;
                        rotationPredictionInput.weights = predictInput.weights;
                        rotationPredictionInput.bias = predictInput.bias }
        this.runPredictionRotation(rotationPredictionInput)

        [TestMethod]
        public void predictUsingTestData_ReturnsExpected()
        {
            var xValues = new List<double>();
            xValues.Add(3.0);
            xValues.Add(4.0);
            var observation0 = new observation(xValues);

            xValues = new List<double>();
            xValues.Add(0.0);
            xValues.Add(1.0);
            var observation1 = new observation(xValues);

            xValues = new List<double>();
            xValues.Add(2.0);
            xValues.Add(5.0);
            var observation2 = new observation(xValues);

            xValues = new List<double>();
            xValues.Add(5.0);
            xValues.Add(6.0);
            var observation3 = new observation(xValues);

            xValues = new List<double>();
            xValues.Add(9.0);
            xValues.Add(9.0);
            var observation4 = new observation(xValues);

            xValues = new List<double>();
            xValues.Add(4.0);
            xValues.Add(6.0);
            var observation5 = new observation(xValues);

            var observations = new List<observation>();
            observations.Add(observation0);
            observations.Add(observation1);
            observations.Add(observation2);
            observations.Add(observation3);
            observations.Add(observation4);
            observations.Add(observation5);


            var weights = new List<double>();
            weights.Add(.0065);
            weights.Add(.0123);

            var bias = -0.0906;

            var predictInput = new predictInput(observations, weights, bias);

            var predictOutput = _perceptron.predict(predictInput);
            Assert.IsNotNull(predictOutput);
        }

When I run all of the unit tests the all run green:

With the Perceptron created, I can now go back and change the code and figure out:

1) Why my weights across the XValues are the same (wrong!)

2) How to implement a more idomatic/recursive way of running rotations so I can remove the mutation

With my unit tests running green, I know I am covered in case I make a mistake

Filed under F#, Machine Learning

Neural Network Part 2: Perceptrons

July 22, 2014 1 Comment

I started working though the second chapter of McCaffrey’s book Neural Networks Using C# Succinctly to see if I could write the examples using F#.

McCaffrey’s code is tough to read though because of its emphasis on loops and global mutable variables. I read though his description and this is how <I think> the Perceptron should be constructed.

The inputs are a series of independent variables (in this case age and income) and the output is a single dependent variable (in this case party affiliation). The values have been encoded and normalized like in this post here.

An example of the input (from page 31 of his book) is:

Or in a more abstract manner:

In terms of data structures, individual inputs (each row) is placed into an array of floats and the output is a single float

I call this single set of inputs an “observation” (my words, not McCaffrey).

Looking at McCaffrey’s example for a perceptron Input-Output,

all of the variables you need are not included. Here is what you need:

Where A0 and B0 are the same as X0 and X1 respectively in his diagram. Also, McCaffrey uses the word “Perceptron” to mean two different concepts: the entire system as a whole and the individual calculation for a given list of X and Bias. I am a big believer of domain ubiquitous languages so I am calling the individual calculation a neuron.

Once you run these values through the neuron for the 1st observation, you might have to alter the Weights and Bias based on the (Y)result. Therefore, the data structure coming out of the Neuron is

These values are feed into the adjustment function to alter the weights and bias with the output as

I am calling this process of taking the a single observation, the xWeights, , and the bias and turning them into a series of weights and bais as a “cycle” (my words, not McCaffrey)

The output of a cycle is then fed with the next observation and the cycle repeats for as many observations as there are fed into the system.

I am calling the process of running a cycle for each observation in the input dataset a rotation (my words, not McCaffrey) and that the perceptron runs rotations for an x number of times to train itself.

Finally, the Perceptron takes a new set of observations where the Y is not known and runs a Rotation once to predict what the Y will be.

So with that mental image in place, the coding became much easier. Basically, there was a 1 to 1 correspondence of F# functions to each step laid out. I started with an individual cycle

type cycleInput = {xValues:float list;yExpected:float;mutable weights:float list;mutable bias:float;alpha:float}
 
let runNeuron (input:cycleInput) = 
    let valuesAndWeights = input.xValues |> List.zip input.weights
    let output = valuesAndWeights 
                    |> List.map(fun (xValue, xWeight) -> xValue*xWeight) 
                    |> List.sumBy(fun x -> x)
    output + input.bias
 
let runActivation input = 
    if input < 0.0 then -1.0 else 1.0

I used record types all over the place in this code just so I could keep things straight in my head. McCaffrey uses ambiguously-named arrays and global variables. Although this makes my code a bit more wordy (esp for functional people), I think the increased readability is worth the trade-off.

In any event, with the Neuron and Activation calc out of the way, I created the functions that adjust the weights and bias:

let calculateWeightAdjustment(xValue, xWeight, alpha, delta) =
    match delta > 0.0, xValue >= 0.0 with
        | true,true -> xWeight – (alpha * delta * xValue)
        | false,true -> xWeight + (alpha * delta * xValue)
        | true,false -> xWeight – (alpha * delta * xValue)
        | false,false -> xWeight + (alpha * delta * xValue)
 
let calculateBiasAdjustment(bias, alpha, delta) =
    match delta > 0.0 with
        | true -> bias – (alpha * delta)
        | false -> bias + (alpha * delta)

This code is significantly different than the for, nested if that McCaffrey uses.

I maintain using this kind of pattern matching makes the intention much easier to comprehend. I also split out the adjustment of the weights and the adjustment of the bias into individual functions.

With these functions ready, I created an input and output record type and implemented the adjustment function

let runAdjustment (input:adjustmentInput) =
    match input.yExpected = input.yActual with
        | true -> {weights=input.weights;bias=input.bias;yActual=input.yActual}
        | false ->
            let delta = input.yActual – input.yExpected
            let valuesAndWeights = input.xValues |> List.zip input.weights
            let weights' =  valuesAndWeights |> List.map(fun (xValue, xWeight) -> calculateWeightAdjustment(xValue,xWeight,input.alpha,delta))
            let bias' = calculateBiasAdjustment(input.bias,input.alpha,delta)
            {weights=weights';bias=bias';yActual=input.yActual}

There is not a corresponding method in McCaffrey’s code, rather he just does some Array.copy and mutates the global variables in the Update method. I am not a fan of side-effect programming so I created a function that explicitly does the modification.

And to wrap up the individual cycle:

let runCycle (cycleInput:cycleInput) =
    let neuronResult = runNeuron(cycleInput)
    let activationResult = runActivation(neuronResult)
    let adjustmentInput = {xValues=cycleInput.xValues;weights=cycleInput.weights;yExpected=cycleInput.yExpected;
                            bias=cycleInput.bias;alpha=cycleInput.alpha;yActual=activationResult}
    runAdjustment(adjustmentInput)

Up next is to run the cycle for each of the observations (called a rotation)

type observation = {xValues:float list;yExpected:float}
type rotationInput = {observations: observation list;mutable weights:float list;mutable bias:float;alpha:float}
type trainingRotationOutput = {weights:float list; bias:float}
type predictionRotationOutput = {observation: observation;yActual:float}
 
let runTrainingRotation(rotationInput: rotationInput)=
    for i=0 to rotationInput.observations.Length do
        let observation = rotationInput.observations.[i]
        let neuronInput = {cycleInput.xValues=observation.xValues;cycleInput.yExpected=observation.yExpected;cycleInput.weights=rotationInput.weights;
                            cycleInput.bias=rotationInput.bias;cycleInput.alpha=rotationInput.alpha}
        let cycleOutput = runCycle(neuronInput)
        rotationInput.weights <- cycleOutput.weights
        rotationInput.bias <- cycleOutput.bias
    {weights=rotationInput.weights; bias=rotationInput.bias}

Again, note the liberal use of records to keep the inputs and outputs clear. I also created a prediction rotation that is designed to be run only once that does not alter the weights and bias.

let runPredictionRotation(rotationInput: rotationInput)=
    let output = new System.Collections.Generic.List<predictionRotationOutput>()
    for i=0 to rotationInput.observations.Length do
        let observation = rotationInput.observations.[i]
        let neuronInput = {cycleInput.xValues=observation.xValues;cycleInput.yExpected=observation.yExpected;cycleInput.weights=rotationInput.weights;
                            cycleInput.bias=rotationInput.bias;cycleInput.alpha=rotationInput.alpha}
        let cycleOutput = runCycle(neuronInput)
        let predictionRotationOutput = {observation=observation;yActual=cycleOutput.yActual}
        output.Add(predictionRotationOutput)   
    output

With the rotations done, the last step was to create the Perceptron to train and then predict:

type perceptronInput = {observations: observation list;weights:float list;bias:float}
type perceptronOutput = {weights:float list; bias:float}
 
let initializeWeights(xValues, randomSeedValue) =
    let lo = -0.01
    let hi = 0.01
    let xWeight = (hi-lo) * randomSeedValue + lo
    xValues |> List.map(fun w -> xWeight)
 
let initializeBias(randomSeedValue) = 
    let lo = -0.01
    let hi = 0.01
    (hi-lo) * randomSeedValue + lo
 
let runTraining(perceptronInput: perceptronInput, maxEpoches:int) =
    let random = System.Random()
    let alpha = 0.001
    let baseObservation = perceptronInput.observations.[0]
    let mutable weights = initializeWeights(baseObservation.xValues,random.NextDouble())       
    let mutable bias = initializeBias(random.NextDouble())
    let rotationList = [0..maxEpoches]
    for i=0 to maxEpoches do
        let rotationInput = {observations=perceptronInput.observations;weights=weights;bias=bias;alpha=alpha}
        let rotationOutput = runTrainingRotation(rotationInput)
        weights <- rotationOutput.weights
        bias <- rotationOutput.bias
    {weights=weights;bias=bias}
 
let runPrediction(perceptronInput: perceptronInput, weights: float list, bias: float) =
    let random = System.Random()
    let alpha = 0.001
    let rotationInput = {observations=perceptronInput.observations;weights=weights;bias=bias;alpha=alpha}
    runPredictionRotation(rotationInput)

Before I go too much further, I have a big code smell. I am iterating and using the mutable keyword. I am not sure how to take the results of a function that is applied to the 1st element in a sequence and then input that into the second. I need to do that with the weights and bias data structures –> each time it is used in a expression, it need to change and feed into the next expression. I think the answer is the List.Reduce, so I am going to pick this up after looking at that in more detail. I also need to implement the shuffle method so that that cycles are not called in the same order across rotations….

Filed under F#, Machine Learning

Neural Networks

July 15, 2014 1 Comment

I picked up James McCaffrey’s Neural Networks Using C# a couple of weeks ago and decided to see if I could rewrite the code in F#. Unfortunately, the source code is not available (as far as I could tell), so I did some C# then F# coding to see if I could get functional equivalence.

My first stop was chapter one. I made the decision to get the F# code working for the sample data that McCaffrey provided first and then refactor it to a more general program that would work with inputs and values of different datasets. My final upgrade will be use Deedle instead of any other data structure. But first things first, I want to get the examples working so I fired up a script file and opened my REPL.

McCaffrey defines a sample dataset like this

string[] sourceData = new string[] { "Sex Age Locale Income Politics", 
    "==============================================", 
    "Male 25 Rural 63,000.00 Conservative", 
    "Female 36 Suburban 55,000.00 Liberal", "Male 40 Urban 74,000.00 Moderate", 
    "Female 23 Rural 28,000.00 Liberal" };

He then creates a parser for the comma-delimited string values into a double[][]. I just created the dataset as a List of tuples.

let chapter1TestData = [("Male",25.,"Rural",63000.00,"Conservative");
                ("Female",36.,"Suburban",55000.00,"Liberal");
                ("Male",40.,"Urban",74000.00,"Moderate");
                ("Female",23.,"Rural",28000.00,"Liberal")]

I did try an implementation using a record type but for reasons below, I am using Tuples. With the equivalent data loaded into the REPL, I tackled the first supporting function: MinMax. Here is the C# code that McCaffrey wrote:

static void MinMaxNormal(double[][] data, int column)
{ 
    int j = column;
    double min = data[0][j];
    double max = data[0][j]; 
    for (int i = 0; i < data.Length; ++i)
    {
        if (data[i][j] < min) min = data[i][j]; 
        if (data[i][j] > max) max = data[i][j];
    }
    double range = max – min; 
    if (range == 0.0) // ugly 
    { for (int i = 0; i < data.Length; ++i)
        data[i][j] = 0.5;
        return; } 
    for (int i = 0; i < data.Length; ++i)
        data[i][j] = (data[i][j] – min) / range;
}

and here is the equivalent F# code.

let minMax (fullSet, i) =
    let min = fullSet |> Seq.min
    let max = fullSet |> Seq.max
    (i-min)/(max-min)

Note that McCaffrey does not have any unit tests but when I ran the dummy data through the F# implementation, the results matched his screen shots so that will work well enough. If you ever need a reason to use F#, consider those 2 code samples. Granted McCaffrey’s code is more abstract because it can be any column in double array, but my counterpoint is that the function is really doing too much and it is trivial in F# to pick a given column. Is there any doubt what the F# code is doing? Is there any certainty of what the C# code is doing?

In any event, moving along to the next functions, McCaffrey created two functions that do all of the encoding of the string values to appropriate numeric ones. Depending on if the value is a X value (independent) or Y value (dependent), there is a different encoding scheme:

 static string EffectsEncoding(int index, int N)
 {
     // If N = 3 and index = 0 -> 1,0. 
     // If N = 3 and index = 1 -> 0,1. 
     // If N = 3 and index = 2 -> -1,-1. 
     if (N == 2)
     // Special case. 
     { if (index == 0) return "-1"; else if (index == 1) return "1"; }
     int[] values = new int[N – 1];
     if (index == N – 1)
     // Last item is all -1s. 
     { for (int i = 0; i < values.Length; ++i) values[i] = -1; }
     else
     {
         values[index] = 1;
         // 0 values are already there.
     } string s = values[0].ToString();
     for (int i = 1; i < values.Length; ++i) s += "," + values[i]; return s;
 }
 
 static string DummyEncoding(int index, int N)
 { 
     int[] values = new int[N]; values[index] = 1;
     string s = values[0].ToString();
     for (int i = 1; i < values.Length; ++i) s += "," + values[i]; 
     return
}

In my F# project, I decided to domain-specific encoding. I plan to refactor this to something more abstract.

//Transform Sex
let testData' = chapter1TestData |> Seq.map(fun (s,a,l,i,p) -> match s with
                                                               | "Male"-> -1.0,a,l,i,p
                                                             | "Female" -> 1.0,a,l,i,p
                                                             | _ -> failwith "Invalid sex")
//Normalize Age
let testData'' = 
    let fullSet =  testData' |> Seq.map(fun (s,a,l,i,p) -> a)
    testData' |> Seq.map(fun (s,a,l,i,p) -> s,minMax(fullSet,a),l,i,p)
 
//Transform Locale
let testData''' = testData'' |> Seq.map(fun (s,a,l,i,p) -> match l with
                                                                | "Rural" -> s,a,1.,0.,i,p
                                                                | "Suburban" -> s,a,0.,1.,i,p
                                                                | "Urban" -> s,a,-1.,-1.,i,p
                                                                | _ -> failwith "Invalid locale")
//Transform and Normalize Income
let testData'''' = 
    let fullSet =  testData''' |> Seq.map(fun (s,a,l0,l1,i,p) -> i)
    testData''' |> Seq.map(fun (s,a,l0,l1,i,p) -> s,a,l0,l1,minMax(fullSet,i),p)
 
//Transform Politics
let testData''''' = testData'''' |> Seq.map(fun (s,a,l0,l1,i,p) -> match p with
                                                                | "Conservative" -> s,a,l0,l1,i,1.,0.,0.
                                                                | "Liberal" -> s,a,l0,l1,i,0.,1.,0.
                                                                | "Moderate" -> s,a,l0,l1,i,0.,0.,1.
                                                                | _ -> failwith "Invalid politics")

When I execute the script:

Which is the same as McCaffrey’s.

Note that he used Gaussian normalization on column 2 and I did Min/Max based on his advice in the book.

Filed under Analytics, F#

Terminator Program: Part 2

July 8, 2014 1 Comment

Following up on my last post, I decided to send the entire photograph to Sky Biometry and have them parse the photograph and identify individual people. This ability is built right into their API. For example, if you pass them this picture, you get the following json back.

I added the red highlight to show that Sky Biometry can recognize multiple people (it is an array of uids) and that each face tag has a center.x and center:y. Reading the API documentation, this point is center of the face tag point and their point is a percentage of the photo width.

So I need to translate the center point of the skeleton from the Kinect to eqiv center point of the sky biometry recognition output and I should be able to identify individual people within the Kinect’s field of vision. Going back to the Kinect code, I ditched the DrawBoxAroundHead method and altered the UpdateDisplay method like so

private void UpdateDisplay(byte[] colorData, Skeleton[] skeletons)
{
    if (_videoBitmap == null)
    {
        _videoBitmap = new WriteableBitmap(640, 480, 96, 96, PixelFormats.Bgr32, null);
    }
    _videoBitmap.WritePixels(new Int32Rect(0, 0, 640, 480), colorData, 640 * 4, 0);
    kinectColorImage.Source = _videoBitmap;
    var selectedSkeleton = skeletons.FirstOrDefault(s => s.TrackingState == SkeletonTrackingState.Tracked);
    if (selectedSkeleton != null)
    {
        var headPosition = selectedSkeleton.Joints[JointType.Head].Position;
        var adjustedHeadPosition =
            _sensor.CoordinateMapper.MapSkeletonPointToColorPoint(headPosition, ColorImageFormat.RgbResolution640x480Fps30);
        var adjustedSkeletonPosition = _sensor.CoordinateMapper.MapSkeletonPointToColorPoint(selectedSkeleton.Position, ColorImageFormat.RgbResolution640x480Fps30);
 
        skeletonCanvas.Children.Clear();
        Rectangle headRectangle = new Rectangle();
        headRectangle.Fill = new SolidColorBrush(Colors.Blue);
        headRectangle.Width = 10;
        headRectangle.Height = 10;
        Canvas.SetLeft(headRectangle, adjustedHeadPosition.X);
        Canvas.SetTop(headRectangle, adjustedHeadPosition.Y);
        skeletonCanvas.Children.Add(headRectangle);
 
        Rectangle skeletonRectangle = new Rectangle();
        skeletonRectangle.Fill = new SolidColorBrush(Colors.Red);
        skeletonRectangle.Width = 10;
        skeletonRectangle.Height = 10;
        Canvas.SetLeft(skeletonRectangle, adjustedHeadPosition.X);
        Canvas.SetTop(skeletonRectangle, adjustedHeadPosition.Y);
        skeletonCanvas.Children.Add(skeletonRectangle);
 
        String skeletonInfo = headPosition.X.ToString() + " : " + headPosition.Y.ToString() + " — ";
        skeletonInfo = skeletonInfo + adjustedHeadPosition.X.ToString() + " : " + adjustedHeadPosition.Y.ToString() + " — ";
        skeletonInfo = skeletonInfo + adjustedSkeletonPosition.X.ToString() + " : " + adjustedSkeletonPosition.Y.ToString();
 
        skeletonInfoTextBox.Text = skeletonInfo;
 
    }
}

Notice that there are two rectangles because I was not sure if the Head.Position or the Skeleton.Position would match SkyBiometry. Turns out that I want the Head.Position for SkyBiometry (besides, the terminator would want head shots only)

So I ditched the Skeleton.Position. I then needed a way to translate the Head.Posotion.X to SkyBiometry.X and Head.Posotion.Y to SkyBiometry.Y. Fortunately, I know the size of each photograph (640 X 480) so calculating the percent is an exercise of altering UpdateDisplay:

private void UpdateDisplay(byte[] colorData, Skeleton[] skeletons)
{
    Int32 photoWidth = 640;
    Int32 photoHeight = 480;
 
    if (_videoBitmap == null)
    {
        _videoBitmap = new WriteableBitmap(photoWidth, photoHeight, 96, 96, PixelFormats.Bgr32, null);
    }
    _videoBitmap.WritePixels(new Int32Rect(0, 0, photoWidth, photoHeight), colorData, photoWidth * 4, 0);
    kinectColorImage.Source = _videoBitmap;
    var selectedSkeleton = skeletons.FirstOrDefault(s => s.TrackingState == SkeletonTrackingState.Tracked);
    if (selectedSkeleton != null)
    {
        var headPosition = selectedSkeleton.Joints[JointType.Head].Position;
        var adjustedHeadPosition =
            _sensor.CoordinateMapper.MapSkeletonPointToColorPoint(headPosition, ColorImageFormat.RgbResolution640x480Fps30);
 
        skeletonCanvas.Children.Clear();
        Rectangle headRectangle = new Rectangle();
        headRectangle.Fill = new SolidColorBrush(Colors.Blue);
        headRectangle.Width = 10;
        headRectangle.Height = 10;
        Canvas.SetLeft(headRectangle, adjustedHeadPosition.X);
        Canvas.SetTop(headRectangle, adjustedHeadPosition.Y);
        skeletonCanvas.Children.Add(headRectangle);
 
        var skyBiometryX = ((float)adjustedHeadPosition.X / photoWidth)*100;
        var skyBioMetryY = ((float)adjustedHeadPosition.Y / photoHeight)*100;
 
        String skeletonInfo = adjustedHeadPosition.X.ToString() + " : " + adjustedHeadPosition.Y.ToString() + " — ";
        skeletonInfo = skeletonInfo + Math.Round(skyBiometryX,2).ToString() + " : " + Math.Round(skyBioMetryY,2).ToString();
 
        skeletonInfoTextBox.Text = skeletonInfo;
 
    }

And so now I have

The next step is to get the Kinect photo to Sky Biometry. I decided to use Azure Blob Storage as my intermediately location. I updated the architectural diagram like so:

At this point, it made sense to move the project over to F# so I could better concentrate on the work that needs to be done and also getting the important code out of the UI code behind. I fired up a F# project in my solution added a couple different implementations of Storing Photos. To keep things consistent, I created a data structure and an interface:

namespace ChickenSoftware.Terminator.Core
 
open System
 
type public PhotoImage (uniqueId:Guid, imageBytes:byte[]) =
    member this.UniqueId = uniqueId
    member this.ImageBytes = imageBytes
 
type IPhotoImageProvider =
    abstract member InsertPhotoImage : PhotoImage -> unit 
    abstract member DeletePhotoImage : Guid -> unit
    abstract member GetPhotoImage : Guid -> PhotoImage 

My 1st stop was to replicate what Miles did with the Save File Dialog box with a File System Provider. It was very much like a C# implementation:

namespace ChickenSoftware.Terminator.Core
 
open System
open System.IO
open System.Drawing
open System.Drawing.Imaging
 
type LocalFileSystemPhotoImageProvider(folderPath: string) = 
 
    member this.GetPhotoImageUri(uniqueIdentifier: Guid) =
        let fileName = uniqueIdentifier.ToString() + ".jpg"
        Path.Combine(folderPath, fileName)
 
    interface IPhotoImageProvider with 
        member this.InsertPhotoImage(photoImage: PhotoImage) = 
            let fullPath = this.GetPhotoImageUri(photoImage.UniqueId)
            use memoryStream = new MemoryStream(photoImage.ImageBytes)
            let image = Image.FromStream(memoryStream)
            image.Save(fullPath)
 
        member this.DeletePhotoImage(uniqueIdentifier: Guid) = 
            let fullPath = this.GetPhotoImageUri(uniqueIdentifier)
            File.Delete(fullPath)        
 
        member this.GetPhotoImage(uniqueIdentifier: Guid) = 
            let fullPath = this.GetPhotoImageUri(uniqueIdentifier)
            use fileStream = new FileStream(fullPath,FileMode.Open)
            let image = Image.FromStream(fileStream)
            use memoryStream = new MemoryStream()
            image.Save(memoryStream,ImageFormat.Jpeg)
            new PhotoImage(uniqueIdentifier, memoryStream.ToArray())

To call the save method, I altered the SavePhoto method in the C# project to use a MemoryStream and not a FileStream:

private void SavePhoto(byte[] colorData)
{
    var bitmapSource = BitmapSource.Create(640, 480, 96, 96, PixelFormats.Bgr32, null, colorData, 640 * 4);
    JpegBitmapEncoder encoder = new JpegBitmapEncoder();
    encoder.Frames.Add(BitmapFrame.Create(bitmapSource));
    using (MemoryStream memoryStream = new MemoryStream())
    {
        encoder.Save(memoryStream);
        PhotoImage photoImage = new PhotoImage(Guid.NewGuid(), memoryStream.ToArray());
 
        String folderUri = @"C:\Data";
        IPhotoImageProvider provider = new LocalFileSystemPhotoImageProvider(folderUri);
 
        provider.InsertPhotoImage(photoImage);
        memoryStream.Close();
    }
    _isTakingPicture = false;
}

And sure enough, it saves the photo to disk:

One problem that took me 20 minutes to uncover is that if you get your file system path wrong, you get the unhelpful exception:

This has been well-bitched about on stack overflow so I won’t comment further.

With the file system up and running, I turned my attention to Azure. Like the File System provider, it is very close to a C# implementation

namespace ChickenSoftware.Terminator.Core
 
open System
open System.IO
open Microsoft.WindowsAzure.Storage
open Microsoft.WindowsAzure.Storage.Blob
 
type AzureStoragePhotoImageProvider(customerUniqueId: Guid, connectionString: string) = 
 
    member this.GetBlobContainer(blobClient:Blob.CloudBlobClient) =
        let container = blobClient.GetContainerReference(customerUniqueId.ToString())
        if not (container.Exists()) then
            container.CreateIfNotExists() |> ignore
            let permissions = new BlobContainerPermissions()
            permissions.PublicAccess <- BlobContainerPublicAccessType.Blob
            container.SetPermissions(permissions)
        container
 
    member this.GetBlockBlob(uniqueIdentifier: Guid) =
        let storageAccount = CloudStorageAccount.Parse(connectionString)
        let blobClient = storageAccount.CreateCloudBlobClient()
        let container = this.GetBlobContainer(blobClient) 
        let photoUri = this.GetPhotoImageUri(uniqueIdentifier)
        container.GetBlockBlobReference(photoUri) 
 
    member this.GetPhotoImageUri(uniqueIdentifier: Guid) =
        uniqueIdentifier.ToString() + ".jpg"
 
    interface IPhotoImageProvider with 
        member this.InsertPhotoImage(photoImage: PhotoImage) =
            let blockBlob = this.GetBlockBlob(photoImage.UniqueId) 
            use memoryStream = new MemoryStream(photoImage.ImageBytes)
            blockBlob.UploadFromStream(memoryStream)
 
        member this.DeletePhotoImage(uniqueIdentifier: Guid) = 
            let blockBlob = this.GetBlockBlob(uniqueIdentifier) 
            blockBlob.Delete()       
 
        member this.GetPhotoImage(uniqueIdentifier: Guid) =
            let blockBlob = this.GetBlockBlob(uniqueIdentifier) 
            if blockBlob.Exists() then
                blockBlob.FetchAttributes()
                use memoryStream = new MemoryStream()
                blockBlob.DownloadToStream(memoryStream)
                let photoArray = memoryStream.ToArray()
                new PhotoImage(uniqueIdentifier,photoArray)
            else
                failwith "photo not found"

And when I pop it into the WPF application,

private void SavePhoto(byte[] colorData)
{
    var bitmapSource = BitmapSource.Create(640, 480, 96, 96, PixelFormats.Bgr32, null, colorData, 640 * 4);
    JpegBitmapEncoder encoder = new JpegBitmapEncoder();
    encoder.Frames.Add(BitmapFrame.Create(bitmapSource));
    using (MemoryStream memoryStream = new MemoryStream())
    {
        encoder.Save(memoryStream);
        PhotoImage photoImage = new PhotoImage(Guid.NewGuid(), memoryStream.ToArray());
 
        Guid customerUniqueId = new Guid("7282AF48-FB3D-489B-A572-2EFAE80D0A9E");
        String connectionString =
            "DefaultEndpointsProtocol=http;AccountName=XXX;AccountKey=XXX";
        IPhotoImageProvider provider = new AzureStoragePhotoImageProvider(customerUniqueId, connectionString);
 
 
        provider.InsertPhotoImage(photoImage);
        memoryStream.Close();
    }
    _isTakingPicture = false;
}

I can now write my images to Azure.

With that out of the way, I can now have SkyBiometry pick up my photo, analyze it, and push the results back. I went ahead and added in the .fs module that I had already created for this blog post. I then added FSharp.Data via NuGet and was ready to roll. In he Save photo event handler,after saving the photo to blob storage, it then calls Sky Biometry to compare against a base image that has already been trained:

private void SavePhoto(byte[] colorData)
{
    var bitmapSource = BitmapSource.Create(640, 480, 96, 96, PixelFormats.Bgr32, null, colorData, 640 * 4);
    JpegBitmapEncoder encoder = new JpegBitmapEncoder();
    encoder.Frames.Add(BitmapFrame.Create(bitmapSource));
    PhotoImage photoImage = UploadPhotoImage(encoder);
 
    String skyBiometryUri = "http://api.skybiometry.com&quot;;
    String uid = "Kinect@ChickenFace";
    String apiKey = "XXXX";
    String apiSecret = "XXXX";
 
    var imageComparer = new SkyBiometryImageComparer(skyBiometryUri, uid, apiKey, apiSecret);
    String basePhotoUri = "XXXX.jpg";
    String targetPhotoUri = "XXXX/" + photoImage.UniqueId + ".jpg";
 
    currentImage.Source = new BitmapImage(new Uri(basePhotoUri));
    compareImage.Source = new BitmapImage(new Uri(targetPhotoUri)); ;
    
    var matchValue = imageComparer.CalculateFacialRecognitionConfidence(basePhotoUri, targetPhotoUri);
    FacialRecognitionTextBox.Text = "Match Value is: " + matchValue.ToString();
    _isTakingPicture = false;
}

And I am getting a result back from Sky Biometry.

Finally, I added in the SkyBiometry X and Y coordinates for the photo and compared to the calculated ones based on the Kinect Skeleton Tracking:

currentImage.Source = new BitmapImage(new Uri(basePhotoUri));
compareImage.Source = new BitmapImage(new Uri(targetPhotoUri)); ;
 
var matchValue = imageComparer.CalculateFacialRecognitionConfidence(basePhotoUri, targetPhotoUri);
 
var selectedSkeleton = skeletons.FirstOrDefault(s => s.TrackingState == SkeletonTrackingState.Tracked);
if (selectedSkeleton != null)
{
    var headPosition = selectedSkeleton.Joints[JointType.Head].Position;
    var adjustedHeadPosition =
        _sensor.CoordinateMapper.MapSkeletonPointToColorPoint(headPosition, ColorImageFormat.RgbResolution640x480Fps30);
 
    var skyBiometryX = ((float)adjustedHeadPosition.X / 640) * 100;
    var skyBioMetryY = ((float)adjustedHeadPosition.Y / 480) * 100;
 
    StringBuilder stringBuilder = new StringBuilder();
    stringBuilder.Append("Match Value is: ");
    stringBuilder.Append(matchValue.Confidence.ToString());
    stringBuilder.Append("Sky Biometry X: ");
    stringBuilder.Append(matchValue.X.ToString());
    stringBuilder.Append("Sky Biometry Y: ");
    stringBuilder.Append(matchValue.Y.ToString());
    stringBuilder.Append("Kinect X: ");
    stringBuilder.Append(Math.Round(skyBiometryX, 2).ToString());
    stringBuilder.Append("Kinect Y: ");
    stringBuilder.Append(Math.Round(skyBioMetryY, 2).ToString());
    FacialRecognitionTextBox.Text = stringBuilder.ToString();
}
 
_isTakingPicture = false;

And the results are encouraging –> it looks like I can use the X and Y to identify different people on the screen:

Match Value is: 53
Sky Biometry X: 10
Sky Biometry Y: 13.33

Kinect X: 47.5
Kinect Y: 39.79

Up next will be pointing the laser and the target…

Filed under C#, F#, Kinect, WPF

TRINUG F# Analytics Prep: Part 2

July 1, 2014 1 Comment

I finished up the second part of the F#/Analytics lab scheduled for August. It is a continuation of going through Astborg’s F# for Quantitative Finance that we started last month. Here is my fist blog post on it.

In this lab, we are going to tackle the more advanced statistical calculations: the Black-Scholes formula, the Greeks, and Monte Carlo simulation. Using the same solution and projects, I started the script file to figure out the Black Scholes formula. Astborg uses a couple of supporting functions which I knocked out first: Power and CumulativeDistribution. I first created his function verbatim like this:

let pow x n = exp(n*log(x))

and then refactored it to make it more readable like this

let power baseNumber exponent = exp(exponent * log(baseNumber))

and then I realized it is the same thing as using pown which is already found in FSharp.Core.

In any event, I then attacked the cumulativeDistribution method. I downloaded the source from his website and then refactored it so that each step is clearly laid out. Here is the refactored function:

let cumulativeDistribution (x) =
        let a1 =  0.31938153
        let a2 = -0.356563782
        let a3 =  1.781477937
        let a4 = -1.821255978
        let a5 =  1.330274429
        let pi = 3.141592654
        let l  = abs(x)
        let k  = 1.0 / (1.0 + 0.2316419 * l)
 
        let a1' = a1*k
        let a2' = a2*k*k
        let a3' = a3*(power k 3.0)
        let a4' = a4*(power k 4.0)
        let a5' = a5*(power k 5.0)
        let w1 = 1.0/sqrt(2.0*pi)
        let w2 = exp(-l*l/2.0)
        let w3 = a1'+a2'+a3'+a4'+a5'
        let w  = 1.0-w1*w2*w3
        if x < 0.0 then 1.0 – w else w

And here is some test values from the REPL:

Finally, the Black Scholes formula. I did create a separate POCO for the input data like this:

type putCallFlag = Put | Call
 
type blackScholesInputData = 
    {stockPrice:float;
    strikePrice:float;
    timeToExpiry:float;
    interestRate:float;
    volatility:float}

And I refactored his code to make it more readable like this:

let blackScholes(inputData:blackScholesInputData, putCallFlag:putCallFlag)=
   let sx = log(inputData.stockPrice / inputData.strikePrice)
   let rv = inputData.interestRate+inputData.volatility*inputData.volatility*0.5
   let rvt = rv*inputData.timeToExpiry
   let vt = (inputData.volatility*sqrt(inputData.timeToExpiry))
   let d1=(sx + rvt)/vt
   let d2=d1-vt
    
   match putCallFlag with
    | Put -> 
        let xrt = inputData.strikePrice*exp(-inputData.interestRate*inputData.timeToExpiry)
        let cdD1 = xrt*cumulativeDistribution(-d2)
        let cdD2 = inputData.stockPrice*cumulativeDistribution(-d1)
        cdD1-cdD2
    | Call ->
        let xrt = inputData.strikePrice*exp(-inputData.interestRate*inputData.timeToExpiry)
        let cdD1 = inputData.stockPrice*cumulativeDistribution(d1)
        let cdD2 = xrt*cumulativeDistribution(d2)
        cdD1-cdD2

And since I was in the script environment, I put in test data that matches the sample that Astborg used in the book:

let inputData = {stockPrice=58.60;strikePrice=60.;timeToExpiry=0.5;interestRate=0.01;volatility=0.3}
let runBSCall = blackScholes(inputData,Call)
let runBSPut = blackScholes(inputData,Put)

And voila, the results match the book:

With the Black-Scholes out of the way, I then implemented the Greeks. Note that I did add helper functions for clarity, and the results match the book:

let blackScholesDelta (inputData:blackScholesInputData, putCallFlag:putCallFlag) =
    let sx = log(inputData.stockPrice / inputData.strikePrice)
    let rv = inputData.interestRate+inputData.volatility*inputData.volatility*0.5
    let rvt = rv*inputData.timeToExpiry
    let vt = (inputData.volatility*sqrt(inputData.timeToExpiry))
    let d1=(sx + rvt)/vt
    match putCallFlag with
    | Put -> cumulativeDistribution(d1) – 1.0
    | Call -> cumulativeDistribution(d1)
 
let deltaPut = blackScholesDelta(inputData, Put)
let deltaCall = blackScholesDelta(inputData, Call)
 
let blackScholesGamma (inputData:blackScholesInputData) =
    let sx = log(inputData.stockPrice / inputData.strikePrice)
    let rv = inputData.interestRate+inputData.volatility*inputData.volatility*0.5
    let rvt = rv*inputData.timeToExpiry
    let vt = (inputData.volatility*sqrt(inputData.timeToExpiry))
    let d1=(sx + rvt)/vt
    normalDistribution.Density(d1)
 
let gamma = blackScholesGamma(inputData)
 
let blackScholesVega (inputData:blackScholesInputData) =
    let sx = log(inputData.stockPrice / inputData.strikePrice)
    let rv = inputData.interestRate+inputData.volatility*inputData.volatility*0.5
    let rvt = rv*inputData.timeToExpiry
    let vt = (inputData.volatility*sqrt(inputData.timeToExpiry))
    let d1=(sx + rvt)/vt   
    inputData.stockPrice*normalDistribution.Density(d1)*sqrt(inputData.timeToExpiry)
 
let vega = blackScholesVega(inputData)
 
let blackScholesTheta (inputData:blackScholesInputData, putCallFlag:putCallFlag) =
    let sx = log(inputData.stockPrice / inputData.strikePrice)
    let rv = inputData.interestRate+inputData.volatility*inputData.volatility*0.5
    let rvt = rv*inputData.timeToExpiry
    let vt = (inputData.volatility*sqrt(inputData.timeToExpiry))
    let d1=(sx + rvt)/vt   
    let d2=d1-vt
    match putCallFlag with
    | Put -> 
        let ndD1 = inputData.stockPrice*normalDistribution.Density(d1)*inputData.volatility
        let ndD1' = ndD1/(2.0*sqrt(inputData.timeToExpiry))
        let rx = inputData.interestRate*inputData.strikePrice
        let rt = exp(-inputData.interestRate*inputData.timeToExpiry)
        let cdD2 = rx*rt*cumulativeDistribution(-d2) 
        -(ndD1')+cdD2
    | Call -> 
        let ndD1 = inputData.stockPrice*normalDistribution.Density(d1)*inputData.volatility
        let ndD1' = ndD1/(2.0*sqrt(inputData.timeToExpiry))
        let rx = inputData.interestRate*inputData.strikePrice
        let rt = exp(-inputData.interestRate*inputData.timeToExpiry)
        let cdD2 = cumulativeDistribution(d2)
        -(ndD1')-rx*rt*cdD2
 
let thetaPut = blackScholesTheta(inputData, Put)
let thetaCall = blackScholesTheta(inputData, Call)
 
let blackScholesRho (inputData:blackScholesInputData, putCallFlag:putCallFlag) =
    let sx = log(inputData.stockPrice / inputData.strikePrice)
    let rv = inputData.interestRate+inputData.volatility*inputData.volatility*0.5
    let rvt = rv*inputData.timeToExpiry
    let vt = (inputData.volatility*sqrt(inputData.timeToExpiry))
    let d1=(sx + rvt)/vt   
    let d2=d1-vt
    match putCallFlag with
    | Put ->
        let xt = inputData.strikePrice*inputData.timeToExpiry
        let rt = exp(-inputData.interestRate*inputData.timeToExpiry)  
        -xt*rt*cumulativeDistribution(-d2)
    | Call -> 
        let xt = inputData.strikePrice*inputData.timeToExpiry
        let rt = exp(-inputData.interestRate*inputData.timeToExpiry)          
        xt*rt*cumulativeDistribution(d2)
 
let rhoPut = blackScholesRho(inputData, Put)
let rhoCall = blackScholesRho(inputData, Call)

Finally, I threw in the Monte Carlo, which also used a POCO:

type monteCarloInputData = 
    {stockPrice:float;
    strikePrice:float;
    timeToExpiry:float;
    interestRate:float;
    volatility:float}
 
let priceAtMaturity (inputData:monteCarloInputData, randomValue:float) =
    let s = inputData.stockPrice
    let rv = (inputData.interestRate-inputData.volatility*inputData.volatility/2.0)
    let rvt = rv*inputData.timeToExpiry
    let vr = inputData.volatility*randomValue
    let t = sqrt(inputData.timeToExpiry)
    s*exp(rvt+vr*t)
    
let maturityPriceInputData = {stockPrice=58.60;strikePrice=60.0;timeToExpiry=0.5;interestRate=0.01;volatility=0.3}
priceAtMaturity(maturityPriceInputData, 10.0)
 
let monteCarlo(inputData: monteCarloInputData, randomValues:seq<float>) = 
    randomValues 
        |> Seq.map(fun randomValue -> priceAtMaturity(inputData,randomValue) – inputData.strikePrice )
        |> Seq.average
 
 
let random = new System.Random()
let rnd() = random.NextDouble()
let data = [for i in 1 .. 1000 -> rnd() * 1.0]
 
let monteCarloInputData = {stockPrice=58.60;strikePrice=60.0;timeToExpiry=0.5;interestRate=0.01;volatility=0.3;}
monteCarlo(monteCarloInputData,data)

One thing I really like about Astborg is that the Monte Carlo function does not new up the array of random numbers, rather they are passed in. This makes the function much more testable and is the right way to right it (IMHO). In fact, I think that seeing “new Random” or “DateTime.Now” hard-coded into functions is an anti-pattern that is all too common.

With the last of the functions done in the script file, I moved them into the .fs file and created covering unit tests based on the sample data that I did in the REPL.

[TestMethod]
public void PowerUsingValidData_ReturnsExpected()
{
    var calculations = new Calculations();
    Double expected = 8;
    Double actual = Math.Round(calculations.Power(2.0, 3.0), 0);
    Assert.AreEqual(expected, actual);
}
 
[TestMethod]
public void CumulativeDistributionUsingValidData_ReturnsExpected()
{
    var calculations = new Calculations();
    Double expected = .84134;
    Double actual = Math.Round(calculations.CumulativeDistribution(1.0),5);
    Assert.AreEqual(expected, actual);
}
 
[TestMethod]
public void BlackScholesCallUsingValidData_ReturnsExpected()
{
    var calculations = new Calculations();
    Double expected = 4.4652;
    var inputData = new BlackScholesInputData(58.6, 60.0, .5, .01, .3);
    Double actual = Math.Round(calculations.BlackScholes(inputData,PutCallFlag.Call), 5);
    Assert.AreEqual(expected, actual);
}
 
[TestMethod]
public void BlackScholesPutUsingValidData_ReturnsExpected()
{
    var calculations = new Calculations();
    Double expected = 5.56595;
    var inputData = new BlackScholesInputData(58.6, 60.0, .5, .01, .3);
    Double actual = Math.Round(calculations.BlackScholes(inputData, PutCallFlag.Put), 5);
    Assert.AreEqual(expected, actual);
}
 
[TestMethod]
public void DaysToYearsUsingValidData_ReturnsExpected()
{
    var calculations = new Calculations();
    Double expected = .08214;
    Double actual = Math.Round(calculations.DaysToYears(30), 5);
    Assert.AreEqual(expected, actual);
}
 
[TestMethod]
public void BlackScholesDeltaCallUsingValidData_ReturnsExpected()
{
    var calculations = new Calculations();
    Double expected = .50732;
    var inputData = new BlackScholesInputData(58.6, 60.0, .5, .01, .3);
    Double actual = Math.Round(calculations.BlackScholesDelta(inputData, PutCallFlag.Call), 5);
    Assert.AreEqual(expected, actual);
}
 
[TestMethod]
public void BlackScholesDeltaPutUsingValidData_ReturnsExpected()
{
    var calculations = new Calculations();
    Double expected = -.49268;
    var inputData = new BlackScholesInputData(58.6, 60.0, .5, .01, .3);
    Double actual = Math.Round(calculations.BlackScholesDelta(inputData, PutCallFlag.Put), 5);
    Assert.AreEqual(expected, actual);
}
 
[TestMethod]
public void BlackScholesGammaUsingValidData_ReturnsExpected()
{
    var calculations = new Calculations();
    Double expected = .39888;
    var inputData = new BlackScholesInputData(58.6, 60.0, .5, .01, .3);
    Double actual = Math.Round(calculations.BlackScholesGamma(inputData), 5);
    Assert.AreEqual(expected, actual);
}
 
[TestMethod]
public void BlackScholesVegaUsingValidData_ReturnsExpected()
{
    var calculations = new Calculations();
    Double expected = 16.52798;
    var inputData = new BlackScholesInputData(58.6, 60.0, .5, .01, .3);
    Double actual = Math.Round(calculations.BlackScholesVega(inputData), 5);
    Assert.AreEqual(expected, actual);
}
 
[TestMethod]
public void BlackScholesThetaCallUsingValidData_ReturnsExpected()
{
    var calculations = new Calculations();
    Double expected = -5.21103;
    var inputData = new BlackScholesInputData(58.6, 60.0, .5, .01, .3);
    Double actual = Math.Round(calculations.BlackScholesTheta(inputData, PutCallFlag.Call), 5);
    Assert.AreEqual(expected, actual);
}
 
[TestMethod]
public void BlackScholesThetaPutUsingValidData_ReturnsExpected()
{
    var calculations = new Calculations();
    Double expected = -4.61402;
    var inputData = new BlackScholesInputData(58.6, 60.0, .5, .01, .3);
    Double actual = Math.Round(calculations.BlackScholesTheta(inputData, PutCallFlag.Put), 5);
    Assert.AreEqual(expected, actual);
}
 
[TestMethod]
public void BlackScholesRhoCallUsingValidData_ReturnsExpected()
{
    var calculations = new Calculations();
    Double expected = 12.63174;
    var inputData = new BlackScholesInputData(58.6, 60.0, .5, .01, .3);
    Double actual = Math.Round(calculations.BlackScholesRho(inputData, PutCallFlag.Call), 5);
    Assert.AreEqual(expected, actual);
}
 
[TestMethod]
public void BlackScholesRhoPutUsingValidData_ReturnsExpected()
{
    var calculations = new Calculations();
    Double expected = -17.21863;
    var inputData = new BlackScholesInputData(58.6, 60.0, .5, .01, .3);
    Double actual = Math.Round(calculations.BlackScholesRho(inputData, PutCallFlag.Put), 5);
    Assert.AreEqual(expected, actual);
}
 
 
[TestMethod]
public void PriceAtMaturityUsingValidData_ReturnsExpected()
{
    var calculations = new Calculations();
    Double expected = 480.36923;
    var inputData = new MonteCarloInputData(58.6, 60.0, .5, .01, .3);
    Double actual = Math.Round(calculations.PriceAtMaturity(inputData, 10.0), 5);
    Assert.AreEqual(expected, actual);
}
 
[TestMethod]
public void MonteCarloUsingValidData_ReturnsExpected()
{
    var calculations = new Calculations();
    var inputData = new MonteCarloInputData(58.6, 60.0, .5, .01, .3);
    var random = new System.Random();
    List<Double> randomData = new List<double>();
    for (int i = 0; i < 1000; i++)
    {
        randomData.Add(random.NextDouble());
    }
 
    Double actual = Math.Round(calculations.MonteCarlo(inputData, randomData), 5);
    var greaterThanFour = actual > 4.0;
    var lessThanFive = actual < 5.0;
 
    Assert.AreEqual(true, greaterThanFour);
    Assert.AreEqual(true, lessThanFive);
}

With all of the tests running green, I then turned my attention to the UI. I created more real state on the MainWindow and added some additional data structures to the results of the analytics that lend themselves to charting and graphing. For example:

public class GreekData
{
    public Double StrikePrice { get; set; }
    public Double DeltaCall { get; set; }
    public Double DeltaPut { get; set; }
    public Double Gamma { get; set; }
    public Double Vega { get; set; }
    public Double ThetaCall { get; set; }
    public Double ThetaPut { get; set; }
    public Double RhoCall { get; set; }
    public Double RhoPut { get; set; }
 
}

And in the code behind of the MainWindow, I added some calcs based on the prior code that was already in it:

var theGreeks = new List<GreekData>();
for (int i = 0; i < 5; i++)
{
    var greekData = new GreekData();
    greekData.StrikePrice = closestDollar – i;
    theGreeks.Add(greekData);
    greekData = new GreekData();
    greekData.StrikePrice = closestDollar + i;
    theGreeks.Add(greekData);
}
theGreeks.Sort((greek1,greek2)=>greek1.StrikePrice.CompareTo(greek2.StrikePrice));
 
foreach (var greekData in theGreeks)
{
    var inputData =
        new BlackScholesInputData(adjustedClose, greekData.StrikePrice, .5, .01, .3);
    greekData.DeltaCall = calculations.BlackScholesDelta(inputData, PutCallFlag.Call);
    greekData.DeltaPut = calculations.BlackScholesDelta(inputData, PutCallFlag.Put);
    greekData.Gamma = calculations.BlackScholesGamma(inputData);
    greekData.RhoCall = calculations.BlackScholesRho(inputData, PutCallFlag.Call);
    greekData.RhoPut = calculations.BlackScholesRho(inputData, PutCallFlag.Put);
    greekData.ThetaCall = calculations.BlackScholesTheta(inputData, PutCallFlag.Call);
    greekData.ThetaPut = calculations.BlackScholesTheta(inputData, PutCallFlag.Put);
    greekData.Vega = calculations.BlackScholesVega(inputData);
 
}
 
this.TheGreeksDataGrid.ItemsSource = theGreeks;
 
 
var blackScholes = new List<BlackScholesData>();
for (int i = 0; i < 5; i++)
{
    var blackScholesData = new BlackScholesData();
    blackScholesData.StrikePrice = closestDollar – i;
    blackScholes.Add(blackScholesData);
    blackScholesData = new BlackScholesData();
    blackScholesData.StrikePrice = closestDollar + i;
    blackScholes.Add(blackScholesData);
}
blackScholes.Sort((bsmc1, bsmc2) => bsmc1.StrikePrice.CompareTo(bsmc2.StrikePrice));
 
var random = new System.Random();
List<Double> randomData = new List<double>();
for (int i = 0; i < 1000; i++)
{
    randomData.Add(random.NextDouble());
}
 
foreach (var blackScholesMonteCarlo in blackScholes)
{
    var blackScholesInputData =
        new BlackScholesInputData(adjustedClose, blackScholesMonteCarlo.StrikePrice, .5, .01, .3);
    var monteCarloInputData =
        new MonteCarloInputData(adjustedClose, blackScholesMonteCarlo.StrikePrice, .5, .01, .3);
 
    blackScholesMonteCarlo.Call = calculations.BlackScholes(blackScholesInputData, PutCallFlag.Call);
    blackScholesMonteCarlo.Put = calculations.BlackScholes(blackScholesInputData, PutCallFlag.Put);
    blackScholesMonteCarlo.MonteCarlo = calculations.MonteCarlo(monteCarloInputData, randomData);
}
 
this.BlackScholesDataGrid.ItemsSource = blackScholes;

And Whammo, the UI.

Fortunately Conrad D’Cruz is a member of TRINUG and an options trader and is going to explain what the heck we are looking at when the SIG gets together again.

Filed under Analytics, F#

← Older posts

Newer posts →

Jamie Dixon's Home

Consuming Azure ML With F#

Sql Saturday and MVP Monday

Fun with Statistics and Charts

Consuming and Analyzing Census Data Using F#

Controlling Servos Using Netdunio and Phidgets

Neural Network Part 3: Perceptrons

Neural Network Part 2: Perceptrons

Neural Networks

Terminator Program: Part 2

TRINUG F# Analytics Prep: Part 2

Categories

Recent Posts

Archives

Blogroll

Meta