Test of R markdown: Code from Feb 26th.
========================================================

This is an R Markdown document. Markdown is a simple formatting syntax for authoring web pages (click the **MD** toolbar button for help on Markdown).

When you click the **Knit HTML** button a web page will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

Gary code
-------------------------------------------------------
I think I lost the original code that produced the plots on my research blot for Feb 3, 2014. 
So I just reconstructed the jist of it, here.

```{r Garychunk}
source("~/Dropbox/Duke/research/code/fn.VM.EM.R")
library(mixtools)
library(movMF)
library(copula); library(VineCopula)
source('~/Dropbox/Duke/research/code/getdata.R')
source('~/Dropbox/Duke/research/code/20131216-functions.R')

vmmarg = 0   # choose which margins to use (von Mises are much better than Gaussian margins)
normmarg = 1

mydata = A.thin
copfam = 30    # 1=Gaussian

Fhat = pobs(mydata)    # pobs comes from 'copula' package
chatpars = BiCopEst(u1=Fhat[,1], u2=Fhat[,2], family=copfam)  # BiCopEst comes from 'VineCopula' package

# Estimate f1hat and f2hat
#      The normalmixEM() function comes from 'mixtools' package
k=3

if (k==2){
     mustart1 = summary(mydata[,1])[c(1,5)]
     mustart2 = summary(mydata[,2])[c(1,5)]
}
if (k==3){
     mustart1 = summary(mydata[,1])[c(1,3,5)]
     mustart2 = summary(mydata[,2])[c(1,3,5)]
}
f1hatpars <-  normalmixEM(x=mydata[,1], k=k, mu=mustart1)
f2hatpars <-  normalmixEM(x=mydata[,2], k=k, mu=mustart2)
f1hatpars.vm <- vm.em(marg=mydata[,1], k=2)
f2hatpars.vm <- vm.em(marg=mydata[,2], k=2)


xgrid = seq(-pi,pi, length.out=100)
d1 = f1hatpars$lambda[1]*dnorm(xgrid, mean=f1hatpars$mu[1], sd=f1hatpars$sigma[1]) 
d2 = f1hatpars$lambda[2]*dnorm(xgrid, mean=f1hatpars$mu[2], sd=f1hatpars$sigma[2]) 
d3 = f1hatpars$lambda[3]*dnorm(xgrid, mean=f1hatpars$mu[3], sd=f1hatpars$sigma[3]) 

hist(mydata[,1],50, freq=FALSE)

lines(xgrid,d1+d2+d3)



d1 = f2hatpars$lambda[1]*dnorm(xgrid, mean=f2hatpars$mu[1], sd=f2hatpars$sigma[1]) 
d2 = f2hatpars$lambda[2]*dnorm(xgrid, mean=f2hatpars$mu[2], sd=f2hatpars$sigma[2]) 
d3 = f2hatpars$lambda[3]*dnorm(xgrid, mean=f2hatpars$mu[3], sd=f2hatpars$sigma[3]) 

hist(mydata[,2],50, freq=FALSE)
#lines(xgrid, d1+d2)
lines(xgrid, d1+d2+d3)
```

