source("cgamma.R")
dpearson4 <- function (x, alpha, m, nu, lambda){
abs(cgamma(m+nu/2*1i)/gamma(m))^2/(alpha*beta(m-1/2, 1/2))*(1+((x-lambda)/alpha)^2)^(-m)*exp(-nu*atan((x-lambda)/alpha))
}

params <- function(alpha, m, nu, lambda){
r <- 2*(m-1);
ave <- lambda - alpha*nu/r;
sig2 <- alpha^2*(r^2+nu^2)/(r^2*(r-1))
skew <- -4*nu/(r-2)*sqrt((r-1)/(r^2+nu^2))
kurt <- 3*(r-1)*((r+6)*(r^2+nu^2)-8*r^2)/((r-2)*(r-3)*(r^2+nu^2)) - 3
c(ave, sig2, skew, kurt)
}

norm2 <- function(x){(t(x)%*%x)[1]}

findparams <- function(m, s2, sk, k){
  f <- function(x){
    norm2(params(x[1],x[2],x[3],x[4])-c(m, s2, sk, k))
  }
  res <- nlm(f, c(sqrt(2), 4, 2, 0));#c(1, 5.3, 0, 1));
  res$estimate
}

t=seq(-5, 5, 0.01);

p <- findparams(0, 1, 0, 0.05)
print(p);

alpha <- p[1]; m <- p[2]; nu <- p[3]; lambda <- p[4];
#print(params(alpha, m, nu, lambda))

plotpearson<-function(alpha, m, nu, lambda){
  plot(t, dnorm(t), type="l", col="red", lwd=3, xlim=c(-5, 5), ylim=c(0, 0.65));
  lines(t, dpearson4(t, alpha, m, nu, lambda), type="l", col="blue", lwd=3)
}


#lambda<-0; nu<-6; alpha<-sqrt(nu); m<-(nu+1)/2;
#print(params(alpha, m, nu, lambda))
#plot(t, dpearson4(t, alpha, m, nu, lambda), type="l")